基于MATLAB的布莱克曼窗FIR数字低通滤波器设计程序

程序：clearall;Wp=1.34*pi/3;%input通带频率;Ws=7.4*pi/30;%input阻带频率;Ar=75;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_ham=(blackman(N))';string=['blackman','N=',num2str(N)];%调用布莱克曼窗函数h=hd.*w_ham;%截取得到实际的单位冲击响应，为理想脉冲响应与窗函数乘积[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:501)));%实际阻带波纹subplot(4,2,1);stem(n,hd);title('理想脉冲响应hd(n)');axis([0,N-1,-0.5,1]);xlabel('n');ylabel('hd(n)')subplot(4,2,2);stem(n,w_ham);axis([0,N-1,0,1.1]);xlabel('n');ylabel('w(n)');text(1.5,1.3,string);subplot(4,2,3);stem(n,h);title('实际脉冲响应h(b)');axis([0,N,-1.4,1.4]);xlabel('n');ylabel('h(n)');subplot(4,2,4);plot(w,pha);title('相频特性');axis([0,3.15,-4,4]);xlabel('频率(pi)');ylabel('相位(Φ)');subplot(4,2,5);plot(w/pi,db);title('幅度特性(dB)');axis([0,1,-100,10]);xlabel('频率(pi)');ylabel('(分贝数)');subplot(4,2,6);plot(w,mag);title('频率特性');axis([0,3.15,0,1.5]);xlabel('频率(rad)');ylabel('幅值');调用子函数1：ideal_hp1functionhd=ideal_hp1(Wc,N)alpha=(N-1)/2;n=0:1:N-1;m=n-alpha+eps;hd=[sin(pi*m)-sin(Wc*m)]./(pi*m);%全通-低通；调用子函数2：ideal_hp1%计算z域频率响应；function[db,mag,pha,w]=freqz_m2(b,a)%频率的取值范围在0—1之间，标准化频率1对应%数字频率为π，对应的模拟频率为采样频率的一半；[H,w]=freqz(b,a,1000,'whole');%b直接型滤波器分子系数;a直接型滤波器分母系数%db为[0:pi]之间的相对振幅;mag为[0pi]区间的绝对值;pha为[0pi]区间的相位响应；H=(H(1:1:501))';w=(w(1:1:501))';mag=abs(H);%幅度db=20*log10((mag+eps)/max(mag));pha=angle(H);%相位%用含多种频率成份离散时间信号作输入信号，得到通过滤波器后与输出信号的对比结果fs=15000;t=(0:100)/fs;x=sin(2*pi*t*750)+sin(2*pi*t*1850)+sin(2*pi*t*2500)+sin(2*pi*t*3350)+sin(2*pi*t*6100);q=filter(h,1,x);[a,f1]=freqz(x);f1=f1/pi*fs/2;[b,f2]=freqz(q);f2=f2/pi*fs/2;subplot(4,2,7);plot(f1,abs(a));title('输入波形频谱图');xlabel('频率/Hz');ylabel('幅度')subplot(4,2,8);plot(f2,abs(b));title('输出波形频谱图');xlabel('频率/Hz');ylabel('幅度')%computationofs-domainfrequencyresponse:modifiedversion%--------------------------------------------------------%[dB,mag,pha,grd,w]=freqz_m(b,a);%dB=RelativemagnitudeindBover[0to2pi]%mag=Absolutemagnitudeover[0to2pi]%pha=Phaseresponseinradiansover[0to2pi]%w=arrayof500frequenysamplesbetween[0to2pi]%b=numeratorpolynomialcoefficientsofH(z)%a=denominatorpolynomialcoefficientsofH(z)%wrangingfrom0to2piradianspersample.%uses1000samplepointsaroundtheentireunitcircletocalculatethefrequencyresponse.%z域数字频率响应计算freqz的改进版本%b直接型滤波器分子系数;a直接型滤波器分母系数%db为[0:pi]之间的相对振幅;mag为[0pi]区间的绝对值;pha为[0pi]区间的相位响应;%grd为[0pi]区间的群延迟;w为[0pi]区间内的501个频率样本数组freqz_m的定义：function[db,mag,pha,grd,w]=freqz_m(b,a)%求取系统的绝对幅度响应、相对的db值幅度响应、相位响应和群延时响应的函数%%db为相对振幅(dB)%mag为绝对振幅%pha为相位响应%grd为群延时%w为频率样本点向量%b为Ha(z)分子多项式系数(对FIR而言，b=h)%a为Hz(z)分母多项式系数(对FIR而言，a=1)%[H,w]=freqz(b,a,1000,'whole');%freqz显示数字滤波器频域中的图形%[H,W]=FREQZ(B,A,N,'whole')usesNpointsaroundthewholeunitcircle.H=(H(1:501))';w=(w(1:501))';mag=abs(H);db=20*log10((mag+eps)/max(mag));%dB=-20log10|H(ejw)|/|H(ejw)|max=0pha=angle(H);grd=grpdelay(b,a,w);frepz_m2的定义为：function[db,mag,pha,w]=freqz_m2(b,a)%滤波器幅值响应（绝对、相对）、相位响应%db:相对幅值响应；%mag:绝对幅值响应；%pha:相位响应；%w:采样频率；%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);fs=15000;t=(0:100)/fs;x=sin(2*pi*t*750)+sin(2*pi*t*2500)+sin(2*pi*t*6100);q=filter(h,1,x);[a,f1]=freqz(x);f1=f1/pi*fs/2;[b,f2]=freqz(q);f2=f2/pi*fs/2;subplot(4,2,7);plot(f1,abs(a));title('输入波形频谱图');xlabel('频率/Hz');ylabel('幅度')subplot(4,2,8);plot(f2,abs(b));title('输出波形频谱图');xlabel('频率/Hz');ylabel('幅度')程序二例2、设计一个FIR数字高通滤波器，使其满足：dB60,4.0dB5.0,6.0ssppAR[MATLAB程序]：%designadigitalFIRhighpassfilterwithBlackmanwindowws=0.4*pi;wp=0.6*pi;tr_width=wp-ws;%过虑带宽度%N=ceil(11*pi/tr_width)+1;N=ceil(11*pi/tr_width);%滤波器长度n=[0:1:N-1];wc=(ws+wp)/2;%理想滤波器的截止频率hd=ideal_lp(pi,N)-ideal_lp(wc,N);%理想滤波器的单位冲击响应w_black=(blackman(N))';h=hd.*w_black;%截取得到的实际的单位脉冲响应[dB,mag,pha,grd,w]=freqz_m(h,[1]);%计算实际滤波器的单位脉冲响应delta_w=2*pi/1000;Rp=-(min(dB(wp/delta_w+1:1:501)));%实际通带波纹As=-round(max(dB(1:1:ws/delta_w+1)));%实际阻带波纹%plotssubplot(2,2,1);stem(n,hd);title('理想脉冲响应');axis([0,N-1,-0.4,0.6]);xlabel('n');%set(gca,'YTickMode','manual','YTick',[0,0.1,0.2,0.3]);subplot(2,2,2);stem(n,w_black);title('BlackmanWindow');axis([0,N-1,0,1.1]);xlabel('n');subplot(2,2,3);stem(n,h);title('实际脉冲响应');axis([0,N-1,-0.4,0.6]);xlabel('n');%set(gca,'YTickMode','manual','YTick',[0,0.1,0.2,0.3]);subplot(2,2,4);plot(w/pi,dB);title('幅值响应indB');grid;axis([0,1,-120,0]);xlabel('pi');ylabel('dB');%set(gca,'XTickMode','manual','XTick',[0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0]);%set(gca,'YTickMode','manual','YTick',[-100,-90,-80,-70,-60,-50,-40,-30,-20,-10,0]);[运行结果]：%计算z域频率响应function[dB,mag,pha,grd,w]=freqz_m(b,a)%computationofs-domainfrequencyresponse:modifiedversion%--------------------------------------------------------%[dB,mag,pha,grd,w]=freqz_m(b,a);%dB=RelativemagnitudeindBover[0to2pi]%mag=Absolutemagnitudeover[0to2pi]%pha=Phaseresponseinradiansover[0to2pi]%w=arrayof500frequenysamplesbetween[0to2pi]%b=numeratorpolynomialcoefficientsofH(z)%a=denomin