Forums

How to create peak/shelf filters, that are symmetric close to nyquist?

Started by jungledmnc November 3, 2011
Hi,

is it possible to modify biquad peak/shelf filters, so that they have
approximately the same shape allover the spectrum. Normally they start to
"narrow" the closer they are to nyquist. I have checked it can more or less
be compensated by combining multiple filters, but that's pretty much an
alchemy :), so I wonder if there isn't some regular approach.

Thanks in advance.
jungledmnc

jungledmnc wrote:

> is it possible to modify biquad peak/shelf filters, so that they have > approximately the same shape allover the spectrum. Normally they start to > "narrow" the closer they are to nyquist. I have checked it can more or less > be compensated by combining multiple filters, but that's pretty much an > alchemy :), so I wonder if there isn't some regular approach.
Q dewarping: x = 2*Pi*F/Fsa; Q_dewarped = Q * sin(x)/x; Vladimir Vassilevsky DSP and Mixed Signal Design Consultant http://www.abvolt.com
>Q dewarping: > >x = 2*Pi*F/Fsa; >Q_dewarped = Q * sin(x)/x;
Thanks Vladimir. But I don't think that's enough. I'm mostly interested in the peak filter. The trouble here is that it will get to 0dB at nyquist, while for symmetric response both side should be equal, so not only the Q must be changed, because the "left" side gets narrower as well, but also the "right" side must be fixed somehow, because at Nyquist the requested response might for example be +3dB. I thought about adding a high shelf filter and it could probably work, but as I said, it's an alchemy... jungledmnc
On 11/3/11 11:00 PM, jungledmnc wrote:
>> Q dewarping: >> >> x = 2*Pi*F/Fsa; >> Q_dewarped = Q * sin(x)/x; > > Thanks Vladimir. But I don't think that's enough. I'm mostly interested in > the peak filter. The trouble here is that it will get to 0dB at nyquist, > while for symmetric response both side should be equal, so not only the Q > must be changed, because the "left" side gets narrower as well, but also > the "right" side must be fixed somehow, because at Nyquist the requested > response might for example be +3dB. I thought about adding a high shelf > filter and it could probably work, but as I said, it's an alchemy... >
do you know about the Orfanidis design? he puts the gain at Nyquist to be the gain of the analog prototype at the same frequency, rather than to be 0 dB. he did a paper following an old one i did in 1994. it's in the AES journal somewhere. i might even have simplified design equations, but i would have to look around. but there is no way to make it perfect because the slope of the gain function at Nyquist must necessarily be 0 for a digital filter (and it isn't for the analog prototype). -- r b-j rbj@audioimagination.com "Imagination is more important than knowledge."
Thanks RBJ, found it, now I need to get into it :)).

Btw. its name is "Digital Parametric Equalizer Design With Prescribed
Nyquist-Frequency Gain"

jungledmnc
>Thanks RBJ, found it, now I need to get into it :)). > >Btw. its name is "Digital Parametric Equalizer Design With Prescribed >Nyquist-Frequency Gain" > >jungledmnc >
I converted the MatLab code to C++ from the Orfanidis' paper, but the results aren't really what's expected. I use these inputs: const double G0 = 1.0; const double G = mpow(10.0, DBGain / -20.0); const double GB = mpow(10.0, DBGain / -40.0); const double w0 = 2 * MPI * Frequency / SampleRate; const double Dw = 2 * w0 * msinh( (msin(w0)/w0) * masinh (1/(2*Q))) ; Also I'd like the filter to have the same parameters as r-b-j's peak filter, so the bandwidth definition is kinda fuzzy. I originally tried something else, but googled this one. What it does is, well, weird :). Any idea what could be wrong? I might have done something wrong when converting the code, but I checked several times, so I'm mostly asking if the inputs make sense. Thanks! jungledmnc
On 11/9/11 9:02 PM, jungledmnc wrote:
> > I converted the MatLab code to C++ from the Orfanidis' paper, but the > results aren't really what's expected. I use these inputs: > > const double G0 = 1.0; > const double G = mpow(10.0, DBGain / -20.0); > const double GB = mpow(10.0, DBGain / -40.0); > const double w0 = 2 * MPI * Frequency / SampleRate;
except for the minus signs on 20.0 and 40.0, these look reasonable. do you want to boost or cut by DBGain?
> const double Dw = 2 * w0 * msinh( (msin(w0)/w0) * masinh (1/(2*Q))) ;
i have to think about this one. i think it should be a little different. from the cookbook: 1/Q = 2*sinh( ln(2)/2 * BW * w0/sin(w0) ) BW in octaves and from any of a number of texts: Q = w0/Dw but, i think if you widen the bandwidth because of bilinear "cramping", then it's Q = w0/( w0/sin(w0) * Dw ) = sin(w0) / Dw
> Also I'd like the filter to have the same parameters as r-b-j's peak > filter, so the bandwidth definition is kinda fuzzy. I originally tried > something else, but googled this one. > > What it does is, well, weird :). Any idea what could be wrong? I might have > done something wrong when converting the code, but I checked several times, > so I'm mostly asking if the inputs make sense. > >
is this the matlab code you mean? % peq.m - Parametric EQ with matching gain at Nyquist frequency % % Usage: [b, a, G1] = peq(G0, G, GB, w0, Dw) % % G0 = reference gain at DC % G = boost/cut gain % GB = bandwidth gain % % w0 = center frequency in rads/sample % Dw = bandwidth in rads/sample % % b = [b0, b1, b2] = numerator coefficients % a = [1, a1, a2] = denominator coefficients % G1 = Nyquist-frequency gain % % Available from: www.ece.rutgers.edu/~orfanidi/intro2sp/mdir/peq.m function [b, a, G1] = peq(G0, G, GB, w0, Dw) F = abs(G^2 - GB^2); G00 = abs(G^2 - G0^2); F00 = abs(GB^2 - G0^2); num = G0^2 * (w0^2 - pi^2)^2 + G^2 * F00 * pi^2 * Dw^2 / F; den = (w0^2 - pi^2)^2 + F00 * pi^2 * Dw^2 / F; G1 = sqrt(num/den); G01 = abs(G^2 - G0*G1); G11 = abs(G^2 - G1^2); F01 = abs(GB^2 - G0*G1); F11 = abs(GB^2 - G1^2); W2 = sqrt(G11 / G00) * tan(w0/2)^2; DW = (1 + sqrt(F00 / F11) * W2) * tan(Dw/2); C = F11 * DW^2 - 2 * W2 * (F01 - sqrt(F00 * F11)); D = 2 * W2 * (G01 - sqrt(G00 * G11)); A = sqrt((C + D) / F); B = sqrt((G^2 * C + GB^2 * D) / F); b = [(G1+ G0*W2 + B), -2*(G1 - G0*W2), (G1 - B + G0*W2)] / (1 + W2 + A); a = [1, [-2*(1 - W2), (1 + W2 - A)] / (1 + W2 + A)]; he says: "Its inputs are the gains G0, G, GB in absolute units, and the digital frequencies ω0, Δω in units of rads/sample." so how are you relating those to the given parameters in the cookbook? even though he mentions using the w0/sin(w0) bandwidth expansion (to compensate the scrunching of BW that the bilinear transform does) in the text, i don't see it in this code. maybe someone can point it out to me. a decade ago i've dealt with all of this. dunno where those equations (or code) is, but i might stumble upon it in the next 48 hours. -- r b-j rbj@audioimagination.com "Imagination is more important than knowledge."
On 10 Nov, 03:02, "jungledmnc" <jungledmnc@n_o_s_p_a_m.gmail.com>
wrote:

> I converted the MatLab code to C++ from the Orfanidis' paper, but the > results aren't really what's expected.
...
> const double Dw = 2 * w0 * msinh( (msin(w0)/w0) * masinh (1/(2*Q))) ;
Never use nested expressions like these when coding C++. The expressions msinh, msin etc are likely template expressions, which is 'partial' library code that lacks a couple of type details to be filled by the user of the library. This means that the code will be completed at compile time, at which time the optimizer is let loose on the expression. This works very well if there is *one* template expression per numerical C++ statement (between two semicolons). For some reason things might not work all that well if there are more than one template expression per statement. I don't know why, but it might be that the optimizer gets confused and tries and mangle elements of several template expressions into one big mess. Try and reorganize the line above as double xx = (msin(w0)/w0); double yy = masinh (1/(2*Q)); const double Dw = 2 * w0 * msinh( xx * yy) ; This way you limit the amount of code fed to the optimizer at any one time, vastly reducing the probability of it creating a mess. Rune
Thanks folks. The thing is, whatever I do with the input parameters, the
results just don't look like a peak filter. I'll have to make Octave work
and check the parameters I guess :).
On 10 Nov, 14:13, "jungledmnc" <jungledmnc@n_o_s_p_a_m.gmail.com>
wrote:
> Thanks folks. The thing is, whatever I do with the input parameters, the > results just don't look like a peak filter. I'll have to make Octave work > and check the parameters I guess :).
Matlab & clones usually take a lot of beating in this forum, so for once I will praise it: This is exactly the kind of thing matlab (or octave or scilab or...) was invented for! You, the user, reads some document about some ethod and want to try it. Instead of diving deep into C or Fortran or assembler or whatever straight away, you first use matlab (or octave or scilab or...) to check that you get the computations right. Only then, when you are happy that you get the results you want, do you tackle the challenges of C++ template programming. Rune