0
$\begingroup$

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!

enter image description here

enter image description here

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

1 Answer 1

0
$\begingroup$

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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