1

I'm trying to use a Butterworth filter in Python as described in this thread with these functions:

def butter_bandpass(lowcut, highcut, fs, order=5):
    nyq = 0.5 * fs
    low = lowcut / nyq
    high = highcut / nyq
    b, a = butter(order, [low, high], btype='band')
    return b, a


def butter_bandpass_filter(data, lowcut, highcut, fs, order=5):
    b, a = butter_bandpass(lowcut, highcut, fs, order=order)
    y = lfilter(b, a, data)
    return y

The output of the FFT of my data without applying the filter gives the following plot:

enter image description here

However, after applying the filter above with:

lowcut = 1.0
highcut = 50.0
x2_Vtcr = butter_bandpass_filter(x_Vtcr, lowcut, highcut, fs, order=4)

where fs is the sampling frequency (1000 in my case) I get as FFT:

enter image description here

It looks like the filter shifts the frequency to the left and I don't get the peak where it should be. Any idea as to why this is happening?

Just in case I place where the data file (second column). The DC is quite strong so after FFT the first element should not be included in the plot.

Thank you!

Community
  • 1
  • 1
Janus Gowda
  • 295
  • 1
  • 4
  • 17

1 Answers1

1

The 2 functions above seem to work fine here. This is an example of a white noise signal sampled at your fs=250, and filtered with the two functions you mention, with a passband between 1Hz and 50 Hz as you do:

from numpy import random, arange
from numpy.fft import rfft, rfftfreq

fs = 250.
t = arange(0., 30., 1 / fs)
x_Vtcr = random.randn(len(t))

hat_x_Vtcr = rfft(x_Vtcr)
freqs = rfftfreq(x_Vtcr.size, 1 / fs)
plt.plot(freqs, abs(hat_x_Vtcr), 'b')

lowcut, highcut = 1.0, 50.0
x2_Vtcr = butter_bandpass_filter(x_Vtcr, lowcut, highcut, fs, order=4)

hat_x2_Vtcr = rfft(x2_Vtcr)
plt.plot(freqs, abs(hat_x2_Vtcr), 'r')

The resulting PSD is fine with me (red is the filtered one): enter image description here

I guess your error is somewhere else? Please compare your code with the above snippet. You may also want to read THIS for next time.

EDIT:reply to the comment.

It does indeed work also with your data. I did:

datafile=loadtxt('V.dat')
x_Vtcr=datafile[:,1]
x_Vtcr-=x_Vtcr.mean()

and then I ran the above script, without the line on which I generate the x_Vtcr data of course. The resulting filtered output follows: enter image description here

Community
  • 1
  • 1
gg349
  • 21,996
  • 5
  • 54
  • 64
  • for me it also works with dummy examples. But it does not happen with the data I have. I have posted the data in the question itself. – Janus Gowda Nov 05 '14 at 15:25
  • My problem is not to get an output. I do get also an output and I posted it in the image above. My problem is to get an output with the peak at the same frequencies. – Janus Gowda Nov 05 '14 at 17:08
  • @JanusGowda: I did not post the code to give you AN output, but to show you how to apply the bandpass filter. In the figure, the red and the blue signals in the 2nd figure above have the peaks exactly at the same frequencies (inside the band of the filter obviously). Please run my script and check for yourself. – gg349 Nov 05 '14 at 18:19