Matlab實現(xiàn)數(shù)字FIR的高通_、帶通、低通、帶阻濾波器的解析

1、利用漢寧窗設(shè)計I型數(shù)字高通濾波器clear all;Wp=0.6*pi;Ws=0.4*pi;tr_width=Wp-Ws; % 過渡帶寬度N=ceil(6.2*pi/tr_width % 濾波器長度n=0:1:N-1;Wc=(Ws+Wp/2; %理想低通濾波器的截止頻率hd=ideal_hp1(Wc,N; %理想低通濾波器的單位沖激響應(yīng) w_ha n=(ha nnin g(N; % 漢寧窗h=hd.*w_ha n; %截取得到實際的單位脈沖響應(yīng)db,mag,pha,w=freqz_m2(h,1; % 計算實際濾波器的幅度響應(yīng) delta_w=2*pi/1000;Ap=-(mi n(db(Wp/

2、delta_w+1:1:501 % 實際通帶紋波As=-rou nd(max(db(1:1:Ws/delta_w+1 % 實際阻帶紋波subplot(221stem( n,hdtitle(理想單位脈沖響應(yīng)hd(nsubplot(222stem (n, w_ha ntitle(漢寧窗 w(nsubplot(223stem( n,htitle(實際單位脈沖響應(yīng)hd(nsubplot(224plot(w/pi,dbtitle(幅度響應(yīng)(dBaxis(0,1,-100,10clear all;Wp=0.6*pi;Ws=0.4*pi;tr_width=Wp-Ws; % 過渡帶寬度N=ceil(6.2*p

3、i/tr_width % 濾波器長度n=0:1:N-1;Wc=(Ws+Wp/2; %理想低通濾波器的截止頻率hd=ideal_hp1(Wc,N; %理想低通濾波器的單位沖激響應(yīng) w_ha n=(ha nnin g(N; % 漢寧窗h=hd.*w_ha n; %截取得到實際的單位脈沖響應(yīng)db,mag,pha,w=freqz_m2(h,1; % 計算實際濾波器的幅度響應(yīng) delta_w=2*pi/1000;Ap=-(mi n(db(Wp/delta_w+1:1:501 % 實際通帶紋波As=-rou nd(max(db(1:1:Ws/delta_w+1 % 實際阻帶紋波 subplot(221 s

4、tem( n,hdtitle(理想單位脈沖響應(yīng)hd(nsubplot(222stem (n ,w_ha ntitle(漢寧窗 w(nsubplot(223stem( n,htitle(實際單位脈沖響應(yīng)hd(nsubplot(224plot(w/pi,dbtitle(幅度響應(yīng)(dBaxis(0,1,-100,10基于切比雪夫一致逼近法設(shè)計FIR數(shù)字低通濾波器clear all;f=0 0.6 0.7 1; %給定頻率軸分點(diǎn)A=1 1 0 0; %給定在這些頻率分點(diǎn)上理想的幅頻響應(yīng)weigh=1 10; %給定在這些頻率分點(diǎn)上的加權(quán)b=remez(32,f,A,weigh; %設(shè)計出切比雪夫最佳一

5、致逼近濾波器h,w=freqz(b,1,256,1;h=abs(h;h=20*log10(h;subplot(211stem(b,.;grid;title(切比雪夫逼近濾波器的抽樣值subplot(212plot(w,h;grid;title(濾波器幅頻特性(dB利用漢寧窗設(shè)計I型數(shù)字帶阻濾波器clear all;Wpl=0.2*pi;Wph=0.8*pi;Wsl=0.4*pi;Wsh=0.6*pi;tr_width=min(Wsl-Wpl,(Wph-Wsh; % 過渡帶寬度N=ceil(6.2*pi/tr_width % 濾波器長度n=0:1:N-1;Wcl=(Wsl+Wpl/2; %理想低

6、通濾波器的截止頻率Wch=(Wsh+Wph/2;hd=ideal_bs(Wcl,Wch,N; %理想低通濾波器的單位沖激響應(yīng)w_ha nn=(ha nnin g(N; %漢寧窗h=hd.*w_ha nn; %截取得到實際的單位脈沖響應(yīng)db,mag,pha,w=freq乙m2(h,1; %計算實際濾波器的幅度響應(yīng) delta_w=2*pi/1000;Ap=-(min(db(1:1:Wpl/delta_w+1 % 實際通帶紋波As=-rou nd(max(db(Wsl/delta_w+1:1:Wsh/delta_w+1 % 實際阻帶紋波subplot(221stem( n,hdtitle(理想單位

7、脈沖響應(yīng)hd(nsubplot(222stem (n, w_ha nntitle(漢寧窗 w(nsubplot(223stem( n,htitle(實際單位脈沖響應(yīng)hd(nsubplot(224plot(w/pi,dbtitle(幅度響應(yīng)(dBaxis(0,1,-100,10利用三角窗設(shè)計川型數(shù)字帶通濾波器clear all;Wpl=0.4*pi;Wph=0.6*pi;Wsl=0.2*pi;Wsh=0.8*pi;tr_width=min(Wpl-Wsl,(Wsh-Wph; % 過渡帶寬度N=ceil(6.1*pi/tr_width % 濾波器長度n=0:1:N-1;Wcl=(Wsl+Wpl/2

8、; %理想低通濾波器的截止頻率Wch=(Wsh+Wph/2;hd=ideal_bp2(Wcl,Wch,N; %理想低通濾波器的單位沖激響應(yīng)w_tri=(triang(N; % 三角窗h=hd.*w_tri; % 截取得到實際的單位脈沖響應(yīng)db,mag,pha,w=freqz_im2(h,1; % 計算實際濾波器的幅度響應(yīng)delta_w=2*pi/1000;Ap=-(min(db(Wpl/delta_w+1:1:Wph/delta_w+1 % 實際通帶紋波As=-rou nd(max(db(Wsh/delta_w+1:1:501 %實際阻帶紋波subplot(221stem( n,hdtitle

9、(理想單位脈沖響應(yīng)hd(nsubplot(222stem( n, w_trititle(三角窗 w(nsubplot(223stem( n,htitle(實際單位脈沖響應(yīng)hd(nsubplot(224plot(w/pi,dbtitle(幅度響應(yīng)(dBaxis(0,1,-100,10利用布拉克曼窗設(shè)計u型數(shù)字帶通濾波器clear all;Wpl=0.4*pi;Wph=0.6*pi;Wsl=0.2*pi;Wsh=0.8*pi;tr_width=min(Wpl-Wsl,(Wsh-Wph; % 過渡帶寬度N=ceil(11*pi/tr_width+1 % 濾波器長度n=0:1:N-1;Wcl=(Wsl

10、+Wpl/2; %理想低通濾波器的截止頻率Wch=(Wsh+Wph/2;hd=ideal_bp1(Wcl,Wch,N; %理想低通濾波器的單位沖激響應(yīng) w_bma n=(blackma n(N; %布拉克曼窗h=hd.*w_bma n; %截取得到實際的單位脈沖響應(yīng)db,mag,pha,w=freq乙m2(h,1; %計算實際濾波器的幅度響應(yīng)delta_w=2*pi/1000;Ap=-(min(db(Wpl/delta_w+1:1:Wph/delta_w+1 % 實際通帶紋波As=-rou nd(max(db(Wsh/delta_w+1:1:501 %實際阻帶紋波 subplot(221ste

11、m( n,hdtitle(理想單位脈沖響應(yīng)hd(nsubplot(222stem (n ,w_bma ntitle(布拉克曼窗w(nsubplot(223stem( n,htitle(實際單位脈沖響應(yīng)hd(nsubplot(224plot(w/pi,dbtitle(幅度響應(yīng)(dBaxis(0,1,-100,10利用海明窗設(shè)計u型數(shù)字低通濾波器clear all;Wp=0.2*pi;Ws=0.4*pi;tr_width=Ws-Wp; % 過渡帶寬度N=ceil(6.6*pi/tr_width+1 % 濾波器長度 n=0:1:N-1;Wc=(Ws+Wp/2; %理想低通濾波器的截止頻率 hd=id

12、eal_lp1(Wc,N; %理想低通濾波器的單位沖激響應(yīng) w_ham=(ham min g(N; % 海明窗h=hd.*w_ham; %截取得到實際的單位脈沖響應(yīng)db,mag,pha,w=freq乙m2(h,1; %計算實際濾波器的幅度響應(yīng) delta_w=2*pi/1000;Ap=-(min(db(1:1:Wp/delta_w+1 % 實際通帶紋波As=-rou nd(max(db(Ws/delta_w+1:1:501 % 實際阻帶紋波 subplot(221 stem( n,hdtitle(理想單位脈沖響應(yīng)hd(nsubplot(222stem (n ,w_hamtitle(海明窗 w(

13、nsubplot(223stem( n,htitle(實際單位脈沖響應(yīng)hd(nsubplot(224plot(w/pi,dbtitle(幅度響應(yīng)(dBaxis(0,1,-100,10%fun ctio ndb,mag,pha,w=freqz_m2(b,a%濾波器的幅值響應(yīng)(相對、絕對、相位響應(yīng)%db:相對幅值響應(yīng)%mag絕對幅值響應(yīng)%pha:相位響應(yīng)%w采樣頻率；%b系統(tǒng)函數(shù)H(z的分子項(對FIR,b=h%a系統(tǒng)函數(shù)H(z的分母項(對FIR,a=1H,w=freqz(b,a,1000,whole：H=(H(1:1:501;w=(w(1:1:501;mag=abs(H; %絕對幅值響應(yīng)db=2

14、0*log10(mag+eps/max(mag; %相對幅值響應(yīng) pha=a ngle(H; %相位響應(yīng)利用模擬Butterworth濾波器設(shè)計數(shù)字低通濾波器% exa4-8_pulseDF for example4-8% using Butterworth analog lowpass filter to design digital lowpass filter %利用模擬Butterworth濾波器設(shè)計數(shù)字低通濾波器%脈沖響應(yīng)不變法wp=0.2*pi;ws=0.3*pi;Rp=1;As=15;T=1;%性能指標(biāo)Rip=10A(-Rp/20;At n=10A(-As/20;OmgP=wp*

15、T;OmgS=ws*T;N,OmgC=buttord(OmgP,OmgS,Rp,As,s;%選取模擬濾波器的階數(shù) cs,ds=butter(N,OmgC,s; %設(shè)計出所需的模擬低通濾波器b,a=impinvar(cs,ds,T; % 應(yīng)用脈沖響應(yīng)不變法進(jìn)行轉(zhuǎn)換 求得相對、絕對頻響及相位、群遲延響應(yīng)db,mag,pha,grd,w=freqz_m(b,a;%下面繪出各條曲線subplot(2,2,1;plot(w/pi,mag;title(幅頻特性；xlabel(w(/pi;ylabel(|H(jw|;axis(0,1,0,1.1;set(gca,XTickMode,manual,XTick,

16、0 0.2 0.3 0.5 1;set(gca,YTickMode,ma nual,YTick,0 Atn Rip 1; gridsubplot(2,2,2;plot(w/pi,db;title(幅頻特性(dB;xlabel(w(/pi;ylabel(dB;axis(0,1,-40,5;set(gca,XTickMode,ma nu al,XTick,0 0.2 0.3 0.5 1;set(gca,YTickMode,manual,YTick,-40 -As -Rp 0; gridsubplot(2,2,3;plot(w/pi,pha/pi;title(相頻特性;xlabel(w(/pi;yl

17、abel(pha(/pi;axis(0,1,-1,1;set(gca,XTickMode,ma nual,XTick,0 0.2 0.3 0.5 1;gridsubplot(2,2,4;plot(w/pi,grd;title(群延遲；xlabel(w(/pi;ylabel(Sample;axis(0,1,0,12;set(gca,XTickMode,manual,XTick,0 0.2 0.3 0.5 1;gridfun ctio ndb,mag,pha,grd,w=freqz_m(b,a%濾波器幅值響應(yīng)(絕對、相對、相位響應(yīng)及群延遲%Usage: db,mag,pha,grd,w=freqz

18、_m(b,a %500點(diǎn) 對應(yīng)0,pi%db相對幅值響應(yīng);mag絕對幅值響應(yīng);pha相位響應(yīng);grd群延遲響應(yīng)%w采樣頻率;b系統(tǒng)函數(shù)H(z的分子項(對FIR,b=h%a系統(tǒng)函數(shù)H(z的分母項(對FIR,a=1H,w=freqz(b,a,500;%500 點(diǎn)的復(fù)頻響應(yīng)mag=abs(H;db=20*log10(mag+eps/max(mag;pha=a ngle(H;grd=grpdelay(b,a,w;基于頻域抽樣法的FIR數(shù)字帶阻濾波器設(shè)計clear all;N=41;T仁0.598;alpha= (N-1/2;l=0:N-1;wl= (2*pi/N*l;Hrs=ones(1,6,T1,z

19、eros(1,7,T1,ones(1,11,T1,zeros(1,7,T1,ones(1,6;理想振幅采樣響應(yīng)Hdr=1,1,0,0,1,1;wdl=0,0.3,0.3,0.7,0.7,1;k1=0:floor(N-1/2;k2=floor(N-1/2+1:N-1;an gH=pi/2-alpha*(2*pi/N*(k1+0.5,-pi/2+alpha*(2*pi/N*(N-k2-0.5; %相位約束 條件Hdk=Hrs.*exp(j*angH; % 構(gòu)成 Hd(kh仁 ifft(Hdk,N;n=0:1:N-1;h=real(h1.*exp(j*pi*n/N; % 實際單位沖激響應(yīng)db,mag

20、,pha,w=freqz_m2(h,1;Hr,ww,a,L=hr_type3(h; % 實際振幅響應(yīng) subplot(221plot(wl/pi+1/N,Hrs,.,wdl,Hdrtitle(頻率樣本 Hd(k :N=41axis(0 1 -0.1 1.2subplot(222stem(l,htitle(實際單位脈沖響應(yīng)h(nsubplot(223plot(ww/pi,Hr,wl/pi+1/N,Hrs,.title(實際振幅響應(yīng)H(waxis(0 1 -0.1 1.2subplot(224plot(w/pi,dbtitle(幅度響應(yīng)(dBaxis(O 1-80 10fun ctio n db,

21、mag,pha,w = freqz_m(b,a;%濾波器的幅值響應(yīng)(相對、絕對、相位響應(yīng)%db:相對幅值響應(yīng)%mag絕對幅值響應(yīng)%pha:相位響應(yīng)%w采樣頻率；%b系統(tǒng)函數(shù)H(z的分子項(對FIR,b=h%a系統(tǒng)函數(shù)H(z的分母項(對FIR,a=1H,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 = an gle(H;% pha = un wrap(a ngle(H;function Hr,w,c丄=hr_type3(h;%計算所設(shè)計的3

22、型濾波器的振幅響應(yīng)%Hr=振幅響應(yīng)%b=3型濾波器的系數(shù)%L=Hr的階次%h=3型濾波器的單位沖擊響應(yīng)M=le ngth(h;L=(M-1/2;c= 2*h(L+1:-1:1;n=0:1:L;w=0:1:500*2*pi/500;Hr=si n(w*n *c;基于頻域抽樣法的FIR數(shù)字帶通濾波器設(shè)計wsl=0.12*pi;%低阻帶邊緣wsh=0.82*pi;%高阻帶邊緣wpl=0.32*pi;%低通帶邊緣wph=0.62*pi;%高通帶邊緣delta=(wpl-wsl;% 過度帶M=ceil(2*pi*3/delta;% 抽樣點(diǎn)數(shù)al=(M-1/2;wl=(2*pi/M; % 抽樣間隔k=0:

23、M-1;T仁0.12; T2=0.6;%過渡帶樣本點(diǎn)Hrs=zeros(1,ceil(0.12*pi/wl+1,T2,T1,o nes(1,ceil(0.3*pi/wl,T1,T2,zeros(1,ceil(0.37 34*pi/wl,T2,T1, on es(1,c eil(0.3*pi/wl,T1,T2,zeros(1,ceil(0.12*pi/wl+1;wdl=0 0.12 0.32 0.62 0.82 1;k1=0:floor(M-1/2;k2=floor(M-1/2+1:M-1;an gH=-al*(2*pi/M*k1,al*(2*pi/M*(M-k2;H=Hrs.*exp(j*a

24、ngH;h=real(ifft(H;%傅立葉反變換figure(1;%沖擊響應(yīng)圖stem(k,h;title(impulse resp on se;xlabel( n;ylabel(h( n;grid;figure(2;%幅頻曲線圖Hf=abs(H;w=k*wl/pi;plot(w,Hf,*b-axis(0 1-0.1 1.1;title(amplitude resp on se;xlabel(freque ncy in pi un its;ylabel(Hr(w; set(gca,xtickmode,ma nu al,xtick,wdl; set(gca,ytickmode,manual,y

25、tick,O 0.12 0.6 1; grid; figure;fs=15000;c,f3=freqz(h,1;f3=f3/pi*fs/2;plot(f3,20*log10(abs(c;title(頻譜特性；xlabel(頻率 /HZ;ylabel(衰減 /dB;grid;t=(0:100/fs;x=si n(2*pi*t*700+si n(2*pi*t*3200+si n(2*pi*t*6200; q=filter(h,1,x;a,f1=freqz(x;f1=f1/pi*fs/2;b,f2=freqz(q;f2=f2/pi*fs/2;figure(4;subplot(2,1,1;plot(f

26、1,abs(a;title(輸入波形頻譜圖；xlabel(頻率;ylabel(幅度subplot(2,1,2;plot(f2,abs(b;title(輸出波形頻譜圖；xlabel(頻率;ylabel(幅度基于漢寧窗的FIR數(shù)字高通濾波器設(shè)計fun ctio n s2Fs=15000;t=(0:100/Fs;x=si n(2*pi*500*t+si n(2*pi*3000*tsubplot(245;stem(x;title(原始信號；axis(0,100,-2,2;Ws=7*pi/30;Wp=13*pi/30;tr_wid=Wp-Ws; %過渡帶寬度N=ceil(11*pi/tr_wid % 濾

27、波器長度n=0:1:N-1;Wc=(Ws+Wp/2; %理想高通濾波器的截止頻率hd=ideal_hp1(Wc,N; %理想高通濾波器的單位沖激響應(yīng) w_bla=(blackma n(N; % 布拉克曼h=hd.*w_bla; %截取得到實際的單位脈沖響應(yīng)db,mag,pha,grd,w=freqz_m(h,1; %計算實際濾波器的幅度響應(yīng) delta_w=2*pi/1000;As=-round(max(db(1:1:Ws/delta_w+1 %實際阻帶紋波,round 是取整函數(shù) y=filter(h,1,xsubplot(246plot(ytitle(濾波后的信號；axis(0,100,-

28、1,1subplot(241stem( n,hdtitle(理想單位脈沖響應(yīng)hd(nsubplot(242stem( n, w_blatitle(布拉克滿窗w(nsubplot(243stem( n,htitle(實際單位脈沖響應(yīng)hd(nsubplot(244plot(w/pi,dbtitle(幅度響應(yīng)(dBaxis(0,1,-100,10fun ctio n db,mag,pha,grd,w = freqz_m(b,a;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 = an gle(H;grd = grpdelay(b,a,w; subplot(247; plot(pha title(相頻響應(yīng)function hd=ideal_hp1(Wc,N alp=(N-1/2