基于频率采样法FIR带通滤波器设计

1、中国*大学电子信息工程学院DSP课程设计报告题 目：基于频率采样法的FIR带通滤波器的设计专业班级：通信工程组 员：*指导老师：*2011年6月13日15目 录一、设计目的3二、设计要求3三、设计原理4四、程序代码6五、调试分析12六、设计总结13七、设计心得13八、参考文献一、 设计目的：1、掌握用频率取样法设计FIR DF的方法，并掌握该方法的MATLAB编程。2、熟悉频率取样理论，熟悉内插函数及其应用。3、了解FIR DF的频率特性和相位特性，观察过渡带取样点对滤波器幅频特性的影响。 二、 设计内容:基于MATLAB结合FFT和IFFT，利用频率采样法设计FIR数字带通滤波器，然后用自己

2、设计的滤波器对采集的加噪后的语音信号进行滤波，并将滤波前后的信号进行比较，回放语音信号。三、设计原理：FIR滤波器的单位脉冲响应是有限长的（），其z变换为的（N-1）阶多项式：可得FIR滤波器的系统差分方程：因此，FIR滤波器又称为卷积滤波器。FIR滤波器的频率响应表达式为： 信号通过FIR滤波器不失真条件是在通带内具有恒定的幅频特性和相位线性特性。理论上可以证明：当FIR滤波器的系数满足下列中心对称条件： 时，滤波器设计在逼近平直幅频特性的同时，还能获得严格的线性相位特性。线性相位FIR滤波器的相位滞后和群延迟在整个频带上是相等且不变的。对于一个N阶的线性相位FIR滤波器，群延迟为常数，即滤

3、波后的信号简单地延迟常数个时间步长。这一特性使通带频域内信号通过滤波器后仍保持原有波形形状而无相位失真。FIR滤波器设计的方法很多，如窗函数法、频率采样法以及其他的各种优化设计方法，本次课程设计使用频率采样法设计FIR带通滤波器。频率采样法是从频域出发。因为有限长序列h(n)又可用其离散傅立叶变换H(k)来唯一表示，H(k)与所要求的FIR滤波器系统函数Hd(z)之间存在着频率取样关系。即Hd(z)在Z平面单位圆上按角度等分的取样值等于Hd(k)的各相应值，就以此Hd(k)值作为实际FIR数字滤波器频率特性的取样值H(k),或者说H(k)正是所要求的频率响应H(ejw)的N各等间隔的取样值。频

4、率采样法就是根据频域采样理论，由滤波特性指标构造希望逼近的滤波器频响函数Hd(ej)，对其在0,2上采样得到:然后，就可求出单位脉冲响应h(n)，或是系统函数H(z)。这样， h(n)或是H(z)就是滤波器的设计结果。频率取样法设计的基本思想：把给出的理想频率响应进行取样，通过IDFT从频谱样点直接求得有限脉冲响应。其设计过程如下频率取样法的关键是正确确定数字频域系统函数H(k)在0，2 内的N 个样点，其约束条件为 H(k)=H(N-k) （m）=- (N-m) 0kN-1频率采样法的优点是可以在频域直接设计，并且适合最优化设计；缺点是采样频率只能等于的整数倍，因而不能确保截止频率的自由取值

5、，要想实现自由地选择截止频率，必须增加采样点数N，但是这又使计算量加大。四、程序代码:语音程序：filename='111's,fs,nbits=wavread(filename);sound(s,fs,nbits); %回放语音信号n = length (s) ; %求出语音信号的长度t=0:1/fs:(n-1)/fs;Y=fft(s); %傅里叶变换figure(1)subplot(2,1,1);plot(s);title('原始信号波形');grid onsubplot(2,1,2);plot(abs(Y);title('原始信号频谱')s

6、ound(s,fs);grid on仿真结果：图一带通滤波器程序：N=40;alfa=(40-1)/2;k=0:N-1;w1=(2*pi/N)*k;T1=0.109021; T2=0.59417456;hrs=zeros(1,5),T1,T2,ones(1,7),T2,T1,zeros(1,9),T1,T2,ones(1,7),T2,T1,zeros(1,4);hdr=0,0,1,1,0,0; wd1=0,0.2,0.35,0.65,0.8,1;k1=0:floor(N-1)/2); k2=floor(N-1)/2)+1:N-1;angH=-alfa*(2*pi)/N*k1,alfa*(2*p

7、i/N*(N-k2);H=hrs.*exp(j*angH);h=real(ifft(H);db,mag,pha,grd,w = freqz_m(h,1);Hr,ww,a,L =Hr_Type2(h);subplot(2,2,1)plot(w1(1:21)/pi,hrs(1:21),'o',wd1,hdr)axis(0,1,-0.1,1.1);title('带通：N=40，T1=0.109021， T2=0.59417456')ylabel('Hr(k)');set(gca,'XTickMode','manual',

8、'XTick',0,0.2,0.35,0.65,0.8,1)set(gca,'YTickMode','manual','YTick',0,0.059,0.109,1);grid %绘制带网格的图像subplot(2,2,2); stem(k,h);axis(-1,N,-0.4,0.4)title('脉冲响应'); ylabel('h(n)'); text(N+1,-0.4,'n')subplot(2,2,3); plot(ww/pi,Hr,w1(1:21)/pi,hrs(1:21),

9、'o');axis(0,1,-0.1,1.1);title('振幅响应')xlabel('频率 （单位：pi）'); ylabel('Hr(w)')set(gca,'XTickMode','manual','XTick',0,0.2,0.35,0.65,0.8,1);set(gca,'YTickMode','manual','YTick',0,0.059,0.109,1);gridsubplot(2,2,4); plot(w/pi,db

10、); axis(0,1,-100,10);gridtitle('幅度响应');xlabel('频率 （单位：pi）'); ylabel('分贝')set(gca,'XTickMode','manual','XTick',0,0.2,0.35,0.65,0.8,1)set(gca,'YTickMode','manual','YTick',-60;0);set(gca,'YTickLabelMode','manual',&#

11、39;YTickLabels',60;0);仿真结果：图二总程序：%设计一个通带为0.35到0.65，阻带下截止为0.2，上截止为0.8Clear;N=40;alfa=(40-1)/2;k=0:N-1;w1=(2*pi/N)*k;T1=0.109021; T2=0.59417456;hrs=zeros(1,5),T1,T2,ones(1,7),T2,T1,zeros(1,9),T1,T2,ones(1,7),T2,T1,zeros(1,4);hdr=0,0,1,1,0,0; wd1=0,0.2,0.35,0.65,0.8,1;k1=0:floor(N-1)/2); k2=floor(N

12、-1)/2)+1:N-1;angH=-alfa*(2*pi)/N*k1,alfa*(2*pi/N*(N-k2);H=hrs.*exp(j*angH);h=real(ifft(H);db,mag,pha,grd,w = freqz_m(h,1);Hr,ww,a,L =Hr_Type2(h);figure(1)subplot(2,2,1)plot(w1(1:21)/pi,hrs(1:21),'o',wd1,hdr)axis(0,1,-0.1,1.1);title('带通：N=40，T1=0.109021， T2=0.59417456')ylabel('Hr(

13、k)');set(gca,'XTickMode','manual','XTick',0,0.2,0.35,0.65,0.8,1)set(gca,'YTickMode','manual','YTick',0,0.059,0.109,1);grid %绘制带网格的图像subplot(2,2,2); stem(k,h);axis(-1,N,-0.4,0.4)title('脉冲响应'); ylabel('h(n)'); text(N+1,-0.4,'n'

14、;)subplot(2,2,3); plot(ww/pi,Hr,w1(1:21)/pi,hrs(1:21),'o');axis(0,1,-0.1,1.1);title('振幅响应')xlabel('频率 （单位：pi）'); ylabel('Hr(w)')set(gca,'XTickMode','manual','XTick',0,0.2,0.35,0.65,0.8,1);set(gca,'YTickMode','manual','YTick&

15、#39;,0,0.059,0.109,1);gridsubplot(2,2,4); plot(w/pi,db); axis(0,1,-100,10);gridtitle('幅度响应');xlabel('频率 （单位：pi）'); ylabel('分贝')set(gca,'XTickMode','manual','XTick',0,0.2,0.35,0.65,0.8,1)set(gca,'YTickMode','manual','YTick',-60;0

16、);set(gca,'YTickLabelMode','manual','YTickLabels',60;0);s,fs,nbits=wavread('sj.wav');%信号de 取样频率为44100HZx=s(:,1);sound(x,fs);L=length(x);f=fs*(0:L-1)/L;t=0:1/fs:(L-1)/fs; %将所加噪声信号的点数调整到与原始信号相同d=0.03*abs(max(x)*cos(2*pi*22000*t)' %噪声为500和3300Hz的余弦信号dz=cos(2*pi*11025

17、*t)'xd=x.*dz;xz=xd+d;sound(xz,fs); %播放加噪声后的语音信号X=fft(x); %求信号的频谱XD=fft(xd); %信号调制后的频谱XZ=fft(xz);figure(2)subplot(3,1,1);plot(t,x)title('未加噪的信号'); xlabel('time s');ylabel('幅度 ');subplot(3,1,2);plot(t,xd)title('调制后的信号'); xlabel('time s');ylabel('幅度 '

18、);subplot(3,1,3);plot(t,xz)title('调制加噪后的信号'); xlabel('time n');ylabel('fuzhi n');figure(3)subplot(3,1,1);plot(f,abs(X);title('原始语音信号频谱');xlabel('频率（单位：Hz）');ylabel('幅度');subplot(3,1,2);plot(f,abs(XD);title('调制后的信号频谱');xlabel('频率（单位：Hz）'

19、);ylabel('幅度');subplot(3,1,3);plot(f,abs(XZ);title('加噪后的信号频谱');xlabel('频率（单位：Hz）');ylabel('幅度');y = fftfilt(h,xd);Y=fft(y);sound(3*y,fs);figure(4)subplot(2,1,1);plot(t,3*y)title('滤波后的信号'); xlabel('time s');ylabel('幅度');subplot(2,1,2);plot(f,abs

20、(Y);title('滤波后的信号频谱'); xlabel('频率（单位：Hz）');ylabel('幅度');仿真结果：图三图四图五图一为原始语音的时域和频域图：语音信号采样频率为44100HZ；(注：为了让信号不失真，采样频率大于信号频率的2倍。)图二为带通滤波器图：通带为0.35到0.65，阻带下截止为0.2，上截止为0.8；图三、图四分别为未加噪声和调制加入噪声后的语音信号时域图和频域图：d=0.03*abs(max(x)*cos(2*pi*22000*t)' 噪声为22000Hz的余弦信号图五为经过滤波器滤除噪声之后的信号时域和

21、频域波形。五、调试分析1、语音部分：我们先利用电脑自带的录音筒录了一段语音，语音的采样频率是44100HZ，然后用MATLAB编写程序进行测试，对其进行时域和频域的分析，得到图一。2、带通滤波器：我们所要设计的基于MATLAB的频率采样法的FIR带通滤波器的参数是通带为0.35到0.65，阻带下截止为0.2，上截止为0.8。在正式设计滤波器之前，我们组员认真、仔细地查阅了很多资料，并参考了滤波器设计等应用图书，最终经过不断调试、修改得到了符合要求的仿真程序，在这个过程中，出现了不少困难和问题：首先，在对话音信号进行时域和频谱分析时没有对话音信号取单通道，由于矩阵大小不一致致使在运算中始终得不到

22、正确结果，调试无法通过。 起初对滤波器的通带和阻带的上下截止频率的设置和加噪后语音信号的中心频率出现了偏差，不能很完美地滤掉噪声。此外，由于滤波器是根据相对频率设计的，再根据=*转换到模拟频率，当语音信号的采样频率变化时，滤波器的通带相应的发生改变，而我们在更换语音信号时忘记了这一点，以至于对有些语音信号不能达到消噪的目的而苦恼不解。当然，在整个过程中，由于不细心出现了各种小问题，最终在队友的齐心努力下一一纠正。3、加噪声过程：根据语音信号的中心频率和频谱分布情况，将话音信号调制到中心频率为11025HZ的载波上，再加入高频噪声，噪声为22000Hz的余弦信号。原先所加的噪声频率在语音信号的的频带内，不利于带通滤波器的滤波处理，这样就会导致滤波器