Implementing a Band-Pass Filter from Scratch

Question

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?

Regards,

Answer 1

A simple way to make a real bandpass filter (assuming you know how to make a low pass filter) is multiply the coefficients of your prototype low pass filter by $cos(2\pi f_c nT)$ where $f_c$ is the desired center frequency of your passband and T is the sampling interval for your coefficients.

Answer 2

The following Matlab codes were based on Steven W. Smith's DSP book,

http://www.gomatlab.de/window-sinc-filter-t19156.html

which I converted into Python;

import numpy as np

def hamm(window_size):
    N = window_size;
    output = np.zeros((N, 1));
    if np.mod(N, 2) == 0 :
        m = np.fix(N / 2)
        n = m
    else:
        m = np.fix(N / 2)+1; 
        n = m-1; 
    window = 0.54 - 0.46 * np.cos(2*np.pi*(np.arange(m)) / (N-1))
    tmp1 = window[:int(m)]
    tmp2 = window[np.arange(int(n)-1,-1,-1)]
    return np.hstack((tmp1,tmp2))

def sinc_filter_low(order, fc1, fs):
    Fc1 = fc1 / np.float(fs) 
    M  = order
    B = np.zeros((M+1, 1))
    window = hamm(M+1)
    for i in range(M+1):
        if 2 * i == M:
            B[i] = 2*np.pi*Fc1
        else:
            tmp1 = 2*np.pi*Fc1 *(i-(M/2.))
            tmp2 = (i-(M/2.))
            B[i] = np.sin(tmp1) / tmp2
        B[i] = B[i] * window[i]
    return B / np.sum(B)

def sinc_filter_high(order, fc1, fs):
    Fc1 = fc1 / np.float(fs) 
    M  = order
    B = np.zeros((M+1, 1))
    window = hamm(M+1)
    for i in range(M+1):
        if 2 * i == M:
            B[i] = 2*np.pi*Fc1
        else:
            tmp1 = 2*np.pi*Fc1 *(i-(M/2.))
            tmp2 = (i-(M/2.))
            B[i] = np.sin(tmp1) / tmp2
        B[i] = B[i] * window[i]
    B = B / np.sum(B)
    B = -B
    B[(M/2)] = B[(M/2)] + 1
    return B

def sinc_filter_band(order, fc1, fc2, fs):
    M = order
    A = sinc_filter_low(order, fc1, fs).T[0]
    B = sinc_filter_high(order, fc2, fs).T[0]
    output = A+B
    output = -output
    output[(M/2)] = output[(M/2)] + 1.
    return output

These work fine.