MATLAB实现数字FIR的高通_和带通等滤波器的源程序.doc

利用汉宁窗设计型数字高通滤波器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_han=(hanning(N); %汉宁窗h=hd.*w_han; %截取得到实际的单位脉冲响应db,mag,pha,w=freqz_m2(h,1); %计算实际滤波器的幅度响应delta_w=2*pi/1000;Ap=-(min(db(Wp/delta_w+1:1:501) %实际通带纹波As=-round(max(db(1:1:Ws/delta_w+1) %实际阻带纹波subplot(221)stem(n,hd)title(理想单位脉冲响应hd(n)subplot(222)stem(n,w_han)title(汉宁窗w(n)subplot(223)stem(n,h)title(实际单位脉冲响应hd(n)subplot(224)plot(w/pi,db)title(幅度响应(dB)axis(0,1,-100,10)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_han=(hanning(N); %汉宁窗h=hd.*w_han; %截取得到实际的单位脉冲响应db,mag,pha,w=freqz_m2(h,1); %计算实际滤波器的幅度响应delta_w=2*pi/1000;Ap=-(min(db(Wp/delta_w+1:1:501) %实际通带纹波As=-round(max(db(1:1:Ws/delta_w+1) %实际阻带纹波subplot(221)stem(n,hd)title(理想单位脉冲响应hd(n)subplot(222)stem(n,w_han)title(汉宁窗w(n)subplot(223)stem(n,h)title(实际单位脉冲响应hd(n)subplot(224)plot(w/pi,db)title(幅度响应(dB)axis(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); %设计出切比雪夫最佳一致逼近滤波器h,w=freqz(b,1,256,1);h=abs(h);h=20*log10(h);subplot(211)stem(b,.);grid;title(切比雪夫逼近滤波器的抽样值)subplot(212)plot(w,h);grid;title(滤波器幅频特性(dB)利用汉宁窗设计型数字带阻滤波器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; %理想低通滤波器的截止频率Wch=(Wsh+Wph)/2;hd=ideal_bs(Wcl,Wch,N); %理想低通滤波器的单位冲激响应w_hann=(hanning(N); %汉宁窗h=hd.*w_hann; %截取得到实际的单位脉冲响应db,mag,pha,w=freqz_m2(h,1); %计算实际滤波器的幅度响应delta_w=2*pi/1000;Ap=-(min(db(1:1:Wpl/delta_w+1) %实际通带纹波As=-round(max(db(Wsl/delta_w+1:1:Wsh/delta_w+1) %实际阻带纹波subplot(221)stem(n,hd)title(理想单位脉冲响应hd(n)subplot(222)stem(n,w_hann)title(汉宁窗w(n)subplot(223)stem(n,h)title(实际单位脉冲响应hd(n)subplot(224)plot(w/pi,db)title(幅度响应(dB)axis(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; %理想低通滤波器的截止频率Wch=(Wsh+Wph)/2;hd=ideal_bp2(Wcl,Wch,N); %理想低通滤波器的单位冲激响应w_tri=(triang(N); %三角窗h=hd.*w_tri; %截取得到实际的单位脉冲响应db,mag,pha,w=freqz_m2(h,1); %计算实际滤波器的幅度响应delta_w=2*pi/1000;Ap=-(min(db(Wpl/delta_w+1:1:Wph/delta_w+1) %实际通带纹波As=-round(max(db(Wsh/delta_w+1:1:501) %实际阻带纹波subplot(221)stem(n,hd)title(理想单位脉冲响应hd(n)subplot(222)stem(n,w_tri)title(三角窗w(n)subplot(223)stem(n,h)title(实际单位脉冲响应hd(n)subplot(224)plot(w/pi,db)title(幅度响应(dB)axis(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(11*pi/tr_width)+1 %滤波器长度n=0:1:N-1;Wcl=(Wsl+Wpl)/2; %理想低通滤波器的截止频率Wch=(Wsh+Wph)/2;hd=ideal_bp1(Wcl,Wch,N); %理想低通滤波器的单位冲激响应w_bman=(blackman(N); %布拉克曼窗h=hd.*w_bman; %截取得到实际的单位脉冲响应db,mag,pha,w=freqz_m2(h,1); %计算实际滤波器的幅度响应delta_w=2*pi/1000;Ap=-(min(db(Wpl/delta_w+1:1:Wph/delta_w+1) %实际通带纹波As=-round(max(db(Wsh/delta_w+1:1:501) %实际阻带纹波subplot(221)stem(n,hd)title(理想单位脉冲响应hd(n)subplot(222)stem(n,w_bman)title(布拉克曼窗w(n)subplot(223)stem(n,h)title(实际单位脉冲响应hd(n)subplot(224)plot(w/pi,db)title(幅度响应(dB)axis(0,1,-100,10)利用海明窗设计型数字低通滤波器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=ideal_lp1(Wc,N); %理想低通滤波器的单位冲激响应w_ham=(hamming(N); %海明窗h=hd.*w_ham; %截取得到实际的单位脉冲响应db,mag,pha,w=freqz_m2(h,1); %计算实际滤波器的幅度响应delta_w=2*pi/1000;Ap=-(min(db(1:1:Wp/delta_w+1) %实际通带纹波As=-round(max(db(Ws/delta_w+1:1:501) %实际阻带纹波subplot(221)stem(n,hd)title(理想单位脉冲响应hd(n)subplot(222)stem(n,w_ham)title(海明窗w(n)subplot(223)stem(n,h)title(实际单位脉冲响应hd(n)subplot(224)plot(w/pi,db)title(幅度响应(dB)axis(0,1,-100,10)%-functiondb,mag,pha,w=freqz_m2(b,a)%滤波器的幅值响应(相对、绝对)、相位响应%db:相对幅值响应%mag：绝对幅值响应%pha: 相位响应%w 采样频率；%b 系统函数H(z)的分子项(对FIR，b=h)%a 系统函数H(z)的分母项(对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); %相位响应利用模拟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=10(-Rp/20);Atn=10(-As/20);OmgP=wp*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,0 0.2 0.3 0.5 1);set(gca,YTickMode,manual,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,manual,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);ylabel(pha(/pi);axis(0,1,-1,1);set(gca,XTickMode,manual,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);gridfunctiondb,mag,pha,grd,w=freqz_m(b,a)%滤波器幅值响应(绝对、相对)、相位响应及群延迟%Usage: db,mag,pha,grd,w=freqz_m(b,a) %500点对应0,pi%db 相对幅值响应；mag 绝对幅值响应；pha 相位响应；grd 群延迟响应%w 采样频率；b 系统函数H(z)的分子项(对FIR，bh)%a 系统函数H(z)的分母项(对FIR，a1)H,w=freqz(b,a,500);%500点的复频响应mag=abs(H);db=20*log10(mag+eps)/max(mag);pha=angle(H);grd=grpdelay(b,a,w);基于频域抽样法的FIR数字带阻滤波器设计clear all;N=41;T1=0.598;alpha= (N-1)/2;l=0:N-1;wl= (2*pi/N)*l;Hrs=ones(1,6),T1,zeros(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;angH=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(k)h1=ifft(Hdk,N);n=0:1:N-1;h=real(h1.*exp(j*pi*n/N); %实际单位冲激响应db,mag,pha,w=freqz_m2(h,1);Hr,ww,a,L=hr_type3(h); %实际振幅响应subplot(221)plot(wl/pi+1/N,Hrs,.,wdl,Hdr)title(频率样本Hd(k) :N=41)axis(0 1 -0.1 1.2)subplot(222)stem(l,h)title(实际单位脉冲响应h(n)subplot(223)plot(ww/pi,Hr,wl/pi+1/N,Hrs,.)title(实际振幅响应H(w)axis(0 1 -0.1 1.2)subplot(224)plot(w/pi,db)title(幅度响应(dB)axis(0 1 -80 10)function db,mag,pha,w = freqz_m(b,a);%滤波器的幅值响应(相对、绝对)、相位响应%db:相对幅值响应%mag：绝对幅值响应%pha: 相位响应%w 采样频率；%b 系统函数H(z)的分子项(对FIR，b=h)%a 系统函数H(z)的分母项(对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);% pha = unwrap(angle(H);function Hr,w,c,L=hr_type3(h);%计算所设计的3型滤波器的振幅响应%Hr=振幅响应%b=3型滤波器的系数%L=Hr的阶次%h=3型滤波器的单位冲击响应M=length(h);L=(M-1)/2;c= 2*h(L+1:-1:1);n=0:1:L;w=0:1:500*2*pi/500;Hr=sin(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:M-1;T1=0.12; T2=0.6;%过渡带样本点Hrs=zeros(1,ceil(0.12*pi/wl)+1),T2,T1,ones(1,ceil(0.3*pi/wl),T1,T2,zeros(1,ceil(0.3734*pi/wl),T2,T1,ones(1,ceil(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;angH=-al*(2*pi)/M*k1,al*(2*pi)/M*(M-k2);H=Hrs.*exp(j*angH);h=real(ifft(H);%傅立叶反变换figure(1);%冲击响应图stem(k,h);title(impulse response);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 response);xlabel(frequency in pi units);ylabel(Hr(w);set(gca,xtickmode,manual,xtick,wdl);set(gca,ytickmode,manual,ytick,0 0.12 0.6 1);grid;figure(3);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=sin(2*pi*t*700)+sin(2*pi*t*3200)+sin(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(f1,abs(a);title(输入波形频谱图);xlabel(频率);ylabel(幅度)subplot(2,1,2);plot(f2,abs(b);title(输出波形频谱图);xlabel(频率);ylabel(幅度)基于汉宁窗的FIR数字高通滤波器设计function s2Fs=15000;t=(0:100)/Fs;x=sin(2*pi*500*t)+sin(2*pi*3000*t)subplot(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) %滤波器长度n=0:1:N-1;Wc=(Ws+Wp)/2; %理想高通滤波器的截止频率hd=ideal_hp1(Wc,N); %理想高通滤波器的单位冲激响应w_bla=(blackman(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,x)subplot(246)plot(y)title(滤波后的信号);axis(0,100,-1,1)subplot(241)stem(n,hd)title(理想单位脉冲响应hd(n)subplot(242)stem(n,w_bla)title(布拉克满窗w(n)subplot(243)stem(n,h)title(实际单位脉冲响应hd(n)subplot(244)plot(w/pi,db)title(幅度响应(dB)axis(0,1,-100,10)function 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 = angle(H);grd = grpdelay(b,a,w);subplot(247);plot(pha)title(相频响应)function hd=ideal_hp1(Wc,N)alp=(N-1)/2