Lowpass Hann filter in Python

Question

I'm trying to design a FIR lowpass filter using Hann windowing in Python but I don't know where to start from... I have fc1,fc2, ft, rp and rs as parameters. This is my Python code:

            import myDSP
        import matplotlib.pyplot as plt
        import scipy.signal as signal
        import numpy as np

        def plotInTime(x,fs):
t=np.arange(0,np.size(x,0)/fs,1/fs)
plt.figure()
plt.plot(t,x)
plt.xlabel('Time (s)')
plt.ylabel('Amplitude')
plt.grid(True)

def plotInFrequency(x,fs):
N=int(np.size(x,0)/2)
if np.size(x,0)==1:
X=np.fft.fft(x,axis=1)
else:
X=np.fft.fft(x,axis=0)
X=np.abs(X)
X=X[:N]
f=np.arange(0,fs/2,fs/2/N)
plt.figure()
plt.plot(f,X)
plt.xlabel('Frequency (Hz)')
plt.ylabel('Magnitude')
plt.grid(True)
        fs, x = myDSP.readWav('Ring05c.wav') #is 22050
        myDSP.plotInTime(x,fs)
        plt.title('Original Signal')
        myDSP.plotInFrequency(x,fs)
        plt.title('Spectrum for the Original Signal')

        fc1=1200
        fc2=3000
        fc= np.arange(fc1,fc2+1)
        ft=150
        rp=0.1
        rs=40
        t = np.linspace(-1, 1, 201)
        #x = (np.sin(2*np.pi*0.75*t*(1-t) + 2.1)+0.1*np.sin(2*np.pi*1.25*t + 1) +0.18*np.cos(2*np.pi*3.85*t))
        wp=(fc-ft/2)/(fs/2)
        wp = wp[0]
        ws=(fc+ft/2)/(fs/2)
        ws = ws[0]
        L, Wn = signal.cheb1ord(wp,ws, rp, rs);
        b, a = signal.cheby1(L, rp, Wn,'low')
        f,H=signal.freqz(b,a,worN=256,plot=None,fs=fs)    

        y=signal.lfilter(b,a,x,axis=0)
        myDSP.plotInTime(y,fs)
        plt.title('Filtered Cheby1 Signal')

        myDSP.plotInFrequency(y,fs)
        plt.title('Spectrum for the Filtered Cheby1 Signal')
        y=y.astype('int16')

        w_fir = signal.firwin(65, fc/(fs/2), window='hann')
    [w,h] = signal.freqz(w_fir)
    plt.plot(w, 20*np.log10(abs(h)))
    plt.plot(w, np.angle(h))

Answer 1

My recommendation would be to generate transfer function for Chebyshev Type I IIR filter and then convert to FIR by Hann windowing. The steps involved in converting IIR to FIR might be a little difficult but this post will help - Manually Windowing an IIR Filter in Matlab

Generating Chebyshev Type 1 IIR Filter: First determine order of filter using cheb1ord(https://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.signal.cheb1ord.html#scipy.signal.cheb1ord) and then generate the transfer function of filter using scipy.signal.cheby1 (https://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.signal.cheby1.html). To achieve 0.1 ripple below unity gain, you need to specify attenuation as 1 (so that $10^{-1/20} \approx 0.9$)

L, Wn = signal.cheb1ord(fc1/(fs/2), fc2/(fs/2), 1, 40);
b, a = signal.cheby1(L, 1, Wn,'low')

Now, from $a$ and $b$, you can get the IIR Impulse response $h$ by convolving with unit impulse response $\delta[n]$ for $N$ samples to get $N$ length filter. You can use scipy.signal.lfilter python command (https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.lfilter.html) for this

x = np.zeros(65)
x[0] = 1    
h = signal.lfilter(b, a, x)
myDSP.plotInFrequency(h,fs)

where $x$ will be unit impulse impulse response $\delta[n]$ of length $N$ (1 followed by N-1 zeros).

Then generate Hann window of size $N$, refer (https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.windows.hann.html?highlight=hann#scipy.signal.windows.hann)

w=signal.hann(N)

$w[n]$, and do point by point multiplication over $h$. That is $h[n]w[n]$ for $0 \le n\le N-1$.

UPDATE: After doing the same above in MATLAB, I have found the optimum $N=65$ so that the Hann windowed FIR filter meets the pass band ripple and stop band attenuation. Bet ween $\omega=0$ and $\omega \approx 0.05\pi$(fc1), the passband ripple within 1dB. At $\omega=0.136\pi$(fc2), the stop band attenuation is more than 40dB meeting the filter specs with FIR Hann Window.

UPDATE: Code to add FIR LPF with Hann window using (https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.firwin.html)

w_fir = signal.firwin(65, fc1/(fs/2), window='hann')
myDSP.plotInFrequency(w_fir,fs)

Answer 2

You should read the documentation when using new functions: scipy.signal.firwin

1. You either provide pass_zero option or two cutoff frequencies. So in your case you don't have to use this option.
2. The first argument corresponds to the length of the filter. So it has to be an odd integer.
3. You either provide the sampling frequency or express the cutoff frequencies as relative frequencies (in your case something like fc1/fs).

In summary, for a a filter of order n, you can get your coefficients as

signal.firwin(n + 1, [fc1, fc2], fs=fs, window='hann')