Matlab实现数字FIR的高通_、带通、低通、带阻滤波器的解析

1、利用汉宁窗设计I型数字高通滤波器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; %理想低通滤波器的单位冲激响应 w_ha n=(ha nnin g(N; % 汉宁窗h=hd.*w_ha n; %截取得到实际的单位脉冲响应db,mag,pha,w=freqz_m2(h,1; % 计算实际滤波器的幅度响应 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(理想单位脉冲响应hd(nsubplot(222stem (n, w_ha ntitle(汉宁窗 w(nsubplot(223stem( n,htitle(实际单位脉冲响应hd(nsubplot(224plot(w/pi,dbtitle(幅度响应(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; %理想低通滤波器的单位冲激响应 w_ha n=(ha nnin g(N; % 汉宁窗h=hd.*w_ha n; %截取得到实际的单位脉冲响应db,mag,pha,w=freqz_m2(h,1; % 计算实际滤波器的幅度响应 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(理想单位脉冲响应hd(nsubplot(222stem (n ,w_ha ntitle(汉宁窗 w(nsubplot(223stem( n,htitle(实际单位脉冲响应hd(nsubplot(224plot(w/pi,dbtitle(幅度响应(dBaxis(0,1,-100,10基于切比雪夫一致逼近法设计FIR数字低通滤波器clear all;f=0 0.6 0.7 1; %给定频率轴分点A=1 1 0 0; %给定在这些频率分点上理想的幅频响应weigh=1 10; %给定在这些频率分点上的加权b=remez(32,f,A,weigh; %设计出切比雪夫最佳一

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利用汉宁窗设计I型数字带阻滤波器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; %理想低通滤波器的单位冲激响应w_ha nn=(ha nnin g(N; %汉宁窗h=hd.*w_ha nn; %截取得到实际的单位脉冲响应db,mag,pha,w=freq乙m2(h,1; %计算实际滤波器的幅度响应 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、脉冲响应hd(nsubplot(222stem (n, w_ha nntitle(汉宁窗 w(nsubplot(223stem( n,htitle(实际单位脉冲响应hd(nsubplot(224plot(w/pi,dbtitle(幅度响应(dBaxis(0,1,-100,10利用三角窗设计川型数字带通滤波器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; %理想低通滤波器的单位冲激响应w_tri=(triang(N; % 三角窗h=hd.*w_tri; % 截取得到实际的单位脉冲响应db,mag,pha,w=freqz_im2(h,1; % 计算实际滤波器的幅度响应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、(理想单位脉冲响应hd(nsubplot(222stem( n, w_trititle(三角窗 w(nsubplot(223stem( n,htitle(实际单位脉冲响应hd(nsubplot(224plot(w/pi,dbtitle(幅度响应(dBaxis(0,1,-100,10利用布拉克曼窗设计u型数字带通滤波器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; %理想低通滤波器的单位冲激响应 w_bma n=(blackma n(N; %布拉克曼窗h=hd.*w_bma n; %截取得到实际的单位脉冲响应db,mag,pha,w=freq乙m2(h,1; %计算实际滤波器的幅度响应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(理想单位脉冲响应hd(nsubplot(222stem (n ,w_bma ntitle(布拉克曼窗w(nsubplot(223stem( n,htitle(实际单位脉冲响应hd(nsubplot(224plot(w/pi,dbtitle(幅度响应(dBaxis(0,1,-100,10利用海明窗设计u型数字低通滤波器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; %理想低通滤波器的单位冲激响应 w_ham=(ham min g(N; % 海明窗h=hd.*w_ham; %截取得到实际的单位脉冲响应db,mag,pha,w=freq乙m2(h,1; %计算实际滤波器的幅度响应 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(理想单位脉冲响应hd(nsubplot(222stem (n ,w_hamtitle(海明窗 w(

13、nsubplot(223stem( n,htitle(实际单位脉冲响应hd(nsubplot(224plot(w/pi,dbtitle(幅度响应(dBaxis(0,1,-100,10%fun ctio ndb,mag,pha,w=freqz_m2(b,a%滤波器的幅值响应(相对、绝对、相位响应%db:相对幅值响应%mag绝对幅值响应%pha:相位响应%w采样频率；%b系统函数H(z的分子项(对FIR,b=h%a系统函数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=2

14、0*log10(mag+eps/max(mag; %相对幅值响应 pha=a ngle(H; %相位响应利用模拟Butterworth滤波器设计数字低通滤波器% exa4-8_pulseDF for example4-8% using Butterworth analog lowpass filter to design digital lowpass filter %利用模拟Butterworth滤波器设计数字低通滤波器%脉冲响应不变法wp=0.2*pi;ws=0.3*pi;Rp=1;As=15;T=1;%性能指标Rip=10A(-Rp/20;At n=10A(-As/20;OmgP=wp*

15、T;OmgS=ws*T;N,OmgC=buttord(OmgP,OmgS,Rp,As,s;%选取模拟滤波器的阶数 cs,ds=butter(N,OmgC,s; %设计出所需的模拟低通滤波器b,a=impinvar(cs,ds,T; % 应用脉冲响应不变法进行转换 求得相对、绝对频响及相位、群迟延响应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%滤波器幅值响应(绝对、相对、相位响应及群延迟%Usage: db,mag,pha,grd,w=freqz

18、_m(b,a %500点 对应0,pi%db相对幅值响应;mag绝对幅值响应;pha相位响应;grd群延迟响应%w采样频率;b系统函数H(z的分子项(对FIR,b=h%a系统函数H(z的分母项(对FIR,a=1H,w=freqz(b,a,500;%500 点的复频响应mag=abs(H;db=20*log10(mag+eps/max(mag;pha=a ngle(H;grd=grpdelay(b,a,w;基于频域抽样法的FIR数字带阻滤波器设计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;理想振幅采样响应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; % 构成 Hd(kh仁 ifft(Hdk,N;n=0:1:N-1;h=real(h1.*exp(j*pi*n/N; % 实际单位冲激响应db,mag

20、,pha,w=freqz_m2(h,1;Hr,ww,a,L=hr_type3(h; % 实际振幅响应 subplot(221plot(wl/pi+1/N,Hrs,.,wdl,Hdrtitle(频率样本 Hd(k :N=41axis(0 1 -0.1 1.2subplot(222stem(l,htitle(实际单位脉冲响应h(nsubplot(223plot(ww/pi,Hr,wl/pi+1/N,Hrs,.title(实际振幅响应H(waxis(0 1 -0.1 1.2subplot(224plot(w/pi,dbtitle(幅度响应(dBaxis(O 1-80 10fun ctio n db,

21、mag,pha,w = freqz_m(b,a;%滤波器的幅值响应(相对、绝对、相位响应%db:相对幅值响应%mag绝对幅值响应%pha:相位响应%w采样频率；%b系统函数H(z的分子项(对FIR,b=h%a系统函数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;%计算所设计的3

22、型滤波器的振幅响应%Hr=振幅响应%b=3型滤波器的系数%L=Hr的阶次%h=3型滤波器的单位冲击响应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数字带通滤波器设计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;% 抽样点数al=(M-1/2;wl=(2*pi/M; % 抽样间隔k=0:

23、M-1;T仁0.12; T2=0.6;%过渡带样本点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;%冲击响应图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数字高通滤波器设计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; %理想高通滤波器的单位冲激响应 w_bla=(blackma n(N; % 布拉克曼h=hd.*w_bla; %截取得到实际的单位脉冲响应db,mag,pha,grd,w=freqz_m(h,1; %计算实际滤波器的幅度响应 delta_w=2*pi/1000;As=-round(max(db(1:1:Ws/delta_w+1 %实际阻带纹波,round 是取整函数 y=filter(h,1,xsubplot(246plot(ytitle(滤波后的信号；axis(0,100,-

28、1,1subplot(241stem( n,hdtitle(理想单位脉冲响应hd(nsubplot(242stem( n, w_blatitle(布拉克满窗w(nsubplot(243stem( n,htitle(实际单位脉冲响应hd(nsubplot(244plot(w/pi,dbtitle(幅度响应(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(相频响应function hd=ideal_hp1(Wc,N alp=(N-1/2