Forums

Problem with fixed point lowpass filter design

Started by Malcolm Haylock August 16, 2003
Hi Everyone,

I'm a newcomer to DSP and am trying to write a Fixed Point 
implementation of the bilinear transform lowpass filter as outlined in 
the Audio EQ Cookbook 
(http://www.harmony-central.com/Computer/Programming/Audio-EQ-Cookbook.txt):

     y[n] = (b0/a0)*x[n] + (b1/a0)*x[n-1] + (b2/a0)*x[n-2]  - 
(a1/a0)*y[n-1] - (a2/a0)*y[n-2]

where for a lowpass filter:

     b0 =  (1 - cos)/2
     b1 =   1 - cos
     b2 =  (1 - cos)/2
     a0 =   1 + alpha
     a1 =  -2*cos
     a2 =   1 - alpha

where :
     cos   = cos(2*pi*frequency/sampleRate)
     alpha = sin/(2*Q)

The trouble is that I can't find a way to get the coefficients within 
the range of my fixed point calculations. For 44.1kHz sample rate things 
are fine with cutoff frequencies down to about 400Hz but below that 
things get tricky.

e.g. for samplerate=44100Hz, Q=1.0, frequency=100Hz:

b0/a0 = 5.03944E-05
b1/a0 = 1.00789E-04
b2/a0 = 5.03944E-05
a1/a0 = -1.98565E+00
a2/a0 = 9.85854E-01

I am working with 24-bit samples and need to keep everything within 
32-bits. Therefore I can afford to multiply the above coefficients by 
256 for fixed point calculations. I need to have coefficients within the 
range 1/256 to 256. In the above example the first three coefficients 
will all be 0 (even after multiplying by 256) so I will not get any 
sound out.

Am I being realistic and is there a way around this using the above 
design or should I be using a different design altogether?

Thanks very much,
Malcolm Haylock
smaugNOSPAM@kagi.com (remove NOSPAM)

Malcolm Haylock wrote:
> Hi Everyone, > > I'm a newcomer to DSP and am trying to write a Fixed Point > implementation of the bilinear transform lowpass filter as outlined in > the Audio EQ Cookbook > (http://www.harmony-central.com/Computer/Programming/Audio-EQ-Cookbook.txt): > > > y[n] = (b0/a0)*x[n] + (b1/a0)*x[n-1] + (b2/a0)*x[n-2] - > (a1/a0)*y[n-1] - (a2/a0)*y[n-2] > > where for a lowpass filter: > > b0 = (1 - cos)/2 > b1 = 1 - cos > b2 = (1 - cos)/2 > a0 = 1 + alpha > a1 = -2*cos > a2 = 1 - alpha > > where : > cos = cos(2*pi*frequency/sampleRate) > alpha = sin/(2*Q) > > The trouble is that I can't find a way to get the coefficients within > the range of my fixed point calculations. For 44.1kHz sample rate things > are fine with cutoff frequencies down to about 400Hz but below that > things get tricky. > > e.g. for samplerate=44100Hz, Q=1.0, frequency=100Hz: > > b0/a0 = 5.03944E-05 > b1/a0 = 1.00789E-04 > b2/a0 = 5.03944E-05 > a1/a0 = -1.98565E+00 > a2/a0 = 9.85854E-01 > > I am working with 24-bit samples and need to keep everything within > 32-bits. Therefore I can afford to multiply the above coefficients by > 256 for fixed point calculations. I need to have coefficients within the > range 1/256 to 256. In the above example the first three coefficients > will all be 0 (even after multiplying by 256) so I will not get any > sound out. > > Am I being realistic and is there a way around this using the above > design or should I be using a different design altogether?
If you have 32 bits then you have plenty of room for scaling. One thing to do is to scale the b coef's as 1, 2, 1, do the sum for x terms, then multiply by 1/a0 to add in the y terms. The other thing to note is that you're assuming the y terms are going to be pretty small, so work with a0*y[n] = b0x[n] + b1x[n-1] + b2x[n-2] + a1/a0*(a0*y[n-1]) + a2/a0*(a0*y[n-2]) This puts all your coefficients in the same range as well as your results. Your final output has to be scaled back by a0. To make life easy, change a0 in my formula to 2^n, and then you can scale by a simple shift. Adjust the coefficients accordingly too! Patience, persistence, truth, Dr. mike -- Mike Rosing www.beastrider.com BeastRider, LLC SHARC debug tools
> b0/a0 = 5.03944E-05 > b1/a0 = 1.00789E-04
Hi,
> b2/a0 = 5.03944E-05
will be 1 i.e. scaler=round(1/5.0E-05) to 2^n. If overflow will be by adding then scaler :2, but values b2/a0 rest 1 if you afraid 0. You don't see difference after division on scaler(>>n). Error of round up isn't reason for absent of sound. Cheers
> > y[n] = (b0/a0)*x[n] + (b1/a0)*x[n-1] + (b2/a0)*x[n-2] - > (a1/a0)*y[n-1] - (a2/a0)*y[n-2] > > > e.g. for samplerate=44100Hz, Q=1.0, frequency=100Hz: > > b0/a0 = 5.03944E-05 > b1/a0 = 1.00789E-04 > b2/a0 = 5.03944E-05 > a1/a0 = -1.98565E+00 > a2/a0 = 9.85854E-01 > > I am working with 24-bit samples and need to keep everything within > 32-bits. Therefore I can afford to multiply the above coefficients by > 256 for fixed point calculations. I need to have coefficients within the > range 1/256 to 256. In the above example the first three coefficients > will all be 0 (even after multiplying by 256) so I will not get any > sound out. >
[b0 b1 b2 a1 a2]/a0*256 = [0.0129 0.0258 0.0129 -508.3264 252.3786] the first 3 terms are not 0. In fixed point, once you chose your fixpt format, in this case, 24bit, should be able to expressed in Q10.14. You said, to have the coefficients in 1/256~256, should multiply 128, [b0 b1 b2 a1 a2]/a0*128 = [0.0065 0.0129 0.0065 -254.1632 126.1893] Why set range on the coefficients, you may scale the singal and leave some space for accurate coefficients. Besides, you may try some multirate filter in this design. Lin. S.