# Time-domain LPF not showing expected behavior

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!

• Hi! It would be better if you could separate the coding from mathematics. The syntax specific issues simply blur the vision... Sep 24, 2018 at 20:25
• That's good advice. I'll try rewrite in a way that's closer to the pure math and edit my original post. Sep 24, 2018 at 20:37

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.