2
$\begingroup$

I am trying to implement a digital filter over a uC (it doesn't really matter which filter and which micro controller because I'm looking forward to learn how to do it in the future with different filters and different microcontrollers). I've been told that you can design, implement and debug a digital filter in python and when everything is ready you can port the code to C without changing anything. How can I do that? I've been searching for a while and I can not find how to do this.

Also, is there a way to plot the transfer function of an implemented fiter?. I mean the actual filter, a function made by me, that takes the input values and performs the calculations. I don't want to plot the transfer function of a filter made with functions like scipy.signal.butter. I want to plot the transfer function of a filter made with a for, some multiplications and sums.

I would really appreciate any help or information that you can provide me.

$\endgroup$
  • 2
    $\begingroup$ well, you gotta know python. i s'pose you can do the same with MATLAB of Octave. i think it might be easier to just write your digital filter in C. i can give you some generic C code if you want (send me an email). for your uC to compute coefficients might be a bitch. is your digital filter one that will change during run time? or do you know all of the parameters of the digital filter at compile time? $\endgroup$ – robert bristow-johnson Jul 22 '19 at 21:34
  • $\begingroup$ Please have a look at our article on how to export filters to Python. The tool itself should help you. $\endgroup$ – ASN Apr 21 at 16:14
3
$\begingroup$

Filter representation and design

A DTLTI IIR filter is characterized by its transfer function $ H(z) = \frac{Y(z)}{X(z)} = \frac{b_0 + b_1 z^{-1} + b_2 z^{-2} + \dots + b_{P} z^{-P}}{a_0 + a_1 z^{-1} + a_2 z^{-2} + \dots + a_{Q} z^{-Q}} $. The transfer function is closely related to the difference equation: $$ y[n] = \frac{1}{a_0} \left( \sum_{i=0}^{P} b_i x[n-i] - \sum_{j=1}^{Q} a_j y[n-j] \right) $$ ($x[n]$ is the input at time step $n$, $y[n]$ is the output at time step $n$).
This is equivalent to $$ \sum_{i=0}^{P} b_i x[n-i] = \sum_{j=0}^{Q} a_j y[n-j] $$ By applying the Z transform to both sides, using the time shift property, and rearranging some factors, you arrive at the formula for the transfer function $H(z)$ mentioned previously.

When you design a Butterworth filter (or any IIR filter) using SciPy, it returns the $b_i$ and $a_j$ coefficients. These coefficients determine the transfer function and therefore the frequency response of the filter. You can use the freqz function to calculate this response (essentially, it evaluates $H(z)$ along the unit circle $z=e^{j\omega}$).

For example, using Python:

from scipy.signal import butter, freqz
import matplotlib.pyplot as plt
from math import pi
import numpy as np

f_s = 360    # Sample frequency in Hz
f_c = 45     # Cut-off frequency in Hz
order = 4    # Order of the butterworth filter

omega_c = 2 * pi * f_c       # Cut-off angular frequency
omega_c_d = omega_c / f_s    # Normalized cut-off frequency (digital)

# Design the digital Butterworth filter
b, a = butter(order, omega_c_d / pi)    
print('Coefficients')
print("b =", b)                           # Print the coefficients
print("a =", a)

w, H = freqz(b, a, 4096)                  # Calculate the frequency response
w *= f_s / (2 * pi)                       # Convert from rad/sample to Hz

# Plot the amplitude response
plt.subplot(2, 1, 1)            
plt.suptitle('Bode Plot')
H_dB = 20 * np.log10(abs(H))              # Convert modulus of H to dB
plt.plot(w, H_dB)
plt.ylabel('Magnitude [dB]')
plt.xlim(0, f_s / 2)
plt.ylim(-80, 6)
plt.axvline(f_c, color='red')
plt.axhline(-3, linewidth=0.8, color='black', linestyle=':')

# Plot the phase response
plt.subplot(2, 1, 2)
phi = np.angle(H)                         # Argument of H
phi = np.unwrap(phi)                      # Remove discontinuities 
phi *= 180 / pi                           # and convert to degrees
plt.plot(w, phi)
plt.xlabel('Frequency [Hz]')
plt.ylabel('Phase [°]')
plt.xlim(0, f_s / 2)
plt.ylim(-360, 0)
plt.yticks([-360, -270, -180, -90, 0])
plt.axvline(f_c, color='red')

plt.show()

The coefficients can be calculated manually, as explained here (do note that the indices of the coefficients are flipped compared to the formulas above), but it's much easier to use filter design tools like SciPy to calculate them.

Implementation

The difference equation can be used directly to implement the filter. Just loop over the previous inputs and outputs, multiply everything with the respective coefficients, and sum all terms.

This means that once you have designed your filter in Python, you just need to copy the $b$ and $a$ coefficients to your microcontroller to use the filter.

A possible C++ implementation could be:

class IIRFilter {
  public:
    template <size_t B, size_t A>
    IIRFilter(const double (&b)[B], const double (&_a)[A]) : lenB(B), lenA(A-1) {
      x = new double[lenB]();
      y = new double[lenA]();
      coeff_b = new double[2*lenB-1];
      coeff_a = new double[2*lenA-1];
      double a0 = _a[0];
      const double *a = &_a[1];
      for (uint8_t i = 0; i < 2*lenB-1; i++) {
        coeff_b[i] = b[(2*lenB - 1 - i) % lenB] / a0;
      }
      for (uint8_t i = 0; i < 2*lenA-1; i++) {
        coeff_a[i] = a[(2*lenA - 2 - i) % lenA] / a0;
      }
    }
    ~IIRFilter() {
      delete[] x;
      delete[] y;
      delete[] coeff_a;
      delete[] coeff_b;
    }
    double filter(double value) {
      x[i_b] = value;
      double b_terms = 0;
      double *b_shift = &coeff_b[lenB - i_b - 1];
      for (uint8_t i = 0; i < lenB; i++) {
        b_terms += x[i] * b_shift[i];
      }
      double a_terms = 0;
      double *a_shift = &coeff_a[lenA - i_a - 1];
      for (uint8_t i = 0; i < lenA; i++) {
        a_terms += y[i] * a_shift[i];
      }
      double filtered = b_terms - a_terms;
      y[i_a] = filtered;
      i_b++;
      if(i_b == lenB)
        i_b = 0;
      i_a++;
      if(i_a == lenA)
        i_a = 0;
      return filtered;
    }
  private:
    const uint8_t lenB, lenA;
    uint8_t i_b = 0, i_a = 0;
    double *x;
    double *y;
    double *coeff_b;
    double *coeff_a;
};

Circular buffers are used to keep track of the previous inputs and outputs. The coefficients are duplicated in a circular fashion as well, in order to simplify the indices in the loop.
If you're not allowed to use dynamic memory, you could easily create a generic class with stack allocated arrays.

For higher order filters, numerical issues may arise, as a small rounding/quantization error on one of the coefficients affects the location of all of the zeros/poles of the filter. This is especially problematic for poles that lie very close to the unit circle.
A solution is to factor the transfer function into a product of second order sections. The entire filter is then implemented as a cascade of biquad filters, each implementing one of the sections.
SciPy can produce the coefficients of these sections as well, just pass the output='sos' parameter to the butter function.

Checking the result of the implementation

Numerically calculating the frequency response from a given filter implementation is not straightforward. However, you can compare the impulse response of your implementation with the impulse response produced by SciPy (using scipy.signal.dimpulse).
The transfer function is the Z transform of the impulse response, so if the impulse responses match, the frequency characteristics will match as well.

| improve this answer | |
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.