采用频率采样法的FIR滤波器

1、 焦阳 吹管乐滤波去噪-使用频率采样法的FIR滤波器 第页 共20页 吹管乐滤波去噪基于频率采样法的FIR滤波器学生姓名：焦阳 指导老师：胡双红摘 要 本课程设计主要内容是设计利用频率采样法设计一个FIR滤波器，对一段吹管乐进行滤波去噪处理并根据滤波前后的波形和频谱分析滤波性能。本课程设计仿真平台为MATLAB7.0，开发工具是M语言编程。首先在网上找到一段笛子独奏，加入一单频噪声，对信号进行频谱分析以确定所加噪声频率，设计滤波器进行滤波去噪处理，比较滤波前后的波形和频谱并进行分析。由分析结果可知，滤波器后的音频信号与原始信号基本一致，即设计的FIR滤波器能够去除信号中所加单频噪声，达到了设计

2、目的。关键词 滤波去噪；FIR滤波器；频率采样法；MATLAB 1 引 言滤波去噪1是信号处理中一种非常基本但十分重要的技术。利用滤波可以从复杂的信号中提取所需的信号，一直不需要的信号。滤波器就是这样一种可以在时域和频域对信号进行滤波处理的系统。通常情况下，有用信号和干扰信号是在不同频段上的，于是通过对滤波器的频率特性精心设计就能达到滤波的目的。本课程设计是采用频率采样法设计频率抽样型滤波器，从而对吹管乐信号滤波去噪。通过对比滤波前后的波形图及回放滤波前后的吹管乐信号，来判断滤波器对噪声信号确实有滤除作用。1.1 课程设计目的（1）熟悉使用MATLAB；（2）了解FIR滤波器原理及结构；（3）

3、利用所学数字信号处理想干知识用MATLAB设计一个FIR滤波器；（4）提高自己动手能力；（5）对加噪声的语音信号进行滤波去噪处理，比较滤波前后的时域波形和频谱并进行分析； 1.2 课程设计要求（1）滤波器指标必须符合工程设计；（2）设计完后应检查其频率响应曲线是否满足指标；（3）处理结果和分析结论应该一致，而且应符合理论；（4）独立完成课程设计并按要求编写课程设计报告；1.3 设计平台本课程设计仿真平台为MATLAB7.0。MATLAB的名称源自Matrix Laboratory，1984年由美工Mathworks公司推向市场。它是一种科学计算软件，专门以矩阵的形式处理数据。MATLAB将高性

4、能的数值计算和可视化集成在一起，并提供了大量的内置函数，从而被广泛地应用于科学计算、控制系统、信号处理等领域的分许、仿真和设计工作。1993年MathWorks公司从加拿大滑铁卢大学购得MAPLE软件的使用权，从而以MAPLE为“引擎”开发了符号数学工具箱（Symbolic Math Toolbox）2。2 设计原理用网上找一段吹管乐，绘制波形并且观察其频谱，给定相应技术指标，用频率采样法设计的一个满足指标的频率采样型FIR滤波器，对该信号进行滤波去噪处理，比较滤波前后的波形和频谱进行分析。2.1 FIR滤波器的设计FIR（Finite Implse Response）3滤波器：有限长单位冲激

5、响应滤波器，又称为非递归型滤波器，是数字信号处理系统中最基本的元件，他可以在保证任意幅频特性的同时具有严格的线性相频特性，同时其单位抽样响应是有限长的，因而滤波器是稳定的系统。因此，FIR滤波器在通信、图像处理、模式识别等领域都有着广泛的应用。有限长单位冲激响应（FIR）滤波器有以下特点：（1） 系统的单位冲激响应h（n）在有限个n值处不为0；（2） 系统函数H（z）在|z|>0处收敛，极点全部在z=0处（因果系统）；（3） 结构上主要是非递归结构，没有输出到输入的反馈，但有些结构中（例如频率抽样结构）也包含有反馈的递归部分。2.2 频率采样型结构把一个有限长序列（长度为N点）的z变换H

6、（z）在单位圆上作N等分抽样，就得到H（k），其主值序列就等于h(n)的离散傅里叶变换H（k）。那里也说到用H（k）表示的H（z）的内插公式为 （2-1）这个公式就为FIR滤波器提供了另外一种结构，这种结构由两部分级联组成。 (2-2)其中级联的第一部分为梳状滤波器，其结构如图所示： （2-3） X(n) 1 r(n) -Z-N图2-1 梳状滤波器结构图第二部分由N各谐振器组成的谐振柜。它是由N个一阶网络并联而成，而这每一个一阶网络都是一个谐振器 （2-4） 其结构如下图所示： H（k） Hk(z)W-A 图2-2 一阶谐振器频率抽样型结构特点：（1）它的系数H（k）直接就是在滤波器在处得频率

7、响应。因此，控制 得频率响应是很直接得。（2）结构有两个主要缺点：a所有的相乘系数及H（k）都是复数，应将它们先化成二阶实数，这样乘起来较复杂，增加乘法次数，存储量。b所有谐振器的极点都是在单位圆上，由决定考虑到系数量化的影响，当系数量化时，极点会移动，有些极点就不能被梳状滤波器的零点所抵消。（零点由延时单元决定，不受量化的影响）系统就不稳定了。（3）将一阶网络合并为二阶网络：a第k和第N-k个谐振器合并为一个实系数的二阶网络，因为h(n)是实数，他的DFT也是圆周共轭对称的。 （2-5）因此，可以将第k和第N-k个谐振器合并为一个二阶网络。 （2-6）其中： b第k和第N-k个谐振器合并为一

8、个二阶网络的极点在单位圆内，而不是在单位圆上，因而从频率响应的几何解释可知，它相当于一个有限Q的谐振器。其谐振频率为 。图2-3 二阶网络结构图除了共轭复根外，还有实根。当N=偶数时，有一对实根，它们分别是k=0，k=N/2两个点。 和 （2-7）当N=奇数时，只有一个实根z=r(k=0)，即只有H0(z)。c修正频率抽样结构流图（N=偶数）图2-4 修正频率抽样结构流图（N=偶数） （2-8）修正频率抽样结构流图（N=奇数）图2-5 修正频率抽样结构流图（N=奇数） （2-9）2.3 频率采样法2.3.1 设计思路：这种设计方法是从频域进行设计的一种方法，首先给定一个希望逼近的频率响应。 （

9、2-10）知道H（k）后，由IDFT定义，可以用这N个采样值H（k）来唯一确定有限长序列h(n)，即 （2-11） （2-12） （2-13）内插公式： （2-15）四种线性相位的FIR滤波器如下表2-1所示：表2-1 四种线性相位的FIR滤波器2.3.2 逼近误差及其改进措施：（1）采样点上滤波器的实际响应是严格地和理想频率响应数值相等的；（2）在采样点之间的频率响应则是由各采样点的加权内插函数的延伸叠加而成的，因而有一定的逼近误差，误差大小取决于理想频率响应曲线形状。（3）理想频率响应特性变化越平缓，则内插值越接近理想值，逼近误差越小；（4）如果采样点之间的理想频率特性变化越陡，则内插值与

10、理想值的误差就越大，因而在理想频率特性的不连续点附近，就会产生尖峰和起伏。2.3.3 滤波器性能的改善：（1） 增加过滤带采样点，它可以大大减少震荡，阻带衰减也可以得到进一步改善。一般一点到二点的过滤带采样即可得到满意的结果。（2） 增加采样点密度，过渡带的宽度与采样点数N成反比。但N值意味着或长度的增加，滤波器运算量必然增大4。3 设计步骤3.1 设计流程图开始在网上找到一段笛子独奏，并用wavread函数采集，加入3200Hz的噪声。进行快速傅里叶变换，绘制加噪前后的频谱对比图，和时域对比图。设定滤波器的性能指标用频率采样法设计频率采样型的FIR滤波器滤波器是否符合标准N Y用设计的滤波器

11、进行滤波处理，回放滤波后的信号，并保存比较滤波前后音频信号的波形和频谱回放音乐结束图3-1 流程图3.2 音频信号获取在网上下载一段吹管乐，并截取适当片段。把音频放到电脑录音机里面，把音频属性设置为8000Hz。如图3-2所示：图3-2 音频信号设置然后在MATLAB软件平台下，利用函数wavread对音频信号进行采样，源程序为：x,fs,bits=wavread(梅花三弄.wav)，记住采样频率和采样点数，MATLAB实现得：fs=8000;bits=8。3.3 音频信号的频谱分析 x,fs,bits=wavread('e:梅花三弄.wav'); % 输入参数为文件的全路径和

12、文件名，输出的第一个参数是每个样本的值，fs是生成该波形文件时的采样率，bits是波形文件每样本的编码位数。sound(x,fs,bits); % 按指定的采样率和每样本编码位数回放n=length(x); % 计算信号x的长度fn=3200; % 单频噪声频率，此参数可改t=0:1/fs:(n-1)/fs; % 计算时间范围，样本数除以采样频率x=x(:,1)' % 将双声道转为单声道y=x+0.1*sin(fn*2*pi*t); % 加噪声sound(y,fs,bits); % 应该可以明显听出有尖锐的单频啸叫声X=abs(fft(x); Y=abs(fft(y); % 对原始信号

13、和加噪信号进行fft变换,取幅度谱X=X(1:n/2); Y=Y(1:n/2); % 截取前半部分deltaf=fs/n; % 计算频谱的谱线间隔f=0:deltaf:fs/2-deltaf; % 计算频谱频率范围运行结果如下图：图3-3 加噪前后时域频域对比图由图3-3可以看出，在频域为3200Hz处加入一个单频噪声，而加入噪声之后，时域的波形图出现了明显失真，通过听取原音频信号x和加噪声音频信号y，可以明显听到y音频信号中有一明显尖锐噪声。3.4 滤波器设计设计一个低通滤波器，将单频信号滤出去，源程序如下所示：fcd=2000; fcu=2500;fs=8000; % 定义设计指标，通带截

14、止频率fcd和阻带截止频率fcuRp=0.5;As=40;wp=fcd/fs*2*pi;ws=fcu/fs*2*pi; % 将模拟指标转换为数字指标delta_w=2*pi/1000;Rp=-(min(db(1:1:wp/delta_w+1);As=-round(max(db(ws/delta_w+1:1:501);T1=0.6025,T2=0.127;T=0.38; % 设置自由样本的最优值TdeltaB=wcu-wcd; % 计算过渡带宽deltaBm=1; % 设置自由样本数mN=(m+1)*2*pi/deltaB+1; % 根据过渡带宽和自由样本数计算采样点数NN=ceil(N+mod

15、(N+1,2); % 将N处理为奇数Ns=fix(fcu-fcd)*2*pi/(2*pi*N); % 计算阻带内样本数Np=ceil(N-Ns); % 计算通带内样本数Ak=ones(1,Np),zeros(1,2*Ns+1),ones(1,Np-1); % 得到理想滤波器幅度特性采样值Ak(Np+1)=T;Ak(2*N-Np+1)=T; % 设置过渡带样本为最优值M=length(Ak); alpha=(M-1)/2;k1=0:floor(M-1)/2);k2=floor(M-1)/2)+1:M-1; angH=-alpha*(2*pi)/M*k1,alpha*(2*pi)/M*(M-k2)

16、; % 计算分段线性相位Hk=Ak.*exp(j*angH); % 幅度和相位相乘得实际滤波器频率特性Hkhn=real(ifft(Hk); % 对HK作IFFT得到脉冲响应hndb,mag,pha,grd,w=freqz_m(hn,1); % 计算hn对应系统的频率响应% 作图检查幅度特性是否满足预定指标程序运行结果如下图3-4所示：图3-4 滤波器参数图由图3-4可以看出，滤波器的衰减大于设定值As=40，满足性能指标，滤波器的衰减可由增加过度带宽来得到。3.5信号滤波处理源程序如下所示：y_fil=fftfilt(hn,y);Y_fil=fft(y_fil);Y_fil=abs(Y_fi

17、l);Y_fil=Y_fil(1:n/2);程序运行结果如下图3-5所示：图3-5 滤波前后时域、频域对比图频谱取前一半由上图3-5可以看出，在滤波之后时域图能得到恢复，频域图中的单频噪声信号也得到滤除，说明了设计的滤波器能滤除加入的噪声信号，因此说达到了设计要求。3.6结果分析开始通过分析决定设计一个带阻的滤波器来滤除加入的单频噪声，根据噪声的频率来设计阻带的范围。在收集音频信号后，按照步骤用频率采样法设计频率采样型滤波器。由图3-4可知，设计的滤波器达到要求。我们观察到图3-5滤波前后音频型号的波形对比图，发现时域波形中加干扰噪声后有明显变化，不过经过滤波后几乎没有变化，说明设计的滤波器达

18、到要求。再通过听取原始语音信号，加噪信号y和滤波之后的信号y_fil。对比之后，发现滤波器确实滤除了噪声。从理想的角度考虑，该带阻滤波器的阻带带宽应该可以变得更窄，让滤波效果更好，但是这样的采样值会变得非常大。考虑到实际的情况，权衡之后，决定牺牲带宽来使得滤波器的阶数降低，因此在上图3-5中我们可以看出，在噪声频谱左右两边的信号也被滤除了。4 出现的问题及解决办法在这次课程设计当中，由于基础不扎实，出现了很多问题，即MATLAB软件操作不当，也有知识掌握程度不够出现的各种问题。1. 在一开始找音频后，没有修改参数值，导致频率抽样过高，在老师提醒下，改为8000Hz;2. 在调用音乐文件时，没有

19、将文件放在MATLAB的工作文件夹下面，到时文件找不到。还有有程序中用到的各种函数没有放入正确的位置;3. 在绘制加噪前后频率对比图时，留白过多，对比不明显，通过axis函数对横纵坐标进行限定；4. 在设计滤波器的过程中，没有准确理解每一条指令代表的含义，导致程序前后不对应，出现更多错误，通过MATLAB中的错误提示，准确找到错误的那一行代码，进行修改；5. 在设计滤波器的过程中，滤波器的衰减小于开始所设置的值，通过牺牲过渡带和调节过渡带的采样值，即T1和T2来使得衰减大于所设定值AS；6. 最后听取滤波后声音过程当中，没有将其保存，请教同学之后，学会了如何保存滤波后的声音文件。5 结束语这是

20、第二次课程设计，在前面的课程设计当中我们学习到如何使用MATLAB，所以对于MATLAB软件的使用并没有那么陌生了，尽管如此，在使用MATLAB的过程中还是出现了很多错误，比如说忘记添加函数文件，参数前后不对应之类的低级错误。此次课程设计，让我更深入的了解到频率抽样法以及频率抽样型的滤波器，由开始的无从下手，再翻阅书上的例子。首先决定使用低通的滤波器将噪声滤除，但由于是加入的是一个单频的信号噪声，就决定使用带阻滤波器的例子，所以就通过书籍信号与系统上高通的例题来进行修改，当然这过程不是一蹴而就的。经历过一次又一次的错误，猜得出来最后的模型。还有这也是我和其他同学一起讨论来的。这次课设让我对滤波

21、器的类型有了一个更加完整了解，在设计中也使我对一些概念有了更深刻的认识。比如：在滤波器分类方面，我深刻的了解了低通，高通滤波器与带通，带阻滤波器的特性区别。还有在课程设计中每一次的数据输入都有其重要意义，用MATLAB编译程序时，可以根据滤波器指标的要求实时知道对滤波器的影响。通过一次次的调试和权衡使滤波器的性能达到最佳。课程设计不仅要求对滤波器理论的研究，更重要的是培养一种遇到问题解决问题的思维。因为有了这次课程设计，我懂得了书本知识只是实际应用的理论指导。如果仅仅只学习书本知识，不去在实践中运用，那只是停留在知识表面，不知其因的层面。比如在数学计算上，可以将噪声完全滤除。而在这次课设中，若

22、要完全滤除噪音，滤波器的阶数就会增高，在现实生活中是很难实现的，所以噪声是不能完全滤除的。课程设计结束了，我相信这次课程设计对今后的学习是很有帮助的，它让我将理论更好的和实践相结合，提高了动手的能力，也填补了自己学习上的一些不足。这次课设的成功，不仅仅是我一个人的努力的结果，更离不开指导老师与同学的帮助，在此向老师和同学们表示衷心的感谢。参考文献1 吴镇扬，数字信号处理M，高等教育出版社，20042 张圣勤，MATLAB7.0实用教程M，北京：机械工程出版社，20063 程佩青，数字信号处理教程M，北京：清华大学出版社，20024 高西全，丁玉美，数字信号处理M，第三版，西安：西安科大出版社，

23、1994附录：源程序x,fs,bits=wavread('e:梅花三弄.wav'); % 输入参数为文件的全路径和文件名，输出的第一个参数是每个样本的值，fs是生成该波形文件时的采样率，bits是波形文件每样本的编码位数。 sound(x,fs,bits); % 按指定的采样率和每样本编码位数回放 n=length(x); % 计算信号x的长度 fn=3200; % 单频噪声频率，此参数可改t=0:1/fs:(n-1)/fs; % 计算时间范围，样本数除以采样频率x=x(:,1)' % 将双声道转为单声道y=x+0.1*sin(fn*2*pi*t); % 加噪声 sou

24、nd(y,fs,bits); % 应该可以明显听出有尖锐的单频啸叫声 X=abs(fft(x); Y=abs(fft(y); % 对原始信号和加噪信号进行fft变换,取幅度谱X=X(1:n/2); Y=Y(1:n/2); % 截取前半部分deltaf=fs/n; % 计算频谱的谱线间隔 f=0:deltaf:fs/2-deltaf; % 计算频谱频率范围figure(1)subplot(2,2,1);plot(t,x);axis(0,3,-1.3,1.3);grid onxlabel('时间（单位:s)');ylabel('幅度');title('原始音

25、乐信号');subplot(2,2,2);plot( f,X(1:n/2);axis(0,5000,0,4000);grid onxlabel('频率（单位:Hz)');ylabel('幅度谱');title('音乐信号幅度谱图');subplot(2,2,3);plot(t,y);axis(0,3,-1.3,1.3);grid onxlabel('时间（单位:s)');ylabel('幅度');title('加入单频干扰的音乐信号');subplot(2,2,4);plot(f,Y);ax

26、is(0,5000,0,4000);grid onxlabel('频率（单位：Hz）');ylabel('幅度谱');title('加入干扰后的音乐信号幅度谱图');fcd=2000; fcu=2500;fs=8000; % 定义设计指标，通带截止频率fcd和阻带截止频率fcuRp=0.5;As=40;wp=fcd/fs*2*pi;ws=fcu/fs*2*pi; % 将模拟指标转换为数字指标delta_w=2*pi/1000;Rp=-(min(db(1:1:wp/delta_w+1);As=-round(max(db(ws/delta_w+1:1

27、:501);T1=0.6025,T2=0.127;T=0.38; % 设置自由样本的最优值TdeltaB=wcu-wcd; % 计算过渡带宽deltaBm=1; % 设置自由样本数mN=(m+1)*2*pi/deltaB+1; % 根据过渡带宽和自由样本数计算采样点数NN=ceil(N+mod(N+1,2); % 将N处理为奇数Ns=fix(fcu-fcd)*2*pi/(2*pi*N); % 计算阻带内样本数Np=ceil(N-Ns); % 计算通带内样本数Ak=ones(1,Np),zeros(1,2*Ns+1),ones(1,Np-1); % 得到理想滤波器幅度特性采样值Ak(Np+1)=

28、T;Ak(2*N-Np+1)=T; % 设置过渡带样本为最优值M=length(Ak); alpha=(M-1)/2;k1=0:floor(M-1)/2);k2=floor(M-1)/2)+1:M-1; angH=-alpha*(2*pi)/M*k1,alpha*(2*pi)/M*(M-k2); % 计算分段线性相位Hk=Ak.*exp(j*angH); % 幅度和相位相乘得实际滤波器频率特性Hkhn=real(ifft(Hk); % 对HK作IFFT得到脉冲响应hndb,mag,pha,grd,w=freqz_m(hn,1); % 计算hn对应系统的频率响应% 作图检查幅度特性是否满足预定指

29、标figure(2)subplot(2,2,1);plot(w/pi,db)grid on;xlabel('w/pi');ylabel('db');title('滤波器db');X_l=0,0,wp/pi,ws/pi;1,1,wp/pi,ws/pi;Y_l=-As,-Rp,-100,-100;-As,-Rp,0,0; % 在wp，ws，Rp，As处画线以更直观判断设计是否达标，每列参数是每个线条的端点坐标line(X_l,Y_l,'Color','r','LineWidth',2

30、,'LineStyle','-')    % 添加线条，红色，线宽为2axis(0 1 -80 5)subplot(2,2,2);plot(w/pi,mag);grid on;xlabel('w/pi');ylabel('mag');title('滤波器幅度响应');axis(0 1 0 1.2)subplot(2,2,3);plot(w/pi,pha);grid on;xlabel('w/pi');ylabel('pha');title(&