tttapa / Arduino-Filters

Arduino Finite Impulse Response and Infinite Impulse Response filter implementations.
GNU General Public License v3.0
115 stars 15 forks source link

Extract coefficients from IIR filter #4

Open sylvanoMTL opened 3 years ago

sylvanoMTL commented 3 years ago

Hi there, I see the implementation require high knowledge of C++ coding. I have been through the example of the Butterworth filter, and I am keen on getting the b and a coefficients to compare with an algorithm I have developed on another platform. Do you know a simple method to extract them?

Also can you confirm the filter is only specified by the cut-off frequency? (i.e. no pass band ripple, stop band attenuation?

Thank you Sylvano

tttapa commented 3 years ago

The butter_coeff function returns the coefficients as a list of second-order sections (SOS). You can then use the sos2tf function to convert these SOS coefficients to a rational transfer function, see the TransferFunction class, it has b and a members you can access.

For example:

#include <Filters.h>

#include <Filters/Butterworth.hpp>

// Sampling frequency
const double f_s = 100; // Hz
// Cut-off frequency (-3 dB)
const double f_c = 25; // Hz
// Normalized cut-off frequency
const double f_n = 2 * f_c / f_s;

// Sixth-order Butterworth filter
auto butter_sos = butter_coeff<6>(f_n);
auto butter_tf = sos2tf(butter_sos);

void setup() {
  Serial.begin(115200);
  while (!Serial);
  Serial << setprecision(7);
  Serial << "fₙ: " << f_n << endl;
  Serial << "b:  " << butter_tf.b << endl;
  Serial << "a:  " << butter_tf.a << endl;
}

void loop() {}

Prints:

fₙ: 0.5000000
b:  0.0295882, 0.1775293, 0.4438233, 0.5917644, 0.4438233, 0.1775293, 0.0295882
a:  1.0000000, 0.0000000, 0.7776959, 0.0000000, 0.1141994, 0.0000000, 0.0017509

You can compare them to the coefficients generated by SciPy, for example:

from scipy.signal import butter
f_s = 100 # Hz, Sampling frequency
f_c = 25 # Hz, Cut-off frequency (-3 dB)
f_n = 2 * f_c / f_s # Normalized cut-off frequency
b, a = butter(6, f_n)
print('b:', ', '.join(map(lambda x: f'{x:10.7f}', b)))
print('a:', ', '.join(map(lambda x: f'{x:10.7f}', a)))

Prints:

b:  0.0295882,  0.1775293,  0.4438234,  0.5917645,  0.4438234,  0.1775293,  0.0295882
a:  1.0000000, -0.0000000,  0.7776960, -0.0000000,  0.1141994, -0.0000000,  0.0017509

A Butterworth filter has no ripple, apart from the filter order, the cut-off frequency is the only parameter, see https://tttapa.github.io/Pages/Mathematics/Systems-and-Control-Theory/Analog-Filters/Butterworth-Filters.html#mjx-eqn-8.

sylvanoMTL commented 3 years ago

Great thank you for your prompt reply and the additional precisions. I have been confused by the MATLAB designfilt GUI. But it seems all make sense to me now. regards Sylvain