0
$\begingroup$

I'm trying to implement 2 lowpass filters and see the difference between them from plots( one FIR filter that needs to be designed with the Hann window method and another IIR filter using Chebyshev type I). The problem is that I can't make 2 correct plots to see the difference between them and I don't know exactly if they are correctly implemented...I have fc1,fc2, ft, rp and rs as parameters. Here it's my python code:

 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 = readWav('Ring05c.wav')

plt.close('all')
plt.figure(figsize=[12,12])
plt.subplot(3,2,1)
plotInTime(x,fs)
plt.title('Original Signal')
plt.subplot(3,2,2)
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)
L, Wn = signal.cheb1ord(fc1/(fs/2), fc2/(fs/2), rp, rs);
b, a = signal.cheby1(L, 1, Wn,'low')
h = signal.lfilter(b, a,x,axis=0)
plotInFrequnecy(h,fs)
w=signal.hann(65) #idk what to do after this
$\endgroup$
8
  • $\begingroup$ to plot the frequency response of your filter once you have the b, a coefficients you can use h, w = signal.freqz(b,a) and then you can plot the magnitude and phase of h versus frequency w easily. $\endgroup$ – Dan Boschen Apr 17 '20 at 21:33
  • $\begingroup$ As for designing FIR filters using the window method, what have you already learned about that? $\endgroup$ – Dan Boschen Apr 17 '20 at 21:36
  • 2
    $\begingroup$ You created a window, but not the filter---what's the procedure to design a filter using the windowing method? Put some details before the code outlining that to show what you know about it and then I can help you where you are getting stuck.(not in the code but just to show what you have studied and then we can better help you where you are stuck) $\endgroup$ – Dan Boschen Apr 18 '20 at 4:07
  • 1
    $\begingroup$ @DanBoschen I ll come back in a few hours with the Details and updated code $\endgroup$ – Robinson Chera Apr 18 '20 at 4:14
  • 1
    $\begingroup$ Ok the hint is that to make a brickwall filter in frequency (rectangular function in frequency) requires an infinitely long function in time given the impulse response in time is the inverse Fourier Transform of the frequency response. You don't have infinite time so you have to window in time-- if you just truncate it in time it would create sharp artifacts in freq so you use a gradual window instead. So you make the ideal impulse response, and then you window that $\endgroup$ – Dan Boschen Apr 18 '20 at 4:17
0
$\begingroup$

Right now I have an error that says: Error: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()

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


        fs, x = myDSP.readWav('Ring05c.wav')



        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=np.array([fc-ft/2])/(fs/2)
ws=np.array([fc+ft/2])/(fs/2)
        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)
        plt.plot(b,'.-b',label='b coefficients')
        plt.plot(a,'.-r',label='a coefficients')
        plt.legend()


        #h = signal.lfilter(b, a,x,axis=0)
        #w=signal.hann(65)
$\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.