Taking another crack at basic filtering

[Music Engineer]

I've switched over to taking the average of the sum of the absolute values of the samples as the magnitude ("a"). It incorporates every sample. I also compute this for the original sine wave ("a0") for the formula

Code: Select all

dB = 20 * log ( a / a0 )

hhmm... well... that may even out some of the errors. but why not just use the exact value as given by the closed form formula (which is also presumably much easier and faster to evaluate)?

[thevinn]

why not just use the exact value as given by the closed form formula (which is also presumably much easier and faster to evaluate)?

Because I don't know which formula. Is it on www.rs-met.com ? Which formula is it?

[Music Engineer]

In my formula I am using -6dB = half power (i.e. 20*log(a/a0)).

-6dB = half amplitude, -3dB = half power

The orange line is a low pass Biquad. I've checked and triple checked the coefficients (even found a bug in an online published version of the cookbook).

wott? a bug in the RBJ formulas as published in:

http://www.musicdsp.org/files/Audio-EQ-Cookbook.txt

i'm using these formulas since ages without any problems.

Is ~20dB drop at the cut-off frequency the expected result from a second order biquad lowpass?

that's a bit much. i'd expect -3 dB but i'm not quite sure at the moment if the RBJ filters adhere to this convention ...the Q/resonance stuff makes things more complicated. i think, at Q=1/sqrt(2) it should be at -3 dB at the cutoff. but not entirely sure.

[thevinn]

I'm starting to wonder if getBiquadMagnitudeAt() is the right formula.

I found this other equation from http://groups.google.com/group/comp.dsp ... 63ac56b686 and I implemented it:

Code: Select all

double getBiquadMagnitudeAt2( double b0, double b1, double b2, double a0, double a1, double a2, double normalizedFrequency )
{
	double a, b;
	double phi;

	phi=sin(M_2_PI*normalizedFrequency);
	phi=phi*phi;
	a=(::pow(b0+b1+b2, 2.0)-4*(b0*b1+4*b0*b2+b1*b2)*phi+16*b0*b2*phi*phi);
	b=(::pow(a0+a1+a2, 2.0)-4*(a0*a1+4*a0*a2+a1*a2)*phi+16*a0*a2*phi*phi);

	return a/b;
}

It gives me different results. The shape of the curve is the same, but the cutoff frequency is closer to where I would expect it to be. But still not where it should be (I think).

[thevinn]

http://www.musicdsp.org/files/Audio-EQ-Cookbook.txt

i'm using these formulas since ages without any problems.

That one is fine but the original PDF has an error in the C-code example:

http://music.ucsd.edu/~tre/biquad.pdf

For case LP (lowpass) they are setting

Code: Select all

a2 = alpha - 1.0F

which clearly, is wrong (and inconsistent with the previous text).

[Music Engineer]

why not just use the exact value as given by the closed form formula (which is also presumably much easier and faster to evaluate)?

Because I don't know which formula. Is it on www.rs-met.com ? Which formula is it?

for a biquad it would be the getBiquadMagnitudeAt function as posted here. for other filters, similar formulas can be worked out. for the 1st order filter, you find them in my pdf (posted here). for the general case of a direct form filter, there's a formula with finite sums in the numerator and denominator which i currently don't have handy.

[thevinn]

for a biquad it would be the getBiquadMagnitudeAt function as posted here.

Oh..I thought you were referring to the function that computes the magnitude response from a given buffer of samples.

[Music Engineer]

for a biquad it would be the getBiquadMagnitudeAt function as posted here.

Oh..I thought you were referring to the function that computes the magnitude response from a given buffer of samples.

the more i think about it, the wronger this approach appears to me. think about an allpass that shifts your sinusoid by a fractional amount of samples. that will certainly alter the the maximum absolute value and quite probably also the average. if you really need a method for measuring frequency responses of filters that are a black box to you, i'd recommend to stick to the FFT of the impulse response

[thevinn]

if you really need a method for measuring frequency responses of filters that are a black box to you, i'd recommend to stick to the FFT of the impulse response

I'm completely exhausted from these efforts. I'm going to rebuild my self-esteem by implementing some more easy (and necessary) features while I recover. I didn't want to start with the FFT so soon but it seems to be the best option now. I will revisit this topic in the future and try the FFT approach.

In the meanwhile, thanks to everyone for their help. Using what I learned I was able to find a bug in my high pass filter, and tweak my equalizer so it produces the results I was looking for from the start (Re: "totally awesome cplusplus DSP code"). My bipolar filter is also improved, and the sliders and knobs feel more useful now with new decibel->gain mappings.

[Music Engineer]

I'm completely exhausted from these efforts. I'm going to rebuild my self-esteem by implementing some more easy (and necessary) features while I recover. I didn't want to start with the FFT so soon but it seems to be the best option now. I will revisit this topic in the future and try the FFT approach.

you probably don't want to implement your own FFT then and you don't have to. there are good libraries out there. ooura, kiss-fft, FFTReal (i think), stephan bernsee's and probably others as well. good luck.

[Music Engineer]

...for the general case of a direct form filter, there's a formula with finite sums in the numerator and denominator which i currently don't have handy...

...update: i digged out the formula for the magnitude (and phase) response for the general direct form filter and updated the pdf document accordingly:

www.rs-met.com/documents/dsp/BasicDigitalFilters.pdf

also, in order to assure correctness, i implemented it in a small matlab/octave function:

Code: Select all

function [m p] = directFormMagnitudeAndPhaseAt(b, a, w);

% Returns the value of the magnitude- and phase-response at the normalized 
% radian frequency w (w = 2*pi*frequency/sampleRate) of the direct form 
% filter defined by the filter coefficient vectors b and a with the 
% difference equation:
%
% y[n] = b[1] x[n] + b2 x[n-1] + ... - a[2] y[n-1] - a[3] y[n-2] - ...
%
% usage:
%  [m p] = directFormMagnitudeAndPhaseAt(b, a, w);
%
% input-variables:
%  -b: vector of the feedforward coefficients
%  -a: vector of the feedback coefficients
%  -w: normalized radian frequency w 
%      (w = 2*pi*frequency/sampleRate) at which the magnitude
%      response id to be evaluated
%
% output-variables:
%  -m: value of the magnitude response at the given w
%  -p: value of the phase response at the given w

%--------------------------------------------------------------------------

M  = length(b);
sb = 0;
cb = 0;
for k=0:M-1
  sb = sb + b(k+1) * sin(k*w);
  cb = cb + b(k+1) * cos(k*w);
end

N  = length(a);
sa = 0;
ca = 1;
for k=1:N-1
  sa = sa + a(k+1) * sin(k*w);
  ca = ca + a(k+1) * cos(k*w);
end

m = sqrt( (cb^2+sb^2) / (ca^2+sa^2) );
p = atan2(sb,cb) - atan2(sa,ca);

% wrap phase into range -pi...pi
if p > pi 
  p = p-2*pi;
elseif p < -pi
  p = p+2*pi;
end

and compared it with an example filter against the output of freqz. note that one normally would avoid loops in octave and would use vector expressions instead for efficiency, but i'm not so much concerned about it in a proof-of-concept implementation and it's also more obvious how to port it to C++ this way. one further significant optimization would be to evaluate the sines and cosines by trigonometric recursion ("fast sine osclillator" alike) but i couldn't be bothered with that right now.

[Music Engineer]

wait...there seems to be a sign inversion error with the phase response...i'll have to investigate this...

[Bluebirdglory]

in the meanwhile i'll try to dig out the formulas....

O.K. here's some pdf where i noted down the formulas for the magnitude responses of some low order filters

www.rs-met.com/documents/dsp/BasicDigitalFilters.pdf

here's the corresponding C++ code for the biquad:

Code: Select all

double getBiquadMagnitudeAt(double b0, double b1, double b2, 
                            double a1, double a2, double omega)
{
  double cosOmega  = cos(omega);
  double cos2Omega = cos(2.0*omega);
  double num       = b0*b0 + b1*b1 + b2*b2 
                     + 2.0*cosOmega*(b0*b1 + b1*b2) 
                     + 2.0*cos2Omega*b0*b2;
  double den       = 1.0   + a1*a1 + a2*a2 
                     + 2.0*cosOmega*(a1    + a1*a2) 
                     + 2.0*cos2Omega*a2;

  return sqrt(num/den);
}

May i know if we want the formula for higher orders,( 3 or 4 ) how the formula will be ? thanks

[thevinn]

May i know if we want the formula for higher orders,( 3 or 4 ) how the formula will be ? thanks

Full source code for what you ask, and much much more:

http://code.google.com/p/dspfilterscpp/

The function you are looking for is called Response(). Use one of the predefined filters, or load a cascade with your own coefficients and then call Response() to get the complex-valued response at angular frequency w. From there you can use abs() to get the magnitude response and phase() to get the phase response.