Search code examples
pythonfilterscipyintegerbutterworth

Scipy filter : force minimal value of SOS coefficient to prepare integer filter


I try to design a low-pass filter (cutting around 40Hz) with Scipy, which seems to work with the following code:

fs = 100000 # Sampling frequency = 100kHz

N, Wn = signal.buttord(40/(fs/2), 50/(fs/2), 0.1, 5)
sos = signal.butter(N, Wn, 'low', output='sos')

This filter must eventually be embedded on a board without FPU so I must convert the SOS into an array of integers.

But here is the issue: some coefficients of the first line of the SOS matrix are too low to be converted into 32 bits integers:

[ 8.00108536e-32  1.60021707e-31  8.00108536e-32  1.00000000e+00
  -9.97022785e-01  0.00000000e+00]

Replacing the almost null value by 0 or 1 (once all other coefficients scaled up) did not work.

Do you know a scipy utility that forces the design of a filter to be compatible with a conversion to integer ?

If not, do you know how can I change the buttord or butter args to increase the coefficients of the first line ? (the sampling freq and the cut frequency cannot be changed).

Thanks in advance


Solution

  • I finally solved it with a post-treatment. Back to the definition of the SOS, the three first columns are the coefficients of the numerator of cascaded biquad filters.

    Therefore the negative powers of the coefficient of the first line can be spread over all the lines. For instance as follow:

    power = int(np.log10(np.abs(sos[0, 0:3]).mean())-1)
    power_by_line = -power // len(sos)
    scale = 10 ** (-power_by_line)
    
    for i in range(1, len(sos)):
        sos[i, 0:3] = sos[i, 0:3] * scale
    
    remaining = (len(sos)-1)*power_by_line
    sos[0, 0:3] = sos[0, 0:3] * 10**remaining
    

    The resulting filter seems to have the same behavior (for the numerical precision that I require).