Tuesday, December 15, 2009

Band Pass Filtering: SCORE!

In the last post, I showed how a simple-minded effort at band-pass filtering failed miserably. Basically the attempt was to apply a simple low-pass exponential weighted rolling average as a low-pass filter, then follow up with differentiation, which cuts out the low-frequency stuff.

The problem was the differentiation basically undid the low-pass filtering of the exponential averaging. What was left was a high pass filter, not a band pass.

Of course I could look this stuff up on Wikipedia or even (gasp!) an actual book (remember those?). But I that would hardly be redemption. Better to dive in and convince myself I still remember something from that expensive college education..

So to retain both high and low-pass characteristics after differentiation what's needed is a better low-pass: it needs to cut off high frequencies faster than the 1/ω behavior of the exponential smoothing.

I previously proposed a truncated sine wave Let's try that (or a truncated cosine wave: same deal). Instead of convolving with an exponential decay, then, we'll convolve with the positive half of a cosine function, with weighting calculated so the integral of the convolution function equals one (so it gives back the original function if that function is sufficiently slowly varying):

f*(t) = (1 / 2τ) ∫ d Δt · cos(Δt / τ) f(t ‒ Δt),

where here the limits of integration are from ‒π/2 τ to +π/2 τ.

Okay, one caveat. This filter is using results of f from future time (Δt < 0). But causality can be trivially restored by inserting a delay of π/2 τ. I'll do the analyis pre-delay to keep the math simpler.

I can then pull the previous trick of representing the cos function in terms of its equivalent complex function:

cos(Δt / τ) = [exp(i Δt / τ) + exp(‒i Δt / τ)] / 2

Then I model for a function (in the last post I used since, but the difference is only in the phase):

f(t) = f0 cos ωt,

for which I will use the complex representation:

f(t) = f0 Re{ exp iωt },

where Re{} denotes the real part of the complex number. .

Okay, simple enough:

f*(t) = (1 / 2τ) ∫ d Δt · cos(Δt / τ) f(t ‒ Δt),

= (f0 / 2τ) Re { ∫ d Δt · [exp(i Δt / τ) + exp(‒i Δt / τ)] exp iωt / 2},

where I recognized the "Re" can be taken outside the integral (the integral of the real part equals the real part of the integral). Consolate the exponentials for integration:

f*(t) = (f0 / 2)) Re { ∫ d Δt · [exp(i Δt / τ + iω (t ‒ Δt) ) + exp(‒i Δt / 2τ + iω(t ‒ Δt))] / 2}.

It's only a bookkeeping challenge to integrate this... (warning: I'm quite poor at bookkeeping):

f*(t) = (f0 / 2) Re { [ exp(i Δt / τ + iω (t ‒ Δt)) / (1 ‒ ωτ) ‒ exp(‒i Δt / τ + iω (t ‒ Δt)) / (1 + ωτ) ] / 2i}.

Consolidate the denominator:

f*(t) = (f0 / 2) Re { [ (1 + ωτ) exp(i Δt / τ + iω (t ‒ Δt)) ‒ (1 ‒ ωτ) exp(‒i Δt / τ + iω (t ‒ Δt)) ] / 2i(1 ‒ ω²τ²)},

which can be brought back to the trigometric form:

f*(t) = (f0 / 2) Re { (sin(Δt / τ) ‒ iωτ cos(Δt / τ)) exp(iω (t ‒ Δt)) / (1 ‒ ω²τ²)}.

Now we need to remember to evaluate this at the limits Δt = ‒π/2 τ and +π/2 τ (we want a definite integral, not the indefinite one):

f*(t) = (f0 / 2) Re { (sin(π/2) ‒ iωτ cos(π/2)) exp(iω (t ‒ π/2 τ)) ‒ (sin(‒π/2) ‒ iωτ cos(‒π/2)) exp(iω (t + π/2 τ)) } / (1 ‒ ω²τ²),

which becomes:

f*(t) = (f0 / 2) Re { exp(iω (t ‒ π/2 τ) + exp(iω (t + π/2 τ)) } / (1 ‒ ω²τ²).

This can e cleaned up a but to to get my final result:

f*(t) = f0 cos(π/2 ωτ) cos(ωt) / (1 ‒ ω²τ²).

For ωτ « 1, this approaches:

f*(t) ≈ f0 cos ωt,

which is huge relief because it what we wanted: for slowly varying functions the smoothing should just yield the original function.

On the other hand, for ωτ » 1, this becomes

f*(t) = ‒f0 cos(π/2 ωτ) cos(ωt) / ω²τ²,

which falls off proportionally to ω²τ². Recall the exponential filter only fell off proportional to ωτ. This is more great news: it means after we differentiate the net result will fall off at high frequency proportional to ωτ, so we have a band pass filter, which is what we wanted. The exponential filter failed to attenuate high frequencies after differentiation.

There is the funny business near ωτ: both the denominator and the numerator go to zero. To see what happens, we can use a polynomial expansion of cosine:

cos (π/2 + ε) ≈ ε,

where ε ≡ (π/2) (ωτ ‒ 1), or equivalently ωτ ≡ 2ε/π + 1. The denominator term (1 ‒ ω²τ²) then becomes 1 ‒ (2ε/π + 1)² ≈ 4ε/π. So the ratio cos((π/2) ωτ) } / (1 ‒ ω²τ²) approaches π/4. So no singularity there after all, despite initial appearances.

Okay, for the record, I suppose it would be good to apply the differentiation explicitly, then multiply by τ to keep the units the same:

τ ∂/∂t f*(t) = ‒f0 cos(π/2 ωτ) sin(ωt) ωτ / (1 ‒ ω²τ²),

which, as promised, falls off both for small and large ωτ.

So what to conclude? Convolution with a truncated cosine wave does a better job of low-pass filtering than an exponentially-weighted rolling average. thats sort of obvious by inspection: there's no hard edges on the truncated cosine wave. But there is a hard corner, so in the future I'll consider an alternative.