Taking another crack at basic filtering

[thevinn]

I was exhausted after my previous efforts to understand the Biquads and tune my filters so they sounded right (never did get that working correctly). Instead I satiated myself with having VST plugins working so they did the filtering job fror me.

Unfortunately in order to apply some of the more advanced techniques like resampling (e.g.44100<->48000<->96000) or transient/beat detection, these darn filters keep appearing! Its like every time I turn my head I see low pass, high pass, band pass, all pass, band reject filters coming up!!

So I am taking on a second round of trying to tackle these things.

[thevinn]

I've decided to take a step back and write a completely new test application that will graph the frequency response of my filters, with parameterized knobs to see the effects of changes in real time.

How does one graph the frequency response (or is it called magnitude response) of a filter? Do I generate a unit sine wave at each frequency, run it through the filter, and then measure the resulting output level?

How much of the sine wave do I need? Is just 2*pi radians enough?

How do I scale the X axis of the graph correctly? My hunch is that a linear scaling is not desired...I've seen graphs and they do something logarithmic with the frequencies. I'm not sure exactly what though.

Any pointers are greatly appreciated.

[bmoviehorror]

What are you using language wise to do this? Matlab has a function called freqz which makes visualisation of your filter's transfer functions really easy. If you don't have access to Matlab then the Scipy module in Python has equivalent functions to do this with minimal effort. If you're tring to test this out in C or something similar then you are much braver than I.

[thevinn]

I'm doing this all in C++!!!!

It was relatively easy for me to start with a new WinMain() and an empty window. I whipped up a graph control and re-used the knobs from my main application. This is what I have so far:



The first 4 knobs control the origin and scale of the graph.

I'm working on the frequency response now.

I prefer to use a C++ application because then I can guarantee that the values I am using in the actual application are exactly the same as the values that I am graphing. So any bugs will be clearly visible. And I will be able to step into the routines if I get a NaN (which seems to be happening more often than I would like).

[docdued]

The frequency response of your filter is the Fourier transform of its output when you feed it with an array x[n] that has x[0]=1 and x[n]=0 otherwise (this is basically "a sine wave at each frequency"). And yes, indeed it's common practice to plot the frequency axis logarithmically (gives you the same length on the axis per octave).

[thevinn]

I'm a little unsure of how to proceed now.

All of my filter objects take in a parameter called "frequencyNormalized":

Code: Select all

		void BiQuad::SetParams(
			double q,
			double frequencyNormalized ); // 0 <= freq <= 0.5

frequencyNormalized is basically the frequency setting of the filter (for example, the cutoff for a low-pass) divided by the sample rate.

I have to produce a sine wave for each frequency I want to test in my response graph. But is it necessary for me to literally create a new sine wave every time? Or can I produce the sine wave once at some given size (lets say, 2*pi radians over 4096 samples), and re-use that wave over and over again, calling the filter with a different value for frequencyNormalized?

[duncanparsons]

really, do what docdued suggested: essentially

Code: Select all

float buffer[4096]

filter->freq = somevalue;
filter->q = somevalue;
filter->reset;
buffer[0] = filter->process(1.f);
for(int i = 1; i<4096; i++)
{ buffer[i] = filter->process(0.f); }

//buffer now has the impulse response
processFFT(buffer);

you then display the results. Laurent de Soras's FFTReal is a good simple one to use.

I have a Delphi program set up to do just this kind of thing

DSP

[thevinn]

I got the thing working, need to make some adjustments now. How do I do the conversion to logarithmic units?

My hunch about the sine() turned out to be correct. I produced a fixed size buffer (256 samples currently) once, and fill it with a full sine wave (2*pi radians).

Every filter usually expresses its frequency parameter as (filterFrequency/sampleRate). I'm sure you are all familiar with this. In order to calculate the apparent sample rate for the fixed length sine buffer, given the frequency on the graph (not to be confused with the frequency parameter to the filter; for example, the cut-off for a low-pass), it is

Code: Select all

sampleRate = frequency * sampleCount

Check it out:



X axis is frequency from 20Hz to 20,000Hz. I used a low pass Biquad. Unfortunatey the x axis is linear...still working on that.

It vaguely looks like what you would expect from a low pass filter and when I tweak the frequency knob and Q paramter it does what you would expect.

[Music Engineer]

FFTing the impulse response works, but you will be limited to the FFT-bins (which are inherently linearly scaled along the frequency axis). if you want to scale your frequency axis logarithmically, you will clearly see the low resolution at low frequencies (check out FFT analyzers such as my signal analyzer). furthermore, there will be artifacts due to truncation of an otherwise infinite impulse response. a more elegant approach would be to just evaluate the magnitude response at each frequency of interest. then, frequencies of interest may be distributed along the axis in whatever way you want (logarithmically, most likely). there are analytic expressions for the magnitude response in terms of the filter coefficients. you may also evaluate the complex frequency response and then take the magnitude.

[thevinn]

a more elegant approach would be to just evaluate the magnitude response at each frequency of interest. then, frequencies of interest may be distributed along the axis in whatever way you want (logarithmically, most likely).

That is exactly what I am doing.

there are analytic expressions for the magnitude response in terms of the filter coefficients. you may also evaluate the complex frequency response and then take the magnitude.

You are saying that it is possible just given the filter coefficients (assuming it is one of the well known filter types) to produce the magnitude response? Well that makes sense. However it won't reveal implementation errors. It is also not generalized so if you want to change the filter, you have to come up with a new closed-form expression for the magnitude response. Plus you lose the flexibility of analyzing other useful output like filters chained together.

What, exactly, is the "magnitude response"? My sine wave has its highest amplitude equal to 1.0. After I filter it, I look through the output to find the highest value and that is what I am plotting. So if the filter were to divide every input sample by two, I would plot 0.5 on the graph.

Is this the right thing to do? Or should I be using some formula, like the area under the curve?

What about the units for the y axis should they be linear? All the diagrams I see talk about the dB response. And there is some magic number -12dB that keeps coming up. I also see -3dB and the word "octave" used together often...can anyone please shed light on these mysteries?

[Music Engineer]

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);
}

omega == 2*pi*frequency/sampleRate here (normalized radian frequency). ...without any warranties. let me know if you find errors.

[thevinn]

maybe read the 'Filter Basics' tutorial on my website. in the meanwhile i'll try to dig out the formulas....

Its making much more sense now (I've read it a few times in the past).

I'm still a little bit confused about magnitude response as it relates to my application.

A filter that passes the signal through unchanged, has a magnitude response of 1. A filter that silences its output completely has a magnitude response of zero.

So, give a block of output samples from a particular filter how do we compute the magnitude response if we know the input? In my case the input is a sine wave with float values between -1.0 and 1.0. Currently I am just going through the output and taking the highest value and assuming that is the magnitude response.

[Music Engineer]

maybe read the 'Filter Basics' tutorial on my website. in the meanwhile i'll try to dig out the formulas....

Its making much more sense now (I've read it a few times in the past).

I'm still a little bit confused about magnitude response as it relates to my application.

well, there's also the phase response. the complex frequency response would then comprise both: magnitude and phase and it's only the complex frequency response that characterizes the filter fully. filters with equal magnitude responses may still differ in their phase response. ... terms like minimum-phase, allpass part, whatever, .... will then crop up

A filter that passes the signal through unchanged, has a magnitude response of 1. A filter that silences its output completely has a magnitude response of zero.

yup. but always think of it as function of frequency, of course.

So, give a block of output samples from a particular filter how do we compute the magnitude response if we know the input? In my case the input is a sine wave with float values between -1.0 and 1.0. Currently I am just going through the output and taking the highest value and assuming that is the magnitude response.

it is. but it's not a very elegant way to evaluate it.

[Music Engineer]

arrg. ...wanted to quote myself and edited my post instead. O.K. - so see above for the formulas and some C++ code. damn!

[thevinn]

Sorry for such a basic question but in your document AmplitudeAndDecibels.pdf you write:

Code: Select all

L = 20 * log10( A / A0 )

In our case, A0 is 1 and A is the magnitude response (correct me if I'm wrong).

Why the number 20? Shouldn't that be 10? Where did 20 come from? I thought a decibel was a tenth of a bel and not a twentieth...