Help with matched bandpass filter *Solved* (Ensuring center frequency is maximum)

[matt42]

Hi, so I decided to have a crack at a matched band pass filter. I'm basing it on the Biquad Fits and Mzti papers.

My math and dsp is self taught, so I have a hard time understanding how in the Biquad Fits paper they arrive at:

Code: Select all

1. Single zero at DC: b0 + b1 + b2 = 0.
2. Maximum at ω0:
d
dω |H(z)|
2
|ω=ω0 = 0.
3. Unity gain at ω0: |H(z)|
2
|ω=ω0 = 1
The first requirement yields B0 = 0. The second and third requirements
imply
B1 + 4(φ0 − φ1)B2 − A0 + A1 + 4(φ0 − φ1)A2 = B1φ1 + B2φ2 A0φ0 + A1φ1 + A2φ2 = 1. (

Specifically, how the second two requirements led to the equation.

Anyway, I figured to try something similar, because a problem with MZT is aliasing of the response, so that for band pass the peak value is off center.

I came up with the following:

1: find poles via MZT
2: add zero at DC
3: calculate gain squared at the center frequency and it's derivative there
4: Find an additional zero using a curve fit that sets the overall derivative of the system to zero at the center frequency
5: Combine the two zeros and adjust gain so center frequency is at unity

Here's the curve fit:

Code: Select all

//sorry for lazy variable naming g = gain squared and d = derivative
(sqrt((g*g - d*d)*sinw*sinw - 2.0*d*g*sinw) + g*sinw - d*cosw) / d

which sets one of the coefficients. The other, for convenience, is 1.0 - overall gain is later compensated and I'm just interested in fixing the derivative at that point. If d ==0 then we avoid division by zero simply because the derivative is correct and we don't need the curve fit.

The problem is the sqrt() function. Everything works great, but once above 10 kHz the input to the sqrt() function quickly becomes massively negative (d outgrows g and due to the squaring the difference quickly becomes very large).

So, is this approached doomed, or should I take a look at adjusting the equation somehow?

Thanks,

Matt

[matt42]

Hi, so don't do filter maths whilst staying up till 4am to watch some stupidity

Anyway I've found a solution that under brief testing seems stable and produces the correct result.

Code: Select all

//current correction filter will initially have gain of 1.0 at f0, 
//probably could bake gain correction in at this stage

targetDerivative = d = bandPassGainAtF0/-bandPassDerivativeAtF0;
//haven't checked if result is minimum phase - possibly need to swap a[0] and a[1]

a1 = sqrt(sqrt(-d*d*sinw*sinw + sinw*sinw + 2.0*d*cosw*sinw)/sinw + (d*cosw/sinw + 1.0)) / 1.414213562373;

a0 = d/(-2.0*sinw*a1);

//The last step are to normalise the gain and combine with the zero at DC for a single second order structure

A shame it needs two sqrt() calls and will need further testing to check if some random parameter settings blow it up.

Anyway, is there a better way? Maybe a simpler way to set f0 as the maximum of the filter? Also I could share full source if any ones interested, though I think the calculation of coefficients will be a little heavy and use at your own risk of course.

Cheers,

Matt

Edit: Decided to add a couple of response curves:



I'm pretty happy with it. You can see some response overshoot at 20 kHz of around a couple of dB on the 10 kHz filter. One more order of correction and I think it'd be near perfect, though not sure if it'd be worth going to a third of fourth order filter for that.

[juha_p]

Martin Vicanek has coded his paper (linked in OP) for Flowstone http://www.dsprobotics.com/support/view ... f=4&t=3710 (MatchedBiquads.fsm there).

[matt42]

Thanks Juha_P,

If I get around to optimizing my implementation it might be interesting to compare with Martin Vicanek's. Also I'm assuming frequency response should be identical (I have used/borrowed his design requirements after all), but would be good to compare the efficiency of the coefficient calculations.

My main question wasn't about implementation of the paper, however. I just don't often follow the logic/mathematics when you get the dreaded sentence "the requirements imply..." then blam - equation! This is probably as I'm not a mathematician. I wonder did Martin Vicanek set the derivative to zero at the center frequency, or perhaps use a smarter method to fix the maximum there? I'm not enough of a mathematician (I'm not a mathematician at all ) to know which methods might be feasible, or to really unravel his equations.

I just went for it with the derivative method and it seems to work nicely (and I had assumed it should for these kinds of relatively simple/low order curves), what other approaches could be used?

Edit: I should use my eyes better. I see that his requirement is for a derivative of zero at the center frequency - my bad

[matt42]

The main difference is that I'm only solving for one zero where as the paper solves both simultaneously. I, maybe, make some savings this way - will have to finalise my version before I can make any kind of comparison. Combining zeros is easy enough. For my initial band pass I place a zero at DC so, for simplicity, the feed forward coefficients are [1,-1]. After solving the correction zero [a0, a1] it's a simple convolution to combine them: a0 stays as it is, a2 = -a1 and a1 -= a0.

[JCJR]

That is interesting matt42. I'm vastly dumber on math.

An idle curiosity if you get around to testing this-- I don't NEED to know, just curious-- If a purpose of designing the beautiful non-warped bandpass response might be to mix it with the original signal to derive a beautiful non-warped Peak/Dip EQ filter-- It would be interesting to know if the resultant peaking filter turns out as equally non-warped as is the bandpass filter?

Am just wondering, because peaking filter response depends as much (or more) on phase reinforcement/cancellation when mixed with the input, as it depends on the bandpass filter's amplitude response. So am curious whether the changed phase response of the un-warped bandpass filter would "mess up" the symmetry of a peaking filter's skirts?

Or is that a silly consideration, a non-issue, because if one wants a non-warped peaking filter then one would just use slightly different solving method customized for a peaking filter, rather than simple-minded mix a bandpass output with the input?

[matt42]

Hi JCJR,

Unfortunately, using the MZT for the poles and curve fitting the zeros doesn't preserve the phase/amplitude relationship like the BLT, so, I think, as you get closer to nyquist the peak filter response (from mixing bandpass with input) will get progressively worse. So, yeah, for a matched response it would be best to roll a specific peak implementation.