学习和掌握用MATLAB程序实现低通、高通、带通和带阻滤波器算法

《数字信号处理》实验报告 实验六  四种数字滤波器设计 班级：计科121              学号：1208060135 姓名： 刘国强                成绩： 日期：2014年11月17日    地点：博学楼 706 一、实验目的 1. 学习和掌握用MATLAB程序实现低通、高通、带通和带阻滤波器算法 2. 用MATLAB实现上述四种滤波器。 二、实验内容 1.用MATLAB程序实现低通、高通、带通和带阻滤波器算法 2.具体实现方法自己确定，但必须有对应的数学公式、 3. 关键参数δp、δs、ωp、 ωs、ωc和N的设定。 4. 探讨N值大小对滤波图曲线影响，并作出该滤波器幅频和相频特性曲线图。 三、实验结果 基于切比雪夫一致逼近法设计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)') 基于汉宁窗的FIR数字高通滤波器设计 function s2 Fs=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; n=0:1:N-1; m=n-alp+eps;    %eps是一个很小很小的数 hd=[sin(pi*m)-sin(Wc*m)]./(pi*m); 基于频域抽样法的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);                %幅频曲线图frequency in pi units Hf=abs(H);w=k*wl/pi; plot(w,Hf,'*b-') axis([0 1 -0.1 1.1]); title(‘impulse  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('幅度') 用matlab设计一个带阻滤波器，信号有5Hz，40Hz，50Hz，要求把40Hz的滤除，程序是这样的，滤波器和产生的信 w1=38/100; w2=42/100; [B,A]=butter(1,[w1,w2],'stop'); [h,w]=freqz(B,A); f=w/pi*100; figure(1), plot(f,20*log10(abs(h))); axis([0,80,-30,10]); grid; xlabel('频率/Hz') ylabel('幅度/dB') title('频率响应曲线') fs=200; t=0:1/fs:1; x1=sin(2*pi*5*t)+sin(2*pi*50*t)+sin(2*pi*40*t); Y1=fft(x1,128); Y2=abs(Y1); f=(0:64)*100/64; figure(2), subplot(2,1,1),plot(f,Y2(1:65)),title('原始信号'); y=filter(B,A,Y2); subplot(2,1,2),plot(y),title('滤波后信号'); 四、实验感想和建议 通过本次的实验，学会了用MATLAB程序实现低通、高通、带通和带阻滤波器算法，并且实现了关键参数δp、δs、ωp、 ωs、ωc和N的设定。 但是未探讨N值大小对滤波图曲线影响，也未作出该滤波器幅频和相频特性曲线图。 继续阅读