Time-domain LPF not showing expected behavior

Question

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!

Answer 1

Your math is correct but your code doesn't match it:

y = (wc*Ts*(x-xprev)+(2-wc*Ts)*yprev)/(2+wc*Ts)

It should be $x+xprev$ not $x-xprev$. The minus sign makes it a highpass and hence the attenuation.