基于MATLAB信号的滤波处理学术论文.doc

基于MATLAB信号的滤波处理吴彦亨，张鹏飞，谢德力，张龙摘要:滤波器设计在数字信号处理中占有极其重要的地位，FIR数字滤波器和IIR滤波器是滤波器设计的重要组成部分。利用MATLAB信号处理工具箱可以快速有效地设计各种数字滤波器。课题基于MATLAB动态信号滤波处理的设计与实现，综合运用数字信号处理的理论知识对动态信号进行时域、频域分析和滤波。通过理论推导得出相应结论，再利用 MATLAB 作为编程工具进行计算机实现。在设计实现的过程中，使用窗函数法来设计FIR数字滤波器，用巴特沃斯、切比雪夫和双线性变法设计IIR数字滤波器，并利用MATLAB 作为辅助工具完成设计中的计算与图形的绘制。通过对对所设计滤波器的仿真和频率特性分析，可知利用MATLAB信号处理工具箱可以有效快捷地设计FIR和IIR数字滤波器，过程简单方便，结果的各项性能指标均达到指定要求。关键词:哈明窗 汉宁窗 巴特沃斯 切比雪夫 双线性变换 频谱分析1 绪论数字信号处理是利用计算机或专用处理设备，以数值计算的方法对信号进行采集、抽样、变换、综合、估值与识别等加工处理，借以达到提取信息和便于应用的目的。它在语音、雷达、图像、系统控制、通信、航空航天、生物医学等众多领域都获得了极其广泛的应用。具有灵活、精确、抗干扰强、设备尺寸小、造价低、速度快等优点。数字滤波器, 是数字信号处理中及其重要的一部分。随着信息时代和数字技术的发展，受到人们越来越多的重视。数字滤波器可以通过数值运算实现滤波，所以数字滤波器处理精度高、稳定、体积小、重量轻、灵活不存在阻抗匹配问题，可以实现模拟滤波器无法实现的特殊功能。数字滤波器种类很多，根据其实现的网络结构或者其冲激响应函数的时域特性，可分为两种，即有限冲激响应( FIR，Finite Impulse Response)滤波器和无限冲激响应( IIR，Infinite Impulse Response)滤波器。FIR滤波器结构上主要是非递归结构，没有输出到输入的反馈，系统函数H (z)在处收敛，极点全部在z = 0处（因果系统），因而只能用较高的阶数达到高的选择性。FIR数字滤波器的幅频特性精度较之于IIR数字滤波器低，但是线性相位，就是不同频率分量的信号经过fir滤波器后他们的时间差不变，这是很好的性质。FIR数字滤波器是有限的单位响应也有利于对数字信号的处理，便于编程，用于计算的时延也小，这对实时的信号处理很重要。 FIR滤波器因具有系统稳定，易实现相位控制，允许设计多通带（或多阻带）滤波器等优点收到人们的青睐。IIR滤波器采用递归型结构，即结构上带有反馈环路。IIR滤波器运算结构通常由延时、乘以系数和相加等基本运算组成，可以组合成直接型、正准型、级联型、并联型四种结构形式，都具有反馈回路。同时，IIR数字滤波器在设计上可以借助成熟的模拟滤波器的成果，如巴特沃斯、契比雪夫和椭圆滤波器等，有现成的设计数据或图表可查，在设计一个IIR数字滤波器时，我们根据指标先写出模拟滤波器的公式，然后通过一定的变换，将模拟滤波器的公式转换成数字滤波器的公式。滤波器的设计可以通过软件或设计专用的硬件两种方式来实现。随着MATLAB软件及信号处理工具箱的不断完善，MATLAB很快成为应用学科等领域不可或缺的基础软件。它可以快速有效地实现数字滤波器的设计、分析和仿真，极大地减轻了工作量,有利于滤波器设计的最优化。2设计内容设计一个工作流程如图3所示的信号滤波演示系统：图1 信号滤波演示系统工作流程数字滤波信号发生器IIR DF(双线型变换法)信号选择FIR DF(窗函数法)滤波后信号滤波器类型选择(IIR、FIR)DFT频谱分析滤波器技术指标基型选择窗口选择滤波器设计测试信号滤波器系数输出信号分析3设计的具体实现3.1 动态信号的实现利用PC 机上的Matlab可以进行输入由多个不同频率正弦信号叠加组合而成的模拟信号公式、采样频率（Hz）以及采样点数，动态生成该信号的采样序列，作为测试信号。3. 2动态信号频谱分析在本次设计中，我们是利用MATLAB生成的模拟信号并对其频谱分析。动态信号频谱分析的主要程序如下：fs=input(频率等于：);%fs=str2num(fs);m=1/fs;f0=100;t=-0.2:m:0.2;sum=0;for n=1:2:9;fn=f0+n*10;hold on;endfor n=1:1:9;sum=sum+100*sin(2*pi*fn*t);x1=sum;endt1=-0.2:m:0.2;y=ones(1,length(t1);t1=-0.2:m:0.2;val=get(hObject,Value);str=get(hObject,String);switch strval case 采样函数%User selectstitle(原始信号)t1=-0.2:m:0.2;y=ones(1,length(t1);stem(t1,y,.),gridaxis(-0.2,0.2,-1000,1000);title(采样信号)t1=-0.2:m:0.2;stem(t1,sum,.),gridtitle(采样后信号)grid on case DFT变换的幅频函数%User selects N=input(请输入N的值 N:);M=input(请输入M的值 M:); X1=fft(sum,M);X2=fft(sum,N);k=0:M-1;wk=2*k/M;plot(wk,abs(X1);title(a)FTx1的幅频特性图);xlabel(w/pi);ylabel(幅度);grid on caseDFT变换的相频函数%User selectsN=input(请输入N的值 N:);M=input(请输入M的值 M:); X1=fft(sum,M);X2=fft(sum,N);k=0:M-1;wk=2*k/M; plot(wk,angle(X1);title(b)FTx1的相频特性图);xlabel(w/pi);ylabel(相位); grid on end图1：采样函数的图形图形2：幅频函数的图形3. 设计FIR和IIR数字滤波器IIR滤波器和FIR滤波器的设计方法完全不同。IIR滤波器设计方法有间接法和直接法，间接法是借助于模拟滤波器的设计方法进行的。其设计步骤是：先设计过渡模拟滤波器得到系统函数H（s），然后将H（s）按某种方法转换成数字滤波器的系统函数H(z)。FIR滤波器比鞥采用间接法，常用的方法有窗函数法、频率采样发和切比雪夫等波纹逼近法。对于线性相位滤波器，经常采用FIR滤波器。 对于数字高通、带通滤波器的设计，通用方法为双线性变换法。可以借助MATLAB中的函数库进行调用程序所要的步骤如下：（1）确定所需类型数字滤波器的技术指标。（2）将所需类型数字滤波器的边界频率转换成相应的模拟滤波器的边界频率，转换公式为=2/T tan(0.5)(3)将相应类型的模拟滤波器技术指标转换成模拟低通滤波器技术指标。（4）调用设计其他滤波器的函数以实现其他滤波器的设计 或者将低通滤波器通过”BH,AH=lp2ph(B,A,wph)”进行低通转高通如前所述，应用于IIR滤波器的设计对于FIR滤波器的设计还是有区别的、一下主要是FIR滤波器的设计步骤： （1）根据对阻带衰减及过渡带的指标要求，选择串窗数类型（矩形窗、三角窗、汉宁窗、哈明窗、凯塞窗等），并估计窗口长度N。先按照阻带衰减选择窗函数类型。原则是在保证阻带衰减满足要求的情况下，尽量选择主瓣的窗函数。（2）构造希望逼近的频率响应函数。（3）计算h(n).。（4）加窗得到设计结果。对于IIR和FIR的滤波器而言，FIR的最大优点在于可以实现线性相位的滤波，因此在这两种不同的滤波器滤波也会使所还原的原信号的失真程度不同，一下介绍IIR与FIR几种特殊的滤波器滤波的比较及其程序：一 IIR低通滤波器 （1）巴特沃斯低通滤波器Fs=input(采样频率等于（Hz）：);rp=input(通带最大衰减（dB）：);rs=input(阻带最小衰减（dB）：);c=input(通带边界频率（Hz）：);d=input(阻带边界频率（Hz）：);T=1/Fs;wp=(2/T)*tan(2*pi*c/Fs)/2);ws=(2/T)*tan(2*pi*d/Fs)/2);N,wc=buttord(wp,ws,rp,rs,s)B,A=butter(N,wc,s);Bz,Az=bilinear(B,A,Fs)wk=0:pi/512:pi;Hz,w=freqz(Bz,Az,wk);Hx=angle(Hz);subplot(2,1,1);plot(w,Hx); hold off;xlabel(omega/pi);ylabel(相位);title(a);subplot(2,1,2);hold off;plot(wk/pi,20*log10(abs(Hz);xlabel(omega/pi);ylabel(幅度（dB）);axis(0,1,-100,5);title(b);hold off;图3 巴特沃斯低通滤波器（2）巴特沃斯带通滤波器wp=(2/T)*tan(2*pi*c/Fs)/2);ws=(2/T)*tan(2*pi*d/Fs)/2); N,wc=buttord(wp,ws,rp,rs,s)B,A=butter(N,wc,high,s);Bz,Az=bilinear(B,A,Fs)wk=0:pi/512:pi;Hz,w=freqz(Bz,Az,wk);Hx=angle(Hz);subplot(2,1,1);plot(w,Hx);grid on;xlabel(omega/pi);ylabel(相位);title(相频特性);subplot(2,1,2);plot(wk/pi,20*log10(abs(Hz);grid on;xlabel(omega/pi);ylabel(幅度(dB);axis(0,1,-100,5);title(幅频特性); grid on图4：巴特沃高通滤波器（3）切比雪夫2带通滤波器Fs=input(采样频率等于（Hz）：);rp=input(通带最大衰减（dB）：);rs=input(阻带最小衰减（dB）：);c=input(通带边界频率（Hz）：);d=input(阻带边界频率（Hz）：);T=1/Fs;omegaph=2*c/Fs;omegash=2*d/Fs;N,wp=cheb1ord(omegaph,omegash,rp,rs)Bz,Az=cheby1(N,rp,0.3,0.5)wk=0:pi/512:pi;Hz,w=freqz(Bz,Az,wk);Hx=angle(Hz);subplot(2,1,1);plot(w,Hx);grid on;xlabel(omega/pi);ylabel(相位);title(相频特性);subplot(2,1,2);plot(wk/pi,20*log10(abs(Hz);grid on;xlabel(omega/pi);ylabel(幅度(dB);axis(0,1,-100,5);title(幅频特性); grid on图5切比雪夫2带通滤波器二.FIR滤波器（1）哈明窗低通滤波器Fs=input(采样频率等于（kHz）：);rp=input(通带最大衰减（dB）：);rs=input(阻带最小衰减（dB）：);T=1/Fs;c=input(通带边界频率（kHz）：);d=input(阻带边界频率（kHz）：); wp=2*pi*c/Fs;ws=2*pi*d/Fs;DB=ws-wp;wc=(wp+ws)/2/pi; N=ceil(6.6*pi/DB) hn=fir1(N-1,wc,hamming(N)M=1024;hk=fft(hn,M);n=0:N-1;subplot(3,1,1);stem(n,hn,.);xlabel(n);ylabel(h(n);grid on;k=1:M/2;w=2*(0:M/2-1)/M;subplot(3,1,2);plot(w,20*log10(abs(hk(k);axis(0,1,-80,5);xlabel(w/pi);ylabel(20lg|Hg(w)|);grid onsubplot(3,1,3)phase=unwrap(angle(hk);phase=fftshift(phase);w=0:1023;w=2*w*pi/1024;plot(w/pi,phase/pi);grid on;xlabel(omega/pi);ylabel(相位/pi);hold off;图6 哈明窗低通（2）哈明窗高通Fs=input(采样频率等于（kHz）：);rp=input(通带最大衰减（dB）：);rs=input(阻带最小衰减（dB）：);T=1/Fs;c=input(通带边界频率（kHz）：);d=input(阻带边界频率（kHz）：); wp=2*pi*c/Fs;ws=2*pi*d/Fs;DB=ws-wp;wc=(wp+ws)/2/pi; N0=ceil(6.6*pi/DB); N=N0+mod(N0+1,2) hn=fir1(N-1,wc,high,hamming(N)M=1024;hk=fft(hn,M);n=0:N-1;subplot(3,1,1);stem(n,hn,.);xlabel(n);ylabel(h(n);grid on;k=1:M/2;w=2*(0:M/2-1)/M;subplot(3,1,2);plot(w,20*log10(abs(hk(k);axis(0,1,-80,5);xlabel(w/pi);ylabel(20lg|Hg(w)|);grid onsubplot(3,1,3)phase=unwrap(angle(hk);phase=fftshift(phase);w=0:1023;w=2*w*pi/1024;plot(w/pi,phase/pi);grid on;xlabel(omega/pi);ylabel(相位/pi);图7 哈明窗高通滤波器（3）汉宁窗带通wls=0.45*pi;wlp=0.55*pi;wup=0.7*pi;wus=0.8*pi;DB=wlp-wls;wp=(wls+wlp)/2/pi,(wus+wup)/2/pi;N0=ceil(6.2*pi/DB);N=N0+mod(N0+1,2)wc=(wlp+wls)/2/pi;hn=fir1(N-1,wp,hanning(N)M=1024;hk=fft(hn,M);n=0:N-1;subplot(3,1,1);stem(n,hn,.);xlabel(n);ylabel(h(n);grid on;k=1:M/2;w=2*(0:M/2-1)/M;subplot(3,1,2);plot(w,20*log10(abs(hk(k);axis(0,1,-80,5);xlabel(w/pi);ylabel(20lg|Hg(w)|);grid onsubplot(3,1,3)phase=unwrap(angle(hk);phase=fftshift(phase);w=0:1023;w=2*w*pi/1024;plot(w/pi,phase/pi);grid on;xlabel(omega/pi);ylabel(相位/pi); grid on图8 汉宁窗带通滤波器3. 用滤波器对动态信号进行滤波用自己设计的各滤波器分别对加噪的语音信号进行滤波，在Matlab中，FIR滤波器利用函数fftfilt对信号进行滤波，IIR滤波器利用函数filter对信号进行滤波。3.比较滤波前后语音信号的波形及频谱以下以IIR中的巴特沃斯低通与高通和FI中的哈明窗的低通与高通的输入与输出波形进行示例分析。IIR低通滤波后的波形Fs=input(采样频率等于：);rp=input(通带最大衰减（dB）：);rs=input(阻带最小衰减（dB）：);T=1/Fs;val=get(hObject,Value);str=get(hObject,String);switch strval case 巴特沃斯低通滤波后的波形 c=input(通带边界频率（Hz）：);d=input(阻带边界频率（Hz）：); wp=(2/T)*tan(2*pi*c/Fs)/2);ws=(2/T)*tan(2*pi*d/Fs)/2);N,wc=buttord(wp,ws,rp,rs,s);B,A=butter(N,wc,s);Bz,Az=bilinear(B,A,Fs);wk=0:pi/512:pi;Hz,w=freqz(Bz,Az,wk);Hx=angle(Hz); fs=input(频率等于：);m=1/fs;f0=100;t=-0.2:m:0.2;sum=0;for n=1:1:9;fn=f0+n*10;endfor n=1:1:9;sum=sum+100*sin(2*pi*fn*t);x1=sum;endtitle(信号叠加后动态型信号);plot(t,sum,linewidth,1.5),gridS=fft(x1,512);subplot(311);plot(t,x1); sf=filter(Bz,Az,x1);subplot(312); plot(t,sf); xlabel(T(s);ylabel(T w); title(滤波后信号的时域波形); SF=fft(sf,512); w=(0:255)/256*(Fs/2); subplot(313); plot(w,abs(S(1:256) SF(1:256); xlabel( (Hz);ylabel(Mag.of Fourier transform);grid on;legend(before,after); title(滤波前后信号对比);低通滤波器图9巴特沃斯低通滤波后的波形FIR低通滤波后的波形 case 哈明窗低通滤波后的波形 c=input(通带边界频率（kHz）：);d=input(阻带边界频率（kHz）：); wp=2*pi*c/Fs;ws=2*pi*d/Fs;DB=ws-wp;N=ceil(6.2*pi/DB);wc=(wp+ws)/2/pi;hn=fir1(N-1,wc,hanning(N);M=1024;hk=fft(hn,M);n=0:N-1;k=1:M/2;w=2*(0:M/2-1)/M;phase=unwrap(angle(hk);phase=fftshift(phase);w=0:1023;w=2*w*pi/1024; fs=input(频率等于（Hz）：);m=1/fs;f0=100;t=-0.2:m:0.2;sum=0;for n=1:1:9;fn=f0+n*10;endfor n=1:1:9;sum=sum+100*sin(2*pi*fn*t);x1=sum;endS=fft(x1,512); subplot(311);plot(t,x1); sf=filter(hn,1,x1);subplot(312); plot(t,sf); xlabel(Time(seconds);ylabel(Time waveform); title(滤波后信号的时域波形); SF=fft(sf,512); w=(0:255)/256*(Fs/2);subplot(313); plot(w,abs(S(1:256) SF(1:256); xlabel(Frequency (Hz);ylabel(Mag.of Fourier transform);grid on;legend(before,after); title(滤波前后信号对比);grid on图10 哈明窗低通滤波后的波形IIR高通滤波后的波形：巴特沃斯高通滤波后的波形： case 巴特沃斯高通滤波后的波形 c=input(通带边界频率（Hz）：);d=input(阻带边界频率（Hz）：); wp=(2/T)*tan(2*pi*c/Fs)/2);ws=(2/T)*tan(2*pi*d/Fs)/2);N,wc=buttord(wp,ws,rp,rs,s);B,A=butter(N,wc,high,s);Bz,Az=bilinear(B,A,Fs);wk=0:pi/512:pi;Hz,w=freqz(Bz,Az,wk);Hx=angle(Hz); fs=input(频率等于：);m=1/fs;f0=100;t=-0.2:m:0.2;sum=0;for n=1:1:9;fn=f0+n*10;endfor n=1:1:9;sum=sum+100*sin(2*pi*fn*t);x1=sum;endtitle(信号叠加后动态型信号);plot(t,sum,linewidth,1.5),gridS=fft(x1,512);subplot(311);plot(t,x1); sf=filter(Bz,Az,x1);subplot(312); plot(t,sf); xlabel(T(s);ylabel(T w); title(滤波后信号的时域波形); SF=fft(sf,512); w=(0:255)/256*(Fs/2); subplot(313); plot(w,abs(S(1:256) SF(1:256); xlabel( (Hz);ylabel(Mag.of Fourier transform);grid on;legend(before,after); title(滤波前后信号对比);grid on 图11 巴特沃斯高通滤波后的波形FIR高通滤波后的波形 c=input(通带边界频率（kHz）：);d=input(阻带边界频率（kHz）：); wp=2*pi*c/Fs;ws=2*pi*d/Fs;DB=ws-wp;wc=(wp+ws)/2/pi;N0=ceil(6.6*pi/DB); N=N0+mod(N0+1,2) ; hn=fir1(N-1,wc,high,hamming(N);M=1024;hk=fft(hn,M);n=0:N-1;k=1:M/2;w=2*(0:M/2-1)/M;phase=unwrap(angle(hk);phase=fftshift(phase);w=0:1023;w=2*w*pi/1024; fs=input(频率等于（Hz）：);m=1/fs;f0=100;t=-0.2:m:0.2;sum=0;for n=1:1:9;fn=f0+n*10;endfor n=1:1:9;sum=sum+100*sin(2*pi*fn*t);x1=sum;endS=fft(x1,512); subplot(311);plot(t,x1); sf=filter(hn,1,x1);subplot(312); plot(t,sf); xlabel(Time(seconds);ylabel(Time waveform); title(滤波后信号的时域波形); SF=fft(sf,512); w=(0:255)/256*(Fs/2);su