I am trying to implement a band-pass filter from scratch. Here is an article which I thought explained the nuts and bolts of how to build one; the author combines a low-pass with a high-pass filter (convolving both) to create a band-pass filter. I applied it to a mixture of sine waves, but the result was not satisfactory.
import numpy as np from numpy import sin, linspace, pi from pylab import figure, plot, show, title, xlabel, ylabel, subplot, savefig, xlim, ylim from scipy import fft, arange def plotSpectrum(y,Fs): n = len(y) # length of the signal k = arange(n) T = n/Fs frq = k/T # two sides frequency range frq = frq[range(n/2)] # one side frequency range Y = fft(y)/n # fft computing and normalization Y = Y[range(n/2)] plot(frq,abs(Y),'r') # plotting the spectrum xlabel('Freq (Hz)') ylabel('|Y(freq)|') Fs = 150.0; # sampling rate Ts = 1.0/Fs; # sampling interval t = arange(0,30,Ts) # time vector ff = 0.2; # frequency of the signal y = sin(2*pi*ff*t) ff = 0.8; # frequency of the signal y = y + sin(2*pi*ff*t) figure() subplot(2,1,1) plot(t,y) xlabel('Time') ylabel('Amplitude') subplot(2,1,2) plotSpectrum(y,Fs) savefig('out1.png') fL = 0.1 fH = 0.3 N = 33 n = np.arange(N) hlpf = np.sinc(2 * fH * (n - (N - 1) / 2.)) hlpf *= np.blackman(N) hlpf = hlpf / np.sum(hlpf) # Compute a high-pass filter with cutoff frequency fL. hhpf = np.sinc(2 * fL * (n - (N - 1) / 2.)) hhpf *= np.blackman(N) hhpf = hhpf / np.sum(hhpf) hhpf = -hhpf hhpf[(N - 1) / 2] += 1 hp = np.convolve(hlpf, hhpf) filtered = np.convolve(y, hp) t2 = np.linspace(0,30.,len(filtered)) # time vector figure() subplot(2,1,1) plot(t2, filtered) xlabel('Time') ylabel('Amplitude') subplot(2,1,2) plotSpectrum(filtered,Fs) savefig('out2.png')
For the original signal:
I get an output like this:
To confirm I could filter this signal, I used a butterworth filter,
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 filtered = butter_bandpass_filter(y, fL, fH, Fs, order=2)
That worked fine (although I had to play around with the order a little bit, higher orders of this filter does not work either). What is wrong with the first filter?