0
$\begingroup$

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))

enter image description here

enter image description here

$\endgroup$
0
$\begingroup$

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.

enter image description here

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)
| improve this answer | |
$\endgroup$
  • $\begingroup$ I don't get a correct plot for the signal... $\endgroup$ – Robinson-Constantin Chera Apr 17 at 17:58
  • $\begingroup$ Please update your code in the questions so that I can see whats wrong. $\endgroup$ – jithin Apr 17 at 17:59
  • $\begingroup$ Updated the code now $\endgroup$ – Robinson-Constantin Chera Apr 17 at 18:00
  • $\begingroup$ The 'x' in my code is different from your 'x'. My 'x' variable is an 1-D array of '1' followed by 64 zeros. So rename your original signal 'x' as 'x1'. Then, you need to multiply 'w' and 'x' point by point (both should be 65 samples long if everything is correct). Let the output be x_hann. The finally you have to do signal.lfilter(x_hann,1,x1) . I will try to edit your code. $\endgroup$ – jithin Apr 17 at 18:05
  • $\begingroup$ first I need to implement using FIR Hann window and then implement an IIR filter using Cheby type 1 in order to compare the results.. $\endgroup$ – Robinson-Constantin Chera Apr 17 at 18:06
0
$\begingroup$

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')

| improve this answer | |
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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