IIR filter coefficients generation

[apar945]

Hello,

I tried to compare the filter coefficients generated by the Butterworth filter class of the library. The values of poles and zeros match to the output from this reference link https://www-users.cs.york.ac.uk/~fisher/mkfilter/trad.html But the numerators, denominators and gain are completely different. What could I be doing wrong when using your library?

regards Ankit

[ar1st0crat]

Hi!

I think, it's just how the results are interpreted. The output generated by the website has the form which is slightly different from NWaves/MATLAB/sciPy functions. I've specified the following parameters:

filtertype = Butterworth
passtype = Lowpass
order = 5
samplerate = 16000
corner1 = 1000

The results are:

Recurrence relation:
y[n] = (  1 * x[n- 5])
     + (  5 * x[n- 4])
     + ( 10 * x[n- 3])
     + ( 10 * x[n- 2])
     + (  5 * x[n- 1])
     + (  1 * x[n- 0])

     + (  0.2777529978 * y[n- 5])
     + ( -1.7411025201 * y[n- 4])
     + (  4.4205122516 * y[n- 3])
     + ( -5.6938879540 * y[n- 2])
     + (  3.7314736649 * y[n- 1])

gain at dc    :   mag = 6.093427809e+03   phase =   0.0000000000 pi

Now, the same thing in in NWaves:

double fc = 1000;
double fs = 16000;
var lpf = new LowPassFilter(fc / fs, 5);

var gain = lpf.Tf.Gain;
var num = lpf.Tf.Numerator; 
var num = lpf.Tf.Denominator;
// ...

And the results are:

gain = 0.000164111241045019

Denominator = [
1
-3.73147366494481
5.6938879539761
-4.42051225162662
1.74110252012946
-0.277752997820681
]

Numerator = [
0.000164111241045019
0.000820556205225094
0.00164111241045019
0.00164111241045019
0.000820556205225094
0.000164111241045019
]

So, denominator coeffs are identical (if the difference equation is written in the form as used on this website, then the sign of coefficients in recursive part of recurrence relation should be changed to opposite).

As for numerator coeffs, they are correct as well. In NWaves (just like in MATLAB and in sciPy) gain is "reflected" in TF numerator coeffs. So, if you divide all coeffs by gain (in this case 0.000164111241045019), you will get [1 5 10 10 5 1].

Finally, the gain is also correct:

gain at dc    :   mag = 6.093427809e+03

// note, in C code on the website the coeffs are divided by this gain,
// i.e. multiplied by 1/gain:

1/gain = 0.00016411124

PS. In MATLAB, for example:

[b,a] = butter(5, 2*1000/16000);
[z,p,k] = tf2zpk(b,a)

you'll get the same results as in NWaves (I wanted this lib to be compliant with standard DSP tools).

[apar945]

Hello Tim,

Thank you for taking the time to explain this in detail. Thank you for sharing this code library. It is very very helpful. Just one more question - Why the does code expect an odd number for the order of the filter?

regards Ankit

[apar945]

Hello Tim,

I tried another comparison based on your review above, but for a Butterworth BandPass filter. Here are the parameters filtertype = Butterworth passtype =Bandpass  order = 9 samplerate = 250 corner1 = 8.1 corner2 =12.9

gain at centre: mag = 1.955889078e+09 phase = -0.9229604531 pi

The gain generated by the NWaves library for the same configuration is Gain 1.2536463403844429E-09

1/gain at centre is not equal to the value above.

This is where things fall apart for me. Just wanted to confirm that this is not a bug and has an explanation.

Regards Ankit

[ar1st0crat]

Hi Ankit,

The odd order is required only for FIR filters. In case of IIR filters (including the Butterworth filter) you can specify any number. In the future I'm going to add code for designing even-order FIR filters as well.

As for the gain of your bandpass filter, I guess in this particular case some round-off errors take place. The numerator coeffs are very close to zero, gain is very small number.:

gain =  0.0000000019938  // NWaves
gain = -0.0000000033447  // MATLAB
gain = 0.000000000007548 // sciPy

Anyway, the results are consistent - if you divide by gain, you'll get correct numerator coeffs:

var bpf = new BandPassFilter(8.1/250, 12.9/250, 9);

var gain = bpf.Tf.Gain;
Console.WriteLine("Gain: " + gain);

foreach (var coeff in bpf.Tf.Numerator)
{
    Console.WriteLine(coeff / gain);
}

// Gain: 1.99388321513846E-09
// 1 0 -9 0 36 0 -84 0 126 0 -126 0 84 0 -36 0 9 0 -1

Regards Tim

[apar945]

Hello Tim,

Thank you for getting back to me. The gain you show above (// Gain: 1.99388321513846E-09) does not match the gain value generated by the code I am using. I am referring to the NWaves.DemoForms project. In the filterform.cs line 480, there is a function called AnalyzeButterworthFilter. I added the following console output to match yours and the gain is not as generated by your code.

filter = new Filters.Butterworth.BandPassFilter(8.1 / 250, 12.9 / 250, 9);
            var gain = _filter.Tf.Gain;
            Console.WriteLine("Gain: " + gain);
            foreach (var coeff in _filter.Tf.Numerator)
            {
                Console.Write(coeff / _filter.Tf.Gain );
                Console.Write("\t");
            }
            Console.Write("\n");

//Gain: 1.25364634038444E-09
//1 0 -9 0  36 0    -84 0 126   0 -126 0 84 0 -36 0 9 0 -1

I believe this is the root cause of my filter not working. The output is not scaled by this gain to reflect the input waveform. Anyway I understand the round off errors are an issue as I am dealing with micro volts signal and will try to troubleshoot this issue.

Regards Ankit

[ar1st0crat]

Hi Ankit!

Take a look at this issue. I think, you're facing the same problem.

Regards Tim

[apar945]

Hello Tim,

Thank you for that suggestion, that helped partially. I tested going up to 14th order of filters and things work for a band pass with a bandwidth of more than 10 Hz. But for smaller bands of 4 Hz bandwidth, the output saturates. I noticed that with TfToSsos, the gain is assigned to the first TF in the array and the rest of the TFs have a gain 1. Any other suggestions to tackle the issue with the small bandwidth bands?

Is cascading an option?

regards Ankit

[ar1st0crat]

Hi Ankit,

Cascading is an option (I mean SOS cascades, which are always better in terms of numeric computations). However, if even second-order-sections are unstable, then the bandwidth is just too small for Butterworth filter design. I quickly checked sciPy results - they are problematic too.

Are FIR filters too expensive for your task?

Regards Tim

[apar945]

Hello Tim,

I am trying to work with in real time. I will experiment with the FIR filters to see if they can keep up with the application. The purpose of filtering is to find the power spectrum of the input signal. FFT is an option but requires a second of data as the sampling rate is not as high as for an audio signal and 1 second is already not real time.

Are there any other methods of calculating the power spectral density in real time that can experimented with?

Regards Ankit

[ar1st0crat]

Hi Ankit,

For achieving 4Hz frequency resolution with signal sampled at 250 Hz, one second of data is required, regardless of the filtering technique. Since there's no way, in general, to pre-fetch future samples, it's theoretically impossible not to be "late" by 250 samples.

Regards Tim

[apar945]

Hello Tim,

Apologies for the delay in the response. I will revisit in a few weeks and let you know what I find.

regards Ankit

[ar1st0crat]

Hi Ankit,

Maybe you'll find it helpful:

33

Regards, Tim

[apar945]

Hello Tim,

Thank you for sharing the ticket. I will try the suggestion and report back.

Regards Ankit