FIR design (or, how was lopass_1000_44100.h calculated?)

Status
Not open for further replies.
Hi there,

I've been trying in vain for the last few days with various Python scripts, as well as http://t-filter.engineerjs.com/ (the page recommended in the Audio System Design Tool), to generate FIR filter coefficients for the AudioFilterFIR object.

Everything I try just produces horribly distorted and aliased audio output!

However, if I copy the filter coefficients from lopass_1000_44100.h (from File > Examples > Audio > Effects > Filter_FIR), this filter sounds really nice, with hardly any distortion, and it works as it should.

I'm wondering if the author of lopass_1000_44100.h is around, would they be able to describe here what tool they used to calculate these coefficients, and how they did it and thought about it more generally?
 
I have no idea which tool Paul Stoffregen used, but the sources include screen captures.

The FIR filter coefficients are in signed fixed-point 16-bit format. The actual coefficient values are 21 17 23 30 37 44 50 55 59 60 59 54 46 33 15 -7 -34 -65 -100 -138 -177 -217 -254 -289 -318 -340 -352 -353 -340 -312 -268 -207 -129 -33 80 208 350 504 666 833 1002 1169 1330 1482 1619 1740 1841 1919 1972 1999 1999 1972 1919 1841 1740 1619 1482 1330 1169 1002 833 666 504 350 208 80 -33 -129 -207 -268 -312 -340 -353 -352 -340 -318 -289 -254 -217 -177 -138 -100 -65 -34 -7 15 33 46 54 59 60 59 55 50 44 37 30 23 17 21, and it has about 63dB of attenuation at higher frequencies. (You can paste those into my FIR Analysis page, which is a trivial FIR filter visualizer in a single Public Domain HTML file, to see its frequency response. Since I created that page just as an illustration on how we can easily create tool web pages, it does not have e.g. any frequency scale on the horizontal axis; it is a demonstration of the web-page-as-a-tool approach only. The number of points determines the fidelity of the graph; I suggest 5000 for this filter so you see all relevant detail.)

Let's say you use whatever tool to generate a set of FIR filter coefficients as decimal numbers. Multiply the coefficients by 32768 (=215) to scale them to 16-bit first, then add 65536 (=216) to each negative value (converting the signed to unsigned using two's complement format, as used here in the filter).

There are many tools for generating FIR filter coefficients. One option is the fir1() generator function in Octave (which is free and available for basically all operating systems), especially because many Matlab FIR filter guides also work in Octave, and you can easily do the scaling and conversion within Octave:
Code:
pkg load signal
filter = fir1(99, 1700/44100, "low");
yields a 99-tap (100 coefficient) lowpass filter with cutoff frequency at 1700 Hz when using 44100 samples per second.
To plot the frequency response as a graph, use
Code:
coeff = round(filter*32768) / 32768;
plot((-0.5:1/4096:0.5-1/4096)*44100,20*log10(abs(fftshift(fft(coeff,4096))))): axis([0 22050 -80 10]);
Finally, use
Code:
mod(round(filter*32768 + 65536), 65536)
to get the unsigned 16-bit coefficients you need to populate the array with.
 
Hi there,

I've been trying in vain for the last few days with various Python scripts, as well as http://t-filter.engineerjs.com/ (the page recommended in the Audio System Design Tool), to generate FIR filter coefficients for the AudioFilterFIR object.

Everything I try just produces horribly distorted and aliased audio output!
Clipping?

Did you normalize the coefficients and scale them to signed 16 bit? If your filter has any gain you need to ensure you reduce
the values to prevent clipping for signals at the peak frequency, or you'll lose headroom.

Normally for a low pass filter (without any peaking) you can just scale so that the _sum_ of the coefficients is close to 32767
(but not greater than 32767), ensuring unity DC gain. Perhaps you scaled the maximum coefficient to 32767 - that's not right.
 
Thanks for the info guys!

Perhaps you scaled the maximum coefficient to 32767 - that's not right.

I did do this initially. But I also tried some much smaller scales (e.g coefficient*0xF) and they all distorted as well. But yes it does seem like clipping, so I will definitely try your approach of normalisation in my Python script and see if that works.

Multiply the coefficients by 32768 (=215) to scale them to 16-bit first, then add 65536 (=216) to each negative value (converting the signed to unsigned using two's complement format, as used here in the filter).
Finally, use
Code:
mod(round(filter*32768 + 65536), 65536)
to get the unsigned 16-bit coefficients you need to populate the array with.

Nominal Animal, are you able to explain the conversion to unsigned to me more fully? This definitely looks like it could be the difference between my coefficients (I just kept any negative short ints) vs the ones from lopass_1000_44100.h which are all hex.

Here is the Python script I wrote:
Code:
# Python script to generate FIR coefficients
from scipy import signal
numtaps = 150
f = 300
coeffs = signal.firwin(numtaps, f, pass_zero='lowpass', fs=44117)
# scale to 16 bit integers
coeffs = [int(x * 0x7fff) for x in coeffs]
# print C array for copy/paste
print("// Low pass filter, cutoff frequency " + str(f) +"Hz.")
print("#define n_coeffs " + str(numtaps))
print("const short coeff_p[n_coeffs] = {", end="")
lineLength = 15
strings = [str(x) for x in coeffs]
print(",\n                                 ".join([", ".join(strings[i:i+lineLength]) for i in range(0,len(strings),lineLength)]), end="")
print("};")
 
Nominal Animal, are you able to explain the conversion to unsigned to me more fully?
Sure.

Why do we multiply by 32768? Because the function that applies the FIR filter is arm_fir_fast_q15(), and that uses Q15 fixed point format, which has a sign bit and 15 fractional bits.
The most negative value a Q15 fixed point number can describe is -1.0 (corresponding to int16_t -32768).
The most positive value a Q15 fixed point number can describe is 0.99996948 (corresponding to int16_t 32767).

We then round to nearest integer. We could also wrap max(-32768, min(32767, ...))) around the expression, to clamp the coefficients to the valid range, but usually it is not necessary.

For some reason, the coefficients in the array are specified as (short)0xHHHH where the 0xHHHH is the coefficient in unsigned 16-bit integer format. ARM (like most desktop computers) use two's complement format to represent negative numbers, so to convert to unsigned, we just need to add 65536 to each negative value. Since we know the values are already -32768 to 32767, inclusive, we can just add 65536 to all values (not just the negative ones) so that they are all positive (or zero); and keep just the least significant 16 bits. We can do that either via binary AND (Python: & 65535) or remainder aka modulo operation (Python: % 65536); I used the Octave modulo/remainder operation.

If you have Python array coeffs, you can use
Code:
values = [ ((max(-32768, min(32767, round(32768.0*c))) + 65536) & 65535) for c in coeffs ]
to generate the array of values, or for example
Code:
for c in coeffs:
    print("\t(short)0x%04x," % ((max(-32768, min(32767, round(32768.0*c))) + 65536) & 65535))
to print the coefficients.

I'm not exactly sure why this should matter, but I guess there is some kind of funky integer promotion stuff going on. I'm not very familiar with this code, since I don't do audio with my Teensys.
 
I generated those coefficients using a program that I wrote many years ago (originally for a Commodore Amiga) by translating the Parks-McLellan filter design program from the original FORTRAN into C. Some years after that I ported it to a PC and added a graphical front end.
Here's a screenshot of the program showing the filter's design parameters and the resulting response.
lopass_1000_cd_myfir.gif
I used the program with several different DSP processors starting with a TI one back around 1995 (I think). I then used an ADSP EZKIT Lite. Each of them required the filter coefficients to be specified in a specific form and I either modified the code to produce that specific format or edited the output. I don't remember why I used the hex format in the Teensy example but "(short)0x0015" is equivalent to "21" and "(short)0xFFF9" is equivalent to "-7", etc.

FYI: A slightly updated form of the original FORTRAN code is here and this page has C code for FIR and IIR filter design.

Pete
 
The scipy.signal library (*), coupled with matplotlib, makes all this filter design and evaluation
relatively painless in Python. Perhaps its time to move on from Fortran!

(*) Gives the similar functionality to MATLAB without having to buy it or learn its scripting language.
 
Status
Not open for further replies.
Back
Top