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

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