I am trying to implement a simple first-order Butterworth Low-Pass filter in Python. I have some code that makes use of scipy.signal.butter and scipy.signal.filtfilt. It works fine, but I wanted to make sure I could build the filter myself as well.
I start with a bilinear transform $s = \frac{2}{T_{s}}\frac{z-1}{z+1}$ on my first order LPF $\frac{\omega_{c}}{\omega_{c}+s}$. Making that substitution and converting from the $z$ domain to the time domain gives me the equation that I'm trying to reproduce:
$y_{n} = \frac{1}{\omega_{c}T_{s}+2}[\omega_{c}T_{s}(x_{n}+x_{n-1})+(2-\omega_{c}T_{s})y_{n-1}]$
I'm confident that this expression is correct, but my attempt at coding it doesn't seem to be working.
EDIT: I changed my implementation so what's actually happening is more clear. As a side effect, I'm now seeing a different type of strange behavior than I was previously. The new code is below:
T = 2 # time length over which samples are collected (s)
f0 = 5 # fundamental frequency (Hz)
fs = 1000 # sampling frequency
Ts = 1/fs # sampling period
fc = 10 # filter cutoff frequency
wc = 2*pi*fc # cutoff frequency in radians
# Create wave input to be sent to filter
t = np.linspace(0,T,T*fs+1)
x_storage = [10 + cos(2*pi*f0*x) for x in t]
plt.plot(t,x_storage)
# Create a place to store output data
y_storage = np.zeros(len(t))
# Initialize input and output data
y = 0 # y_(n)
yprev = 0 # y_(n-1)
x = 0 # x_(n)
xprev = 0 # x_(n-1)
# Begin calculation
for idx in range(len(t)):
# One timestep of the calculation
xprev = x # move the old x_(n) value to x_(n-1)
x = x_storage[idx] # bring in the current datapoint to be the new x_(n)
yprev = y # move the old y_(n) value to y_(n-1)
y = (wc*Ts*(x-xprev)+(2-wc*Ts)*yprev)/(2+wc*Ts) # calculate the new y_(n) according to the filter equation
y_storage[idx] = y # store the y_(n) value for later plotting
plt.figure()
plt.plot(t,y_storage)
My input is a wave with a DC offset and a frequency component below the cutoff threshold, so I would expect the output wave to look more-or-less the same. Instead, I see the DC bias vanish and the signal itself get hugely attenuated (see below). Any thoughts would be much appreciated!