A one-step band-pass
So far, I've been doing band pass filtering by a low pass followed by differentiation. This involves first a convolution to realize the smoothing, then a difference function for the differentiation. But now that I think about it, it's better to avoid the difference function and accomplish the full band pass filtering with a single convolution.
Okay, you can stop laughing now. I know this is all super-basic-DSP-stuff. I'm afraid I never actually took a class in DSP. So laugh if you want. I can take it.
Anyway, the low-pass filter used last time was the following:
f*(t) = (2 / πτ) ∫ d Δt · cos²(Δt / τ) f(t ‒ Δt),
where here the limits of integration are again from ‒π/2 τ to +π/2 τ.
This can be turned into a band pass filter by simply multiplying the convolution function by Δt / τ:
f*(t) = (K / τ) ∫ d Δt · (‒Δt / τ) cos²(Δt / τ) f(t ‒ Δt).
A plot, after applying the delay by π/2 τ:
I suppose we should calculate K. The integral of the absolute value of the convolution function should be one. This is the first time I've use "absolute value of the integral" rather than just integral. It was a non-issue before since the convolution function was strictly non-negative. So I need to calculate twice the integral from 0 to π/2 of θ cos²(θ) dθ. This calls out for integration by parts, but honestly, I'm getting tired of integrating, so will seek professional assistance. The result is (π² / 4 - 1) / 2. The magic of the internet.... I could have done it myself. Honest. K is the reciprical of this: 8 / (π² - 4).
So then:
f*(t) = (8 / (π² - 4) τ) ∫ d Δt · (‒Δt / τ) cos²(Δt / τ) f(t ‒ Δt),
Again, the model for the function:
f(t) = f0 cos ωt,
I'll skip complex representations and the subsequent, again, integration by parts, which is guaranteed to yield an algebra error, and instead rely on my friend Stephen. Stephen doesn't do definite integrals for free, unfortunately, so I need to handle the limits myself. First, I'll simplify the input:
cos ω(t - Δt) = cos(ωt) cos(ω Δt) + sin(ωt) sin(ω Δt).
Then I recognize that by symmetry the first part yields zero in the integral. For the purposes of the code, I use a change in variables: w ≡ ωτ, x ≡ Δt / τ. Here's the result:
The integal is evaluated at the end points ‒π/2 to +π/2 for x:
1/2 [‒2 sin((π/2)ωτ) / ω²τ² + 2 (π/2) cos((π/2)ωτ) / ωτ ‒ sin((π/2)(ωτ ‒ 2)) / (ωτ ‒ 2)² ‒ sin((π/2)(ωτ + 2)) / (ωτ + 2)² + (π/2) cos((π/2)(ωτ ‒ 2)) / (ωτ ‒ 2) + (π/2) cos((π/2)(ωτ + 2))/ (ωτ + 2)].
This can be further simplified:
sin((π/2)(ωτ ± 2)) = ‒sin((π/2) ωτ),
cos((π/2)(ωτ ± 2)) = ‒cos((π/2) ωτ),
yielding:
1/2 [‒2 sin((π/2)ωτ) / ω²τ² + 2 (π/2) cos((π/2)ωτ) / ωτ + sin((π/2) ωτ) / (ωτ ‒ 2)² + sin((π/2)ωτ) / (ωτ + 2)² ‒ (π/2) cos((π/2)ωτ)/ (ωτ ‒ 2) ‒ (π/2) cos((π/2)ωτ) / (ωτ + 2)].
Then we can use
1/(ωτ ‒ 2)² + 1/(ωτ + 2)² = [(ωτ ‒ 2)² + (ωτ + 2)²] / (ω²τ² ‒ 4)² = 2 (ω²τ² + 4) / (ω²τ² ‒ 4)²
and
1/(ωτ ‒ 2) + 1/(ωτ + 2) = [(ωτ ‒ 2) + (ωτ + 2)] / (ω²τ² ‒ 4) = 2 ωτ / (ω²τ² ‒ 4)
Then the integral becomes:
sin((π/2)ωτ) [‒1 / ω²τ² + (ω²τ² + 4) / (ω²τ² ‒ 4)²] +
(π/2) cos((π/2)ωτ) [1 / ωτ ‒ ωτ / (ω²τ² ‒ 4)] =
sin((π/2)ωτ) [‒(ω²τ² ‒ 4)² + ω²τ² ω²τ² + 4)] / [ω²τ²(ω²τ² ‒ 4)²] +
(π/2) cos((π/2)ωτ) [(ω²τ² ‒ 4) ‒ ω²τ²] / [ωτ (ω²τ² ‒ 4)] =
sin((π/2)ωτ) (‒ω⁴τ⁴ + 8 ω²τ² ‒ 16 + ω⁴τ⁴ + 4 ω²τ²) / [ω²τ² (ω²τ² ‒ 4)²] ‒
2π cos((π/2)ωτ) ) / [ωτ (ω²τ² ‒ 4)] =
4 sin((π/2)ωτ) (3 ω²τ² ‒ 4) / [ω²τ² (ω²τ² ‒ 4)²] ‒
2 π cos((π/2)ωτ) ) / [ωτ(ω²τ² ‒ 4)] .
So then, putting this into the full expression:
f*(t) = (8 / (π² - 4)) [4 sin((π/2)ωτ) (3 ω²τ² ‒ 4) / [ω²τ² (ω²τ² ‒ 4)²] ‒ 2 π cos((π/2)ωτ) / [ωτ(ω²τ² ‒ 4)]] sin(ωt).
I admit I had to stare long and hard at this to convince myself it goes to zero as ωτ→0. But a Taylor series expansion of the tricky part, sin((π/2) ωτ) / ω²τ² + (π/2) cos((π/2) ωτ) / ωτ, is indeed proportional to ωτ for small values of ωτ, despite the apparent poles at zero due to the denominators.
So what about at high ωτ? The sine part is proportional to 1/(ωτ)⁴. The cosine part is proportional to 1/(ωτ)³. Therefore at high ωτ the cosine part will dominate, and the function will fall off proportional to (ωτ)³.
Wow -- that was a lot of work considering I "cheated" by using Mathematica for the integral. Maybe I would have been better off biting the bullet and doing the integral directly. Maybe, but I doubt it. But the end result is exciting: this filter is even more effective at suppressing high-frequency "noise" than the cosine squared filter was when that was combined with differentiation. That only fell off proportional to (ωτ)². And the icing is that no differentiation step is required.
I like this filter.
Okay, you can stop laughing now. I know this is all super-basic-DSP-stuff. I'm afraid I never actually took a class in DSP. So laugh if you want. I can take it.
Anyway, the low-pass filter used last time was the following:
f*(t) = (2 / πτ) ∫ d Δt · cos²(Δt / τ) f(t ‒ Δt),
where here the limits of integration are again from ‒π/2 τ to +π/2 τ.
This can be turned into a band pass filter by simply multiplying the convolution function by Δt / τ:
f*(t) = (K / τ) ∫ d Δt · (‒Δt / τ) cos²(Δt / τ) f(t ‒ Δt).
A plot, after applying the delay by π/2 τ:
I suppose we should calculate K. The integral of the absolute value of the convolution function should be one. This is the first time I've use "absolute value of the integral" rather than just integral. It was a non-issue before since the convolution function was strictly non-negative. So I need to calculate twice the integral from 0 to π/2 of θ cos²(θ) dθ. This calls out for integration by parts, but honestly, I'm getting tired of integrating, so will seek professional assistance. The result is (π² / 4 - 1) / 2. The magic of the internet.... I could have done it myself. Honest. K is the reciprical of this: 8 / (π² - 4).
So then:
f*(t) = (8 / (π² - 4) τ) ∫ d Δt · (‒Δt / τ) cos²(Δt / τ) f(t ‒ Δt),
Again, the model for the function:
f(t) = f0 cos ωt,
I'll skip complex representations and the subsequent, again, integration by parts, which is guaranteed to yield an algebra error, and instead rely on my friend Stephen. Stephen doesn't do definite integrals for free, unfortunately, so I need to handle the limits myself. First, I'll simplify the input:
cos ω(t - Δt) = cos(ωt) cos(ω Δt) + sin(ωt) sin(ω Δt).
Then I recognize that by symmetry the first part yields zero in the integral. For the purposes of the code, I use a change in variables: w ≡ ωτ, x ≡ Δt / τ. Here's the result:
The integal is evaluated at the end points ‒π/2 to +π/2 for x:
1/2 [‒2 sin((π/2)ωτ) / ω²τ² + 2 (π/2) cos((π/2)ωτ) / ωτ ‒ sin((π/2)(ωτ ‒ 2)) / (ωτ ‒ 2)² ‒ sin((π/2)(ωτ + 2)) / (ωτ + 2)² + (π/2) cos((π/2)(ωτ ‒ 2)) / (ωτ ‒ 2) + (π/2) cos((π/2)(ωτ + 2))/ (ωτ + 2)].
This can be further simplified:
sin((π/2)(ωτ ± 2)) = ‒sin((π/2) ωτ),
cos((π/2)(ωτ ± 2)) = ‒cos((π/2) ωτ),
yielding:
1/2 [‒2 sin((π/2)ωτ) / ω²τ² + 2 (π/2) cos((π/2)ωτ) / ωτ + sin((π/2) ωτ) / (ωτ ‒ 2)² + sin((π/2)ωτ) / (ωτ + 2)² ‒ (π/2) cos((π/2)ωτ)/ (ωτ ‒ 2) ‒ (π/2) cos((π/2)ωτ) / (ωτ + 2)].
Then we can use
1/(ωτ ‒ 2)² + 1/(ωτ + 2)² = [(ωτ ‒ 2)² + (ωτ + 2)²] / (ω²τ² ‒ 4)² = 2 (ω²τ² + 4) / (ω²τ² ‒ 4)²
and
1/(ωτ ‒ 2) + 1/(ωτ + 2) = [(ωτ ‒ 2) + (ωτ + 2)] / (ω²τ² ‒ 4) = 2 ωτ / (ω²τ² ‒ 4)
Then the integral becomes:
sin((π/2)ωτ) [‒1 / ω²τ² + (ω²τ² + 4) / (ω²τ² ‒ 4)²] +
(π/2) cos((π/2)ωτ) [1 / ωτ ‒ ωτ / (ω²τ² ‒ 4)] =
sin((π/2)ωτ) [‒(ω²τ² ‒ 4)² + ω²τ² ω²τ² + 4)] / [ω²τ²(ω²τ² ‒ 4)²] +
(π/2) cos((π/2)ωτ) [(ω²τ² ‒ 4) ‒ ω²τ²] / [ωτ (ω²τ² ‒ 4)] =
sin((π/2)ωτ) (‒ω⁴τ⁴ + 8 ω²τ² ‒ 16 + ω⁴τ⁴ + 4 ω²τ²) / [ω²τ² (ω²τ² ‒ 4)²] ‒
2π cos((π/2)ωτ) ) / [ωτ (ω²τ² ‒ 4)] =
4 sin((π/2)ωτ) (3 ω²τ² ‒ 4) / [ω²τ² (ω²τ² ‒ 4)²] ‒
2 π cos((π/2)ωτ) ) / [ωτ(ω²τ² ‒ 4)] .
So then, putting this into the full expression:
f*(t) = (8 / (π² - 4)) [4 sin((π/2)ωτ) (3 ω²τ² ‒ 4) / [ω²τ² (ω²τ² ‒ 4)²] ‒ 2 π cos((π/2)ωτ) / [ωτ(ω²τ² ‒ 4)]] sin(ωt).
I admit I had to stare long and hard at this to convince myself it goes to zero as ωτ→0. But a Taylor series expansion of the tricky part, sin((π/2) ωτ) / ω²τ² + (π/2) cos((π/2) ωτ) / ωτ, is indeed proportional to ωτ for small values of ωτ, despite the apparent poles at zero due to the denominators.
So what about at high ωτ? The sine part is proportional to 1/(ωτ)⁴. The cosine part is proportional to 1/(ωτ)³. Therefore at high ωτ the cosine part will dominate, and the function will fall off proportional to (ωτ)³.
Wow -- that was a lot of work considering I "cheated" by using Mathematica for the integral. Maybe I would have been better off biting the bullet and doing the integral directly. Maybe, but I doubt it. But the end result is exciting: this filter is even more effective at suppressing high-frequency "noise" than the cosine squared filter was when that was combined with differentiation. That only fell off proportional to (ωτ)². And the icing is that no differentiation step is required.
I like this filter.
Comments