KVR Audio

mystran · Post by **mystran** » Sat Apr 18, 2009 6:38 am

Uhm.. exactly how did you manage to get ripple like that from a 2nd order filter?

I'd guessed it'd be an artifact of short FFT window but it seems to get more dense (rather than the other way around) at lower frequencies.. so what the hell?

Are you sure your graphing is correct?

thevinn · Post by **thevinn** » Sat Apr 18, 2009 7:01 am

mystran wrote:Uhm.. exactly how did you manage to get ripple like that from a 2nd order filter?

I'd guessed it'd be an artifact of short FFT window but it seems to get more dense (rather than the other way around) at lower frequencies.. so what the hell?

Are you sure your graphing is correct?

There's no FFT involved.

I'm using a low pass biquad

Code: Select all

void LowPass::SetParams( double q, double frequencyNormalized )
{
	double w0, cs, sn, al;

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

	w0 = M_2_PI * frequencyNormalized;
	cs = cos(w0);
	sn = sin(w0);

	al = cs / ( 2 * q );
	
	b0 = ( 1 - cs ) / 2;
	b1 =   1 - cs;
	b2 = ( 1 - cs ) / 2;
	a0 =   1 + al;
	a1 =  -2 * cs;
	a2 =   1 - al;

	Optimize();
}

void BiQuad::Optimize( void )
{
	// optimize filter coefficients
	a1/=a0;
	a2/=a0;
	b0/=a0;
	b1/=a0;
	b2/=a0;
	a0 =1;
}

I have a high degree of confidence in the y axis coordinates. When I use the "3 band equalizer" and set the low gain to 0.25 I get a result of -12dB cut in the lows. Here is a graph:

At 0.5 gain I get -6 dB. So I am very confident about the y axis.

The x-axis I am also very confident about. When I use the following code to return the ideal magnitude (green line on graphs);

Code: Select all

	if( frequency<=params->f )
		*idealMagnitude=1;
	else
		*idealMagnitude=0;

The result is a perfect vertical line at the expected frequency. You can see this in image with the ripple.

Here is the algorithm for calculating the frequency response graph:

Code: Select all

	int i;
	double p=log10( freqHi-freqLo );
	for( i=0;i<count;i++ )
	{
		double freq=freqLo+::pow( 10, p*i/(count-1) );
		double sampleRate=freq*double(m_nsine);

		::memcpy( m_out, m_sine, m_nsine*sizeof(m_sine[0]) );

		double ideal=0;
		if( m_proc!=0 )
			m_proc( m_nsine, m_out, sampleRate, freq, &ideal, m_opaque );

		if( ideal<0.00001 )
			m_ideal[i]=-1000;
		else
			m_ideal[i]=20*log10(ideal);

		double a=calcmax( m_nsine, m_out );
		
		if( a<0.00001 )
			m_level[i]=-1000;
		else
			m_level[i]=20*log10(a/m_a0);
	}

where m_nsine is the number of samples in the precalculated sine wave buffer. freqHi and freqLo are the min and max frequencies on the graph (20 and 20,000 respectively). "count" is the width of the graph in pixels. When the routine is done, m_ideal is left with the ideal filter response in dB and m_level is the actual filter response in dB.

thevinn · Post by **thevinn** » Sat Apr 18, 2009 7:08 am

I am beginning to suspect that one or more of the formula I am using are vulnerable to numerical instability. Possibly getBiquadMagnitudeAt() for certain settings, and perhaps my biquad coefficients.

mystran · Post by **mystran** » Sat Apr 18, 2009 7:19 am

So you are filtering a sinewave and measuring the amplitude of the output to get the graph?

I can't exactly see how that could possibly work except maybe if you used complex sinusoids but picked the complex magnitude to eliminate phase.. but at least I think it'd explain the ripple..

In any case, I'm pretty sure the ripple can't possibly be the result of a biquad, because there simply isn't enough degrees of freedom for more than one bump and more than one notch.

thevinn · Post by **thevinn** » Sat Apr 18, 2009 10:51 am

mystran wrote:So you are filtering a sinewave and measuring the amplitude of the output to get the graph?

Yes.

mystran wrote:I can't exactly see how that could possibly work except maybe if you used complex sinusoids but picked the complex magnitude to eliminate phase.. but at least I think it'd explain the ripple..

Well if it doesn't work, then it is producing graphs that are shockingly similar to what you would expect from a working version, for the case of the 3-band eq. Its only the biquad thats acting funny. 3-band eq is looking like a champ.

Others said that my approach should work. It sure seems like it should work. I mean, I took the concepts of magnitude response to different frequencies quite literally and applied them. What makes you say it shouldn't work?

mystran · Post by **mystran** » Sat Apr 18, 2009 11:09 am

thevinn wrote:
Others said that my approach should work. It sure seems like it should work. I mean, I took the concepts of magnitude response to different frequencies quite literally and applied them. What makes you say it shouldn't work?

Ah, maybe I'm not understanding what exactly is it that you do... i guess if you run the sine continuously for a while such that it's steady state and then pick a peak amplitude that'd be workable.. but you need a few cycles minimum to stabilize it, otherwise you can get ripple depending on time-domain properties of the filter (since it's not linear-phase) [edit: and even if it was linear-phase you'd need to stabilize it into steady state, but then it'd be trivial as you'd just make sure to fill all the taps from the sine]

thevinn · Post by **thevinn** » Sat Apr 18, 2009 11:18 am

mystran wrote:Ah, maybe I'm not understanding what exactly is it that you do... i guess if you run the sine continuously for a while such that it's steady state and then pick a peak amplitude that'd be workable.. but you need a few cycles minimum to stabilize it, otherwise you can get ripple depending on time-domain properties of the filter (since it's not linear-phase) [edit: and even if it was linear-phase you'd need to stabilize it into steady state, but then it'd be trivial as you'd just make sure to fill all the taps from the sine]

Its like, I understand each word individually but then when you put them together I am left scratching my head.

Time-domain properties? Linear-phase? Steady state? Taps?

Why would the filter's output be any different for the first copy of the sine wave than the Nth copy? Its a first order Biquad. Those things only remember the last two samples. So the phase can only differ by 2 samples at most (I think). The sine wave is much larger than that (2048 samples).

Alright I'm reading:

http://en.wikipedia.org/wiki/Linear_phase

mystran · Post by **mystran** » Sat Apr 18, 2009 12:59 pm

There is no first order biquad. Biquad is a second order IIR filter.

Recursive filters have state, and they are called infinite-impulse response (IIR) filters because they are (theoretically, in a world where all values have infinite precision) affected by infinitely many samples in the past.

The way IIR filters work, is basically by delaying different frequencies by different amounts, such that when the result is combined with another copy of the signal without the delay (or with different delay), some frequencies (the pass-band) will be in phase, while others (the stop band) will be out of phase and cancel. It doesn't necessarily look like that though in a structure such as biquad, as such a structure could be seen as a sort of optimized setup and whatever summing there is happens as part of the filters operation.. that gotta be the most confusing possible explanation but the important thing is that as a result you will have different phase-shifts (delays) at different frequencies, and whatever way you measure your signals should not rely on particular phase unless that exact phase is known.

What I was wondering was whether you might not have taken this into account; if you looked at a single sample of the output (where there was a peak in the input) the peak might have shifted (depending on frequency) which might give you a lot of ripple in the measurement. If you're doing this, you could run two filters instead, and feed one with sine, one with cosine (that is, 90 degrees phase difference) such that you can get the amplitude response at any particular sample as the pythagorean norm (sqrt(s^2+c^2) where s=filtered sine, c=filtered cosine). That's basically the same as feeding the filter with a complex sinusoid, but avoids having to modify the filter to work with complex numbers.

Steady-state means "once the response has stabilized" (as opposed to "transient response" which would be what you get immediately after something changes suddenly, such as an input being turned on) because when you start feeding a filter with a sine, you aren't actually feeding the filter with just the sine: the sine also has an envelope. If you just start with a sine without any fading and the filter was reset to zero state (the steady-state of a biquad in response to zero signal), this envelope will be the step function (instant attack) which has a spectra which falls off at 6dB/octave; once it's multiplied with the sine, you have a spectral peak at the sine-frequency which falls down at 6dB/octave both above and below (on frequency axis) of the center frequency (the sine frequency), and the filter will respond to this transient as well. If you let the sine run for a while, the response will stabilize into a "steady state" as the filter will gradually forget the transient... but I admit that for low-order filters you don't need a hugely long signal to stabilize the response and your 2k samples should be enough unless your cutoff frequency is very low.

Oh and "taps" was referring to a FIR filters, since linear-phase (all frequencies delayed by the same amount) implies finite impulse response if the filter is supposed to be causal (as linear-phase impulse response has to be symmetric around some point in time, and a infinite symmetric response would have to see infinitely far into the future).

thevinn · Post by **thevinn** » Sat Apr 18, 2009 1:14 pm

mystran wrote:The way IIR filters work, is basically by delaying different frequencies by different amounts, such that when the result is combined with another copy of the signal without the delay (or with different delay), some frequencies (the pass-band) will be in phase, while others (the stop band) will be out of phase and cancel.

Great explanation. I was wondering how these things worked. I read up on phase response and I get it now. I see what you mean about phase differences being used to cancel things out.

mystran wrote:What I was wondering was whether you might not have taken this into account; if you looked at a single sample of the output (where there was a peak in the input) the peak might have shifted (depending on frequency) which might give you a lot of ripple in the measurement.

Its definitely been in the back of my mind. There is a noticable discontinuity when the sine wave starts up (first sample==1.0). The filter gets re-initialized for each test frequency and that means its sample history has zeroes. Using the sine, I get -24dB falloff at the cutoff frequency. With a cosine I am getting -19dB falloff. There should not be a difference but there is. So I can confirm that my method of measurement is not independent of phase.

but I admit that for low-order filters you don't need a hugely long signal to stabilize the response and your 2k samples should be enough unless your cutoff frequency is very low.

Things start getting funky at the low frequencies. I will try only calculating the magnitude resposne of the last half of the buffer and see if that changes things.

thevinn · Post by **thevinn** » Sat Apr 18, 2009 1:24 pm

So anyway, apparently I'm not the first person with the idea of using a sine wave over and over again. This guy has an algorithm almost identical to mine:

http://www.musicdsp.org/archive.php?classid=3#251

thevinn · Post by **thevinn** » Sat Apr 18, 2009 1:37 pm

Can someone please give me a good value for Q for a low pass?

Music Engineer · Post by **Music Engineer** » Sat Apr 18, 2009 1:52 pm

mystran wrote: i guess if you run the sine continuously for a while such that it's steady state and then pick a peak amplitude that'd be workable.. but you need a few cycles minimum to stabilize it, otherwise you can get ripple depending on time-domain properties of the filter (since it's not linear-phase)

that was what i was first thinking too. but there's still a catch, even when the the filter has stabilized: simply measuring the highest absolute value of the output sinusoid will still have an error because you can't be sure that the peaks of the (input- or output) sinusiod will coincide with sample instants. that is to say: the maximum amplitude may occur in between samples as in intersample peaks.

Music Engineer · Post by **Music Engineer** » Sat Apr 18, 2009 1:55 pm

thevinn wrote:Can someone please give me a good value for Q for a low pass?

for a biquad lowpass a-la RBJ, a Q value of 1/sqrt(2)=0.707... will give the fastest transition between passband and stopband without overshooting unity gain (without a resonant bump). ...if that's good

thevinn · Post by **thevinn** » Sat Apr 18, 2009 1:57 pm

Robin from www.rs-met.com wrote:peaks of the (input- or output) sinusiod will coincide with sample instants. that is to say: the maximum amplitude may occur in between samples as in intersample peaks.

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 )

thevinn · Post by **thevinn** » Sat Apr 18, 2009 2:03 pm

I'm using the sum of the absolute value of the samples as the magnitude. I have 4800 samples in my sine buffer. The green line is the return value from getBiquadMagnitudeAt().

I am extremely confident about the x and y axis coordinates. I've stepped through the code, and I've done test cases and the mappings always put the frequencies where you would expect. In my formula I am using -6dB = half power (i.e. 20*log(a/a0)).

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). I have filter frequency = 7999 and Q=0.71:

Why is getBiquadMagnitudeAt() not giving me the expected drop near the cutoff frequency?

Why is my biquad rolling off early? Should it be?

Whats with the bump?

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

Does this look like an acceptible frequency response?

Taking another crack at basic filtering