Forum Rule: Always post complete source code & details to reproduce any issue!
Results 1 to 7 of 7

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

  1. #1
    Junior Member
    Join Date
    Oct 2020
    Location
    Melbourne, Australia
    Posts
    12

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

    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?

  2. #2
    Senior Member
    Join Date
    Feb 2015
    Location
    Finland
    Posts
    229
    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.

  3. #3
    Senior Member
    Join Date
    Jul 2020
    Posts
    1,118
    Quote Originally Posted by JeremiahRose View Post
    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.

  4. #4
    Junior Member
    Join Date
    Oct 2020
    Location
    Melbourne, Australia
    Posts
    12
    Thanks for the info guys!

    Quote Originally Posted by MarkT View Post
    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.

    Quote Originally Posted by Nominal Animal View Post
    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).
    Quote Originally Posted by Nominal Animal View Post
    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("};")

  5. #5
    Senior Member
    Join Date
    Feb 2015
    Location
    Finland
    Posts
    229
    Quote Originally Posted by JeremiahRose View Post
    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.

  6. #6
    Senior Member
    Join Date
    Nov 2012
    Posts
    1,693
    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.
    Click image for larger version. 

Name:	lopass_1000_cd_myfir.gif 
Views:	18 
Size:	8.7 KB 
ID:	25035
    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

  7. #7
    Senior Member
    Join Date
    Jul 2020
    Posts
    1,118
    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.

Posting Permissions

  • You may not post new threads
  • You may not post replies
  • You may not post attachments
  • You may not edit your posts
  •