<jupyter_start><jupyter_text>Design and test a Butterworth lowpass filterThis document describes how to design a Butterworth lowpass filter with a cutoff frequency $\omega_c$ and compute the discrete coefficients so that it can be implemented on hardware.<jupyter_code># Packages and adjustments to the figures
from scipy import signal
import matplotlib.pyplot as plt
import numpy as np
import math
plt.rcParams["figure.figsize"] = 10,5
plt.rcParams["font.size"] = 16
plt.rcParams.update({"text.usetex": True,"font.family": "sans-serif","font.sans-serif": ["Helvetica"]})<jupyter_output><empty_output><jupyter_text>1. Generate a test signal * A simple test signal $\boldsymbol{y} = \{ y_i \}$ is generated with a fixed sampling frequency using the function:$$y(t) = m_0 \sin(2\pi f_0 t) + m_1 \sin(2\pi f_1 t)$$* The power spectrum is plotted as the magnitude of the discrete fourier transform (DFT): $|\hat{\boldsymbol{y}}|$<jupyter_code># Generate a signal
samplingFreq = 1000; # sampled at 1 kHz = 1000 samples / second
tlims = [0,1]        # in seconds
signalFreq = [2,50]; # Cycles / second
signalMag = [1,0.2]; # magnitude of each sine
t = np.linspace(tlims[0],tlims[1],(tlims[1]-tlims[0])*samplingFreq)
y = signalMag[0]*np.sin(2*math.pi*signalFreq[0]*t) + signalMag[1]*np.sin(2*math.pi*signalFreq[1]*t)

# Compute the Fourier transform
yhat = np.fft.fft(y);
fcycles = np.fft.fftfreq(len(t),d=1.0/samplingFreq); # the frequencies in cycles/s

# Plot the signal
plt.figure()
plt.plot(t,y);
plt.ylabel("$y(t)$");
plt.xlabel("$t$ (s)");
plt.xlim([min(t),max(t)]);

# Plot the power spectrum
plt.figure()
plt.plot(fcycles,np.absolute(yhat));
plt.xlim([-100,100]);
plt.xlabel("$\omega$ (cycles/s)");
plt.ylabel("$|\hat{y}|$");<jupyter_output><empty_output><jupyter_text>2. Butterworth low-pass filter transfer functionThis document does not derive the formula for a Butterworth filter. Instead, it uses the standard form with DC gain $G=1$.* A cutoff frequency $\omega_c$ is selected* The Butterworth low-pass filter transfer function with $\omega_c = 1$ can be written as (see https://en.wikipedia.org/wiki/Butterworth_filter)$$H(s) = \frac{1}{\sum_1^{n} a_k s^k}$$where $n$ is the order of the filter. The coefficients are given by the recursion formula:$$a_{k+1} = \frac{\cos( k \gamma )}{\sin((k+1)\gamma)}$$with $a_0 = 1$ and $\gamma = \frac{\pi}{2n}$.* Because the Butterworth polynomial is $$B_n(s) = \sum_{k=0}^n a_k s^k$$and we want to set a new cutoff frequency of $\omega_c$, substitute$$B_n = \sum_{k=0}^n a_k \left(\frac{s}{\omega_c}\right)^k = \sum_{k=0}^n \frac{a_k}{{\omega_c}^k} s^k$$for convenience set $$B_n(s) = \sum_{k=0}^n c_k s^k$$ with $c_k = \frac{a_k}{{\omega_c}^k}$<jupyter_code># Butterworth filter
wc = 2*np.pi*5; # cutoff frequency (rad/s)
n = 2; # Filter order

# Compute the Butterworth filter coefficents
a = np.zeros(n+1);
gamma = np.pi/(2.0*n);
a[0] = 1; # first coef is always 1
for k in range(0,n):
    rfac = np.cos(k*gamma)/np.sin((k+1)*gamma);
    a[k+1] = rfac*a[k]; # Other coefficients by recursion

print("Butterworth polynomial coefficients a_i:                " + str(a))

# Adjust the cutoff frequency
c = np.zeros(n+1);
for k in range(0,n+1):
    c[n-k] = a[k]/pow(wc,k)

print("Butterworth coefficients with frequency adjustment c_i: " + str(c))
# Low-pass filter
w0 = 2*np.pi*5; # pole frequency (rad/s)
num = [1];      # transfer function numerator coefficients
den = c;        # transfer function denominator coefficients
lowPass = signal.TransferFunction(num,den) # Transfer function

# Generate the bode plot
w = np.logspace( np.log10(min(signalFreq)*2*np.pi/10), np.log10(max(signalFreq)*2*np.pi*10), 500 )
w, mag, phase = signal.bode(lowPass,w)

# Magnitude plot
plt.figure()
plt.semilogx(w, mag)
for sf in signalFreq:
    plt.semilogx([sf*2*np.pi,sf*2*np.pi],[min(mag),max(mag)],'k:')
plt.ylabel("Magnitude ($dB$)")
plt.xlim([min(w),max(w)])
plt.ylim([min(mag),max(mag)])

# Phase plot
plt.figure()
plt.semilogx(w, phase)  # Bode phase plot
plt.ylabel("Phase ($^\circ$)")
plt.xlabel("$\omega$ (rad/s)")
plt.xlim([min(w),max(w)])
plt.show()<jupyter_output><empty_output><jupyter_text>3. Discrete transfer functionTo implement the low-pass filter on hardware, you need to compute the discrete transfer function using the signal's sampling frequency.* The time step is $\Delta t = 1/f_s$* Compute the discrete transfer function using Tustin's method by setting $s = \frac{2}{\Delta t} \left( \frac{1-z^{-1}}{1+z^{-1}} \right)$* Why do it yourself? The to_discrete method computes the bilinear transform (Tustin's method when $\alpha = 1/2$)<jupyter_code># Compute the discrete low pass with delta_t = 1/samplingFrequency
dt = 1.0/samplingFreq;
discreteLowPass = lowPass.to_discrete(dt,method='gbt',alpha=0.5)
print(discreteLowPass)<jupyter_output>TransferFunctionDiscrete(
array([0.00024132, 0.00048264, 0.00024132]),
array([ 1.        , -1.95558189,  0.95654717]),
dt: 0.001
)<jupyter_text>4. Filter coefficientsWe want to find the filter coefficients for the discrete update:$$y[n] = a_1 y[n-1] + a_2 y[n-2] + ... + b_0 x[n] + b_1 x[n-1] + ...$$The coefficients can be taken directly from the discrete transfer function of the filter in the form:$$H(z) = \frac{b_0 + b_1 z^{-1} + b_2 z^{-2} + \ldots}{1 - a_1 z^{-1} - a_2 z^{-2} + \ldots}$$(This is a result of taking the Z-transform which is not shown here)Compare this to a transfer function with coefficientsnum = [b_0, b_1, b_2]den = [1, a_1, a_2]is $$H(z) = \frac{b_0 z^2 + b_1 z + b_2}{z^2 + a_1 z + a_2}$$which is equivalent to$$H(z) = \frac{b_0 + b_1 z^{-1} + b_2 z^{-2}}{1 + a_1 z^{-1} + a_2 z^{-2}}$$So you can take the coefficients in the same order that they are defined in the numerator and denominator of the transfer function object. The only difference is that the **coefficients in the denominator need a negative sign**.* To filter the signal, apply the filter using the discrete update* The filtered signal and filtered signal power spectrum are plotted alongside the unfiltered signal<jupyter_code># The coefficients from the discrete form of the filter transfer function (but with a negative sign)
b = discreteLowPass.num;
a = -discreteLowPass.den;
print("Filter coefficients b_i: " + str(b))
print("Filter coefficients a_i: " + str(a[1:]))

# Filter the signal
Nb = len(b)
yfilt = np.zeros(len(y));
for m in range(3,len(y)):
    yfilt[m] = b[0]*y[m];
    for i in range(1,Nb):
        yfilt[m] += a[i]*yfilt[m-i] + b[i]*y[m-i];
    
# View the result
# Plot the signal
plt.figure()
plt.plot(t,y);
plt.plot(t,yfilt);
plt.ylabel("$y(t)$")
plt.xlim([min(t),max(t)]);

# Generate Fourier transform
yfilthat = np.fft.fft(yfilt)
fcycles = np.fft.fftfreq(len(t),d=1.0/samplingFreq)

plt.figure()
plt.plot(fcycles,np.absolute(yhat));
plt.plot(fcycles,np.absolute(yfilthat));
plt.xlim([-100,100]);
plt.xlabel("$\omega$ (cycles/s)");
plt.ylabel("$|\hat{y}|$");<jupyter_output>Filter coefficients b_i: [0.00024132 0.00048264 0.00024132]
Filter coefficients a_i: [ 1.95558189 -0.95654717]
StarcoderdataJupyterStructuredCleanDedup