FIR Filter Design in Matlab code. Please use the following code to run it. Specification: 1. Filt Filter er Type Type:: Low Low Pass Pass 2. PassB PassBand and Frequ Frequenc ency(F y(Fp) p) = 1KHz 1KHz 3. StopBa StopBand nd Freq Frequen uency( cy(Fs) Fs) = 4.3 4.3 KHz KHz 4. Sampl Sampling ing Frequ Frequenc ency y (Fsam (Fsam)) = 10KHz 10KHz Transition Width(Δw) Normalized Normalized = (4300-1000)/10,000 = 0.33 Cut off frequency(Normalized) frequency(Normalized) = 1000/(10,000/2)=0.2 1000/(10,000/2)=0.2 Approximation: Order of Filter = 3.3/Δw = 10; N = 10;a1=1;a2=1;Fs=10000;Wc= 10;a1=1;a2=1;Fs=10000;Wc=0.2; 0.2; b1 = fir1(N-1,Wc,'low' fir1(N-1,Wc,'low' ,hamming(N)); ,hamming(N)); [db1, w1] = freqz(b1, a1, 512, Fs); mag1 = 20*log10(abs(db1)); subplot(2,2,1);plot(w1, subplot(2,2,1);plot(w1, mag1), grid on; xlabel('Frequency xlabel('Frequency (Hz)');ylabel('Magnitude (Hz)');ylabel('Magnitude Response (dB)'); title('FIR Low Pass using Hamming Window'); subplot(2,2,2);stem(b1);grid subplot(2,2,2);stem(b1);grid on; xlabel('n+1');ylabel('b(n)');title('F xlabel('n+1');ylabel('b(n)');title('Filter ilter Coefficient'); Coefficient'); b2 = fir1(N-1,Wc,'low' fir1(N-1,Wc,'low' ,boxcar(N)); ,boxcar(N)); [db2, w2] = freqz(b2, a2, 512, Fs); mag2 = 20*log10(abs(db2)); subplot(2,2,3);plot(w2, subplot(2,2,3);plot(w2, mag2), grid on; xlabel('Frequency xlabel('Frequency (Hz)');ylabel('Magnitude (Hz)');ylabel('Magnitude Response (dB)'); title('FIR Low Pass using Rectangular window'); subplot(2,2,4);stem(b2);grid subplot(2,2,4);stem(b2);grid on; xlabel('n+1');ylabel('b(n)');title('F xlabel('n+1');ylabel('b(n)');title('Filter ilter Coefficient'); Coefficient');
Specification: 5. Filt Filter er Type Type:: High High Pass Pass 6. Other Other filter filter specifi specificati cation on same same as that that of Low Low Pass Pass filter. filter. N = 10;a1=1;a2=1;Fs=10000;Wc= 10;a1=1;a2=1;Fs=10000;Wc=0.2; 0.2; b1 = fir1(N-1,Wc,'high' fir1(N-1,Wc,'high' ,hamming(N+1)); ,hamming(N+1)); [db1, w1] = freqz(b1, a1, 512, Fs); mag1 = 20*log10(abs(db1)); subplot(2,2,1);plot(w1, subplot(2,2,1);plot(w1, mag1), grid on; xlabel('Frequency xlabel('Frequency (Hz)');ylabel('Magnitude (Hz)');ylabel('Magnitude Response (dB)'); title('FIR High Pass using Hamming Window'); subplot(2,2,2);stem(b1);grid subplot(2,2,2);stem(b1);grid on; xlabel('n+1');ylabel('b(n)');title('F xlabel('n+1');ylabel('b(n)');title('Filter ilter Coefficient'); Coefficient'); b2 = fir1(N-1,Wc,'high' fir1(N-1,Wc,'high' ,boxcar(N+1)); ,boxcar(N+1)); [db2, w2] = freqz(b2, a2, 512, Fs); mag2 = 20*log10(abs(db2)); subplot(2,2,3);plot(w2, subplot(2,2,3);plot(w2, mag2), grid on; xlabel('Frequency xlabel('Frequency (Hz)');ylabel('Magnitude (Hz)');ylabel('Magnitude Response (dB)');
title('FIR High Pass using Rectangular window'); subplot(2,2,4);stem(b2);grid on; xlabel('n+1');ylabel('b(n)');title('Filter Coefficient'); Specification: 7. Filter Type: Band Pass 8. PassBand Frequency(Fp1&Fp2) = 9KHz,17KHz. 9. StopBand Frequency(Fs1,Fs2) = 8 KHz,16KHz 10. Sampling Frequency (Fsam) = 48 KHz Transition Width(Δw) Normalized = (9000-8000)/48,000 = 0.020833 Cut off frequency(Normalized) : Wn1 = 9000/24000 = 0.375 ; Wn2 = 17000/24000=0.7083; Wn = [0.375,0.7083] Approximation: Order of Filter = 3.3/Δw = 159; N = 159;Fs=48000;a1=1; Wn = [0.375,0.7083]; b1 = fir1(N-1,Wn,'band' ,hamming(N)); [db1, w1] = freqz(b1, a1, 512, Fs); mag1 = 20*log10(abs(db1)); subplot(2,2,1);plot(w1, mag1), grid on; xlabel('Frequency (Hz)');ylabel('Magnitude Response (dB)'); title('FIR Band Pass using Hamming Window'); subplot(2,2,2);stem(b1);grid on; xlabel('n+1');ylabel('b(n)');title('Filter Coefficient'); a2=1; b2 = fir1(N-1,Wn,'band' ,boxcar(N)); [db2, w2] = freqz(b2, a2, 512, Fs); mag2 = 20*log10(abs(db2)); subplot(2,2,3);plot(w2, mag2), grid on; xlabel('Frequency (Hz)');ylabel('Magnitude Response (dB)'); title('FIR Band Pass using Rectangular window'); subplot(2,2,4);stem(b2);grid on; xlabel('n+1');ylabel('b(n)');title('Filter Coefficient');
Specification: 11. Filter Type: Band Stop 12. Other Filter Specification same as that of Band Pass as given above. N = 159;Fs=48000;a1=1; a2=1; Wn = [0.375,0.7083]; b1 = fir1(N-1,Wn,'stop' , kaiser(N)); [db1, w1] = freqz(b1, a1, 512, Fs); mag1 = 20*log10(abs(db1)); subplot(2,2,1);plot(w1, mag1), grid on; xlabel('Frequency (Hz)');ylabel('Magnitude Response (dB)'); title('FIR Band Stop using Kaiser Window'); subplot(2,2,2);stem(b1);grid on; xlabel('n+1');ylabel('b(n)');title('Filter Coefficient');
b2 = fir1(N-1,Wn,'stop' ,hanning(N)); [db2, w2] = freqz(b2, a2, 512, Fs); mag2 = 20*log10(abs(db2)); subplot(2,2,3);plot(w2, mag2), grid on; xlabel('Frequency (Hz)');ylabel('Magnitude Response (dB)'); title('FIR Band Stop using Hanning window'); subplot(2,2,4);stem(b2);grid on; xlabel('n+1');ylabel('b(n)');title('Filter Coefficient');
IIR Filter Design: Utilities Scripts:: LP_butterworth.script function [b,a] = LP_butterworth(Wp,Ws,Rp,As); %N=Order of filter %b,a= Numerator & Denominator Coefficient of the filter. %Wc,Wp,Ws,Rp,As= Cuttoff frequency,Passband frequency %Stopband frequency,Passband Ripple & Stopband attenuation. N = ceil((log10((10^(Rp/10)-1)/(10^(As/10)-1)))/(2*log10(Wp/Ws))); fprintf('\n Butterworth Filter Order = %2.0f \n',N) Wc = Wp/((10^(Rp/10)-1)^(1/(2*N))); [b,a]= unnorm_buttap(N,Wc); unnorm_buttap.script function [b,a] = unnorm_buttap(N, Wc); % Unnormalized Butterworth Analog Lowpass Filter [z,p,k] = buttap(N); p = p* Wc; k = k* Wc ^N; B = real(poly(z)); b0 = k; b = k*B; a = real(poly(p)); LP_elliptic.script function [b,a] = LP_elliptic (Wp,Ws,Rp,As); %N=Order of filter %b,a= Numerator & Denominator Coefficient of the filter. %Wc,Wp,Ws,Rp,As= Cuttoff frequency,Passband frequency %Stopband frequency,Passband Ripple & Stopband attenuation. ep = sqrt(10^(Rp/10)-1); A = 10^(As/20); Wc = Wp; k = Wp/Ws; k1 = ep/sqrt(A*A-1); capk = ellipke([k.^2 1-k.^2]); % Version 4.0 code
capk1 = ellipke([(k1 .^2) 1-(k1 .^2)]); % Version 4.0 code N = ceil(capk(1)*capk1(2)/(capk(2)*capk1(1))); fprintf('Elliptic Filter Order = %2.0f \n',N) [b,a]= unnorm_elliptic (N,Rp,As, Wc); unnorm_elliptic.script function [b,a] = unnorm_elliptic(N,Rp,As,Wc); % Unnormalized Elliptic Analog Lowpass Filter. [z,p,k] = ellipap(N,Rp,As); a = real(poly(p)); aNn = a(N+1); p = p* Wc; a = real(poly(p)); aNu = a(N+1); b = real(poly(z)); M = length(b); bNn = b(M); z = z* Wc; b = real(poly(z)); bNu = b(M); k = k*(aNu*bNn)/(aNn*bNu); b0 = k; b = k*b; LP_chbeshyv1.script function [b,a] = LP_chbeshyv1(Wp,Ws,Rp,As); %N=Order of filter %b,a= Numerator & Denominator Coefficient of the filter. %Wc,Wp,Ws,Rp,As= Cuttoff frequency,Passband frequency %Stopband frequency,Passband Ripple & Stopband attenuation. ep = sqrt(10^(Rp/10)-1); A = 10^(As/20); Wc = Wp; Wr = Ws/Wp; g = sqrt(A*A-1)/ep; N = ceil(log10(g+sqrt(g*g-1))/log10(Wr +sqrt(Wr * Wr -1))); fprintf('Chebyshev-1 Filter Order = %2.0f \n',N) [b,a]= unnorm_chbeshyv1(N,Rp, Wc); unnorm_chbeshyv1.script function [b,a] = unnorm_chbeshyv1 (N,Rp,Wc); % Unnormalized Chebyshev-1 Analog Lowpass Filter. [z,p,k] = cheb1ap(N,Rp); a = real(poly(p)); aNn = a(N+1); p = p* Wc; a = real(poly(p)); aNu = a(N+1); k = k*aNu/aNn; b0 = k; B = real(poly(z)); b = k*B;
LP_chbeshyv2.script function [b,a] = LP_chbeshyv2 (Wp,Ws,Rp,As); %N=Order of filter %b,a= Numerator & Denominator Coefficient of the filter. %Wc,Wp,Ws,Rp,As= Cuttoff frequency,Passband frequency %Stopband frequency,Passband Ripple & Stopband attenuation. ep = sqrt(10^(Rp/10)-1); A = 10^(As/20); OmegaC = Wp; OmegaR = Ws/Wp; g = sqrt(A*A-1)/ep; N = ceil(log10(g+sqrt(g*g-1))/log10(OmegaR+sqrt(OmegaR*OmegaR-1))); fprintf('Chebyshev-2 Filter Order = %2.0f \n',N) [b,a]= unnorm_chbeshyv2 (N,As,Ws); unnorm_chbeshyv2.script function [b,a] = unnorm_chbeshyv2 (N,As,Wc); % Unnormalized Chebyshev-2 Analog Lowpass. [z,p,k] = cheb2ap(N,As); a = real(poly(p)); aNn = a(N+1); p = p* Wc; a = real(poly(p)); aNu = a(N+1); b = real(poly(z)); M = length(b); bNn = b(M); z = z* Wc; b = real(poly(z)); bNu = b(M); k = k*(aNu*bNn)/(aNn*bNu); b0 = k; b = k*b; freqs_m.script function [db,mag,pha,w] = freqs_m(b,a,wmax); % Computation of s-domain frequency response: % -----------------------------------------------------------% [db,mag,pha,w] = freqs_m(b,a,wmax); % db = Relative magnitude in db over [0 to wmax] % mag = Absolute magnitude over [0 to wmax] % pha = Phase response in radians over [0 to wmax] % w = array of 500 frequency samples between [0 to wmax] % b = Numerator polynomial coefficents of Ha(s) % a = Denominator polynomial coefficents of Ha(s) % wmax = Maximum frequency in rad/sec over which response is desired % w = [0:1:500]*wmax/500; H = freqs(b,a,w); mag = abs(H); db = 20*log10((mag+eps)/max(mag)); pha = angle(H); imp_invr.script
function [b,a] = imp_invr(c,d,T) % Impulse Invariance Transformation from Analog to Digital Filter % --------------------------------------------------------------% [b,a] = imp_invr(c,d,T) % b = Numerator polynomial in z^(-1) of the digital filter % a = Denominator polynomial in z^(-1) of the digital filter % c = Numerator polynomial in s of the analog filter % d = Denominator polynomial in s of the analog filter % T = Sampling (transformation) parameter % [R,p,k] = residue(c,d); p = exp(p*T); [b,a] = residuez(R,p,k); b = real(b'); a = real(a'); freqz_m.script function [db,mag,pha,grd,w] = freqz_m(b,a); % Modified version of freqz subroutine % -----------------------------------% [db,mag,pha,grd,w] = freqz_m(b,a); % db = Relative magnitude in dB computed over 0 to pi radians % mag = absolute magnitude computed over 0 to pi radians % pha = Phase response in radians over 0 to pi radians % grd = Group delay over 0 to pi radians % w = 501 frequency samples between 0 to pi radians % b = numerator polynomial of H(z) (for FIR: b=h) % a = denominator polynomial of H(z) (for FIR: a=[1]) % [H,w] = freqz(b,a,1000,'whole'); H = (H(1:1:501))'; w = (w(1:1:501))'; mag = abs(H); db = 20*log10((mag+eps)/max(mag)); pha = angle(H); grd = grpdelay(b,a,w);
Impseq.script function [x,n] = impseq(n0,n1,n2) % Generates x(n) = delta(n-n0); n1 <= n,n0 <= n2 % ---------------------------------------------% [x,n] = impseq(n0,n1,n2) % if ((n0 < n1) | (n0 > n2) | (n1 > n2)) error('arguments must satisfy n1 <= n0 <= n2') end n = [n1:n2]; %x = [zeros(1,(n0-n1)), 1, zeros(1,(n2-n0))]; x = [(n-n0) == 0]; Specification: 13. Filter Type: Low Pass(Butterworth) 14. PassBand Frequency(Wp) = 30 rad/sec. 15. StopBand Frequency(Ws) = 40 rad/sec.
16. Passband Ripple= 1. 17. Stopband Attunation=30.
Wp=30;Ws=40;Rp=1;As=30; [b,a] = LP_butterworth(Wp,Ws,Rp,As); Wmax=100; [db,mag,pha,w]=freqs_m(b,a,Wmax); pha = unwrap(pha); [b,a]= LP_butterworth(Wp/10,Ws/10,Rp,As); [ha,x,t]=impulse(b,a); t=t/10;ha=ha/10; %Plots subplot(2,2,1);plot(w,mag);axis([0,Wmax,0,1.1]);grid on; xlabel('Analog frequency in rad/sec(w)');ylabel('Absolute Magnitude(mag)'); title('Analog Butterworth LowPass filter.'); subplot(2,2,2);plot(w,db);axis([0,Wmax,-100,0]); grid on; xlabel('Analog frequency in rad/sec'); ylabel('Magnitude in db'); subplot(2,2,3);plot(w,pha/pi);axis([0,Wmax,floor(min(pha/pi)),floor(ceil(max(pha/pi)))]); grid on; xlabel('Analog frequency in rad/sec'); ylabel('Phase in pi units'); subplot(2,2,4);plot(t,ha);title('Impulse Response'); grid on; xlabel('t(sec) ');ylabel('ha(t) '); Specification: 18. Filter Type: Low Pass(Elliptic) 19. PassBand Frequency(Wp) = 30 rad/sec. 20. StopBand Frequency(Ws) = 40 rad/sec. 21. Passband Ripple= 1. 22. Stopband Attunation=30. Wp=30;Ws=40;Rp=1;As=30; [b,a] = LP_elliptic(Wp,Ws,Rp,As); a0=a(1);b=b/a0;a=a/a0; Wmax=100; [db,mag,pha,w]=freqs_m(b,a,Wmax);pha=unwrap(pha); [ha,x,t]=impulse(b,a); title('Analog Low Pass Elliptic Filter'); subplot(2,2,1);plot(w,mag);axis([0,Wmax,0,1]);grid on; title('Analog Low Pass Elliptic Filter'); xlabel('Analog frequency in rad/sec(w)');ylabel('Absolute Magnitude(mag)'); subplot(2,2,2);plot(w,db);axis([0,Wmax,-100,0]);grid on; xlabel('Analog frequency in rad/sec(w)');ylabel('Log-Magnitude in dB'); minpha=floor(min(pha/pi));maxpha=ceil(max(pha/pi)); subplot(2,2,3);plot(w,pha/pi);axis([0,Wmax,minpha,maxpha]);grid on; xlabel('Analog frequency in rad/sec(w)'):ylabel('Phase in pi units'); subplot(2,2,4);plot(t,ha);title('Impulse Response');grid on; xlabel('t (sec) ');ylabel('ha(t) '); Specification: 23. Filter Type: Low Pass(Chebyshev-I) 24. PassBand Frequency(Wp) = 30 rad/sec. 25. StopBand Frequency(Ws) = 40 rad/sec. 26. Passband Ripple= 1. 27. Stopband Attunation=30.
Wp=30;Ws=40;Rp=1;As=30; [b,a]= LP_chbeshyv1(Wp,Ws,Rp,As); Wmax=100; [db,mag,pha,w]=freqs_m(b,a,Wmax);pha=unwrap(pha); subplot(2,1,1);plot(w/(2*pi),mag);grid on; title('Analog Low Pass(Chebyshev-I)'); xlabel('Analog frequency in Hz(w)');ylabel('Magnitude(mag)'); subplot(2,1,2);plot(w/(2*pi),db);grid on; xlabel('Analog frequency in Hz');ylabel('Decibels(db)'); title('Log-Magnitude Response'); Specification: 28. Filter Type: Low Pass(Chebyshev-II) 29. PassBand Frequency(Wp) = 30 rad/sec. 30. StopBand Frequency(Ws) = 40 rad/sec. 31. Passband Ripple= 1. 32. Stopband Attunation=30. Wp=30;Ws=40;Rp=1;As=30; [b,a]= LP_chbeshyv2(Wp,Ws,Rp,As); Wmax=100; [db,mag,pha,w]=freqs_m(b,a,Wmax);pha=unwrap(pha); subplot(2,2,1);plot(w/(2*pi),mag);grid on; title('Analog Low Pass Frequency(Chebyshev-II)'); xlabel('Analog frequency in Hz(w)');ylabel('Absolute Magnitude(mag)'); title('Absolute magnitude Response'); subplot(2,2,2);plot(w/(2*pi),db);grid on; xlabel('Analog frequency in Hz(w)');ylabel('Relative Magnitude in Decibels(db)'); title('Log-Magnitude Response'); subplot(2,2,3);plot(w/(2*pi),pha/pi);grid on; xlabel('Analog frequency in Hz');ylabel('Phase in pi units'); title('Phase Response'); subplot(2,2,4);plot(t,ha);title('Impulse Response');grid on; xlabel('t (sec) ');ylabel('ha(t) '); Specification: 33. Filter Type: Low Pass Digital(ButterWorth) Using Impulse Invariance Transformation, T=2. 34. PassBand Frequency(Wp) = 0.4*pi. 35. StopBand Frequency(Ws) = 0.6*pi. 36. Passband Ripple= 0.5. 37. Stopband Attunation=50. Wp=0.4*pi ; Ws=0.6*pi ; Rp=0.5;As=50; T=2; OmegaP=Wp/T; OmegaS=Ws/T; [cs,ds] = LP_butterworth (OmegaP,OmegaS,Rp,As); [b,a] = imp_invr(cs,ds,T); [db,mag,pha,grd,w]=freqz_m(b,a); subplot(2,1,1);plot(w/pi,db);grid on; title('Digital ButterWorth Low Pass Filter using Impulse Invariance,T=2'); xlabel('Frequency in Hz');ylabel('Magnitude in Decibel'); Nmax=100;t=0:T/10:Nmax; [ha,x,t]=impulse(cs,ds,t); subplot(2,1,2);plot(t,ha);grid on; [x,n] = impseq(0,0,Nmax);
h=filter(b,a,x); stem(n*T,h); xlabel('time in seconds');ylabel('ha(t) ');
Specification: 38. Filter Type: Low Pass Digital(ButterWorth) Using Billinear Transformation,T=2. 39. PassBand Frequency(Wp) = 0.4*pi. 40. StopBand Frequency(Ws) = 0.6*pi. 41. Passband Ripple= 0.5. 42. Stopband Attunation=50. Wp=0.4*pi ; Ws=0.6*pi ; Rp=0.5;As=50; T=2; OmegaP = (2/T)*tan(Wp/2); OmegaS = (2/T)*tan(Ws/2); [cs,ds] = LP_butterworth (OmegaP,OmegaS,Rp,As); [b,a] = bilinear(cs,ds,1/T); [db,mag,pha,grd,w] = freqz_m(b,a); subplot(2,1,1);plot(w/pi,db);grid on; xlabel('Frequency in Hz');ylabel('Decibel'); title('Digital Low Pass Butterworth filter,Billinear Transformation,T=2'); Nmax=100;t=0:T/10:T*Nmax; [ha,x,t]=impulse(cs,ds,t); subplot(2,1,2);plot(t,ha);grid on; xlabel('time in seconds');ylabel('ha(t) '); title('Impulse Response'); [x,n]=impseq(0,0,Nmax);h=filter(b,a,x); stem(n*T,h);grid on;
