使用最优等波纹法设计FIR录波器

1、音乐信号滤波去噪使用最优等波纹法设计的FIR滤波器学生姓名：王俊 指导老师：曹敦摘要本课程设计是设计和使用FIR滤波器滤除信号中的噪声。本次设计是在MATLAB软件平台上进行的，采用最优等波纹法设计FIR滤波器。原始信号是录音采集的音乐信号，加入一个单频噪声后当作干扰后的信号，通过设计的FIR滤波器还原成原始信号。本次课程设计中设计的滤波器起到了很好的滤波作用，滤波后的信号和原始信号基本一样，达到了设计目的。关键词音乐信号；FIR滤波器；最优等波纹法；MATLAB1 引言滤波器作为通信系统的基本模块，是指执行信号处理功能的电子系统，它专门用于去除信号中不想要的或者增强所需成分。根据性质，滤波器

2、可以分为非线性的、线性的、时不变的、时变的（自适应的），连续的、离散的（数字的）、无限脉冲响应（IIR）的、有限脉冲响应（FIR）的等。选择和过滤信号是滤波器的重要功能。从频率域上看，就是将有用的信号频率成分选择出来，而阻止其他频率成分的信号或干扰。根据信号过滤的频域特征，又可将滤波器分为低通、带通、高通、带通、带阻、全通以及梳状滤波器等类型。能够通过滤波器的信号频率部分称为通带，而被阻止的频率部分称为阻带。1.1课程设计的目的（1） 巩固和扩大所学数字信号处理，信号与系统及通信原理等专业课程的知识。（2） 综合使用这些知识和技术完成通信系统中的数字信号处理算法的设计，激发学习兴趣和培养独立思

3、考和动手能力。1.2课程设计的要求进一步熟悉数字信号处理算法和数字信号处理系统。（1） 将课程设计内容用MATLAB实现。（2） 设计过程中，要求独立思考，刻苦钻研，按时、全面的完成规定的设计任务。（3） 设计报告必须由设计者独立编写，不得抄袭。报告内容简明扼要，条例清楚。（4） 按时完成设计任务，准时上交设计报告的打印件。1.3设计平台本次课程设计的设计平台是MATLAB。MATLAB是由美国mathworks公司发布的主要面对科学计算、可视化以及交互式程序设计的高科技计算环境。它将数值分析、矩阵计算、科学数据可视化以及非线性动态系统的建模和仿真等诸多强大功能集成在一个易于使用的视窗环境中，

4、为科学研究、工程设计以及必须进行有效数值计算的众多科学领域提供了一种全面的解决方案，并在很大程度上摆脱了传统非交互式程序设计语言（如C、Fortran）的编辑模式，代表了当今国际科学计算软件的先进水平。MATLAB 的应用范围非常广，包括信号和图像处理、通讯、控制系统设计、测试和测量、财务建模和分析以及计算生物学等众多应用领域。附加的工具箱（单独提供的专用 MATLAB 函数集）扩展了 MATLAB 环境，以解决这些应用领域内特定类型的问题。特点1) 高效的数值计算及符号计算功能,能使用户从繁杂的数学运算分析中解脱出来; 2) 具有完备的图形处理功能,实现计算结果和编程的可视化; 3) 友好的

5、用户界面及接近数学表达式的自然化语言,使学者易于学习和掌握; 4) 功能丰富的应用工具箱(如信号处理工具箱、通信工具箱等) ,为用户提供了大量方便实用的处理工具。2 设计原理2.1 FIR滤波器FIR(Finite Impulse Response)滤波器：有限长单位冲激响应滤波器，是数字信号处理系统中最基本的元件，它可以在保证任意幅频特性的同时具有严格的线性相频特性，同时其单位抽样响应是有限长的，因而滤波器是稳定的系统。因此，FIR滤波器在通信、图像处理、模式识别等领域都有着广泛的应用。有限长单位冲激响应（FIR）滤波器有以下特点： (1) 系统的单位冲激响应h (n)在有限个n值处不为零；

6、 (2) 系统函数H(z)在|z|>0处收敛，极点全部在z = 0处（因果系统）； (3) 结构上主要是非递归结构，没有输出到输入的反馈，但有些结构中（例如频率抽样结构）也包含有反馈的递归部分。 设FIR滤波器的单位冲激响应h (n)为一个N点序列，0 n N 1，则滤波器的系统函数为 H(z)=h(n)*z-n（2.1）就是说，它有（N1）阶极点在z = 0处，有（N1）个零点位于有限z平面的任何位置。FIR滤波器的基本结构有：横截型（卷积型、直接型）、级联型、频率抽样型、快速卷积结构。FIR滤波器的基本设计方法有：窗口设计法、频率采样设计法、最优等波纹设计法。2.2最优等波纹设计法窗

7、口法设计和频率采样设计都存在某些缺陷。首先，在设计中不能将边缘频率和精确地给定。其次，不能够同时标定波纹因子1和2。最后，近似误差（即理想响应和实际响应之间的差）在频率区间上不是均匀分布的。而最优等纹波设计法能解决上面三个问题。对于线性相位FIR滤波器来说，有可能导得一组条件，对这组条件能够证明，在最大近似误差最小化的意义下这个设计解是最优的（有时就称最大值最小或切比雪夫（Chebyshev）误差）。具有这种性质的滤波器就是通过最优等波纹设计法设计得到的等波纹滤波器。切比雪夫近似问题现在能定义为：确定这组系数或或或或等效为a(n)或b(n)或c(n)或d(n)以使在通带和阻带内E(w)的最大绝

8、对值最小，即（2.2）交错点定理：设S是闭区间0,内任意闭合子集，为使P（）是在S上对的唯一最大值最小近似，其必要与充分条件是E（）在S内至少呈现出（L+2）个“交错点”或极值频率；这就是说，在S内一定存在（L+2）个频率使之有E()=-E()=（2.3） 将这个定理与前面的结论结合在一起，表明最优等波纹滤波器在S内它的误差函数不是有（L+2）个就是有（L+3）个交错点。Parks-McClellan算法：假定滤波器长度M（或L）和比值已知，选取加权函数，正确的选定阶M，当这个解得到时就有。M和是互为关联的；M愈大，就愈小。近似的M可由下面这个公式得到。（2.4）Parks-McClellan

9、算法从估算（L+2）个极值频率开始并估计出在这些频率上的最大误差。然后通过（2.3）式给出的点拟合一个L阶多项式。在一个很细的密度上确定局部最大误差，并在这些新的极值上调整极值频率。通过这些新的极值频率又拟合出一个新的L阶多项式，这个过程一直重复下去。这一迭代过程一直持续到最优一组频率和全局最大误差被找到为止。最后脉冲响应被计算出来。由于是对M的近似，最大误差可能不等于。若是这样，那么必须增大M（若>），或者减小M（若<），再求出一个新的。重复这个过程直到<为止。这样，最优等波纹FIR滤波器就确定了。Parks-McClellan算法在MATLAB中作为一个称为firpm的函

10、数是可以得到的，这个函数一种的句法是h=firpm(N,f,m,weights)；设计一N阶（注意，滤波器长度是M=N+1）FIR数字滤波器，其频率响应由数组f和m给定。数组f包含以为单位的频带边缘频率，也即0.0<f<1.0，这些频率必须是是以递增的次序从0.0开始，并以1.0终结。数组m包含有在f所给定频率上的期望幅度响应。数组weights给出了再每个频带内的加权函数。为了估计出滤波器的阶N，SP工具箱提供了函数firpmord。这个函数也能估计在firpm函数中用到的其他参数。这个基本语法是N,f0,m0,weights=firpmord(f,m,delta)；这个函数计算

11、窗的阶N，在f0中的归一化频带，m0中的振幅响应，以及在weights中的频带加权值。向量f是归一化频带边缘频率向量，m是由f定义的频带上给定期望振幅值的向量。F的长度小于两倍m的长度，也即f不包括0或1。向量delta给定在每个频带内的容度。估计出的参数现在能用于firpm函数。3设计步骤3.1 设计流程图图3.1总体设计流程图图3.2最优等波纹法设计FIR滤波器流程图本课程设计的总体设计流程如图3.1所示。利用最优等波纹法设计FIR滤波器的流程图如图3.2所示。3.2 录制语音信号在电脑上自己录制一段语音信号，且为单声道，再将此.格式语音控制在3秒内，以减少设计中的误差。然后在Matlab

12、软件平台下，利用函数wavread对语音信号进行采样，记住采样频率和采样点数，采集完成后在信号中加入一个单频噪声，绘制原语音信号和加噪后的语音信号的时域和频域的波形图。具体调用如下：x,fs,bits=wavread('e:wangjun.wav'); % 输入参数为文件的全路径和文件名，输出的第一个参数是每个样本的值，fs是生成该波形文件时的采样率，bits是波形文件每样本的编码位数。sound(x,fs,bits); % 按指定的采样率和每样本编码位数回放N=length(x); % 计算信号x的长度fn=2100; % 单频噪声频率，此参数可改t=0:1/fs:(N-1)

13、/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); % 对原始信号和加噪信号进行fft变换,取幅度谱X=X(1:N/2) ; Y=Y(1:N/2); % 截取前半部分deltaf=fs/N; % 计算频谱的谱线间隔f=0:deltaf:fs/2-deltaf; % 计算频谱频率范围得到的波形图如图3.2所示：图3-2 原始语音信号与加噪后的语音信号时域、频域波形图3.3 滤波器设计（1） 设计数字滤

14、波器因为噪声的频率是1500Hz，所以设计各滤波器的性能指标如下所示：，As50dB，Rp1 dBfp1=1400;fs1=1450;fs2=1550;fp2=1600; Rp=1;As=50;% 带阻滤波器设计指标wp1=fp1/Fs*2*pi;ws1=fs1/Fs*2*pi;ws2=fs2/Fs*2*pi;wp2=fp2/Fs*2*pi;%将Hz为单位的模拟频率换算为rad为单位的数字频率delta1,delta2=db2delta(Rp,As);%把以dB计的通带波纹和阻带衰减转换成相应的容度（或波纹）f1=wp1,ws1,ws2,wp2/pi;m=1,0,1;delta=delta1,

15、delta2,delta1;%设置firpmord的输入参数N,f1,m,weights=firpmord(f1,m,delta);N%估猜窗的阶NN = 288 h=firpm(N,f1,m,weights);%通过Parks-McClellan算法计算N阶FIR数字滤波器的频率响应db,mag,pha,grd,w=freqz_m(h,1); % 调用自编函数计算滤波器的频率特性 delta_w=2*pi/1000; ws1i=floor(ws1/delta_w)+1;wp1i=floor(wp1/delta_w)+1;ws2i=floor(ws2/delta_w)+1;wp2i=floor

16、(wp2/delta_w)+1; Asd=-max(db(ws1i:1:ws2i)Asd = 43.7257%校核最小阻带衰减并与As比较，然后增大（或降低）这个滤波器系数%- N=N+1;h=firpm(N,f1,m,weights);db,mag,pha,grd,w=freqz_m(h,1); % 调用自编函数计算滤波器的频率特性Asd=-max(db(ws1i:1:ws2i)Asd = 42.3696%校核最小阻带衰减并与As比较，然后增大（或降低）这个滤波器系数%-中间省略了一些重复的步骤N=N+1;h=firpm(N,f1,m,weights);db,mag,pha,grd,w=fr

17、eqz_m(h,1); % 调用自编函数计算滤波器的频率特性Asd=-max(db(ws1i:1:ws2i)Asd = 50.1366%一直比较判断Asd和As的大小，当Asd第一次大于等于As时完成%此时Asd>As，完成 NN = 413 M=N+1M = 414 %滤波器长度M=N+1（2） 画图分析滤波器的各参数 n=0:M; subplot(2,2,1);stem(n,h);axis(0,M,-0.1,1)xlabel('n');ylabel('h(n)');title('滤波器脉冲响应图');grid on subplot(2,

18、2,2);plot(w/pi,db); axis(0,1,-60,5); set(gca,'Xtick',0，(wp1+ws1)/(2*pi),(wp2+ws2)/(2*pi)，1,'Ytick',-50,-1);grid on xlabel('w/pi');ylabel('dB');title('滤波器幅度响应图'); subplot(2,2,3);plot(w/pi,pha); axis(0,1,-4,4);grid onxlabel('w/pi');ylabel('相位');t

19、itle('滤波器相位响应图'); subplot(2,2,4);plot(w/pi,mag); axis(0,1,0,1.1);grid onxlabel('w/pi');ylabel('幅度mag');title('滤波器幅度响应图');图 3.4 最优等波纹的滤波器用最优等波纹法设计的FIR滤波器的各参数特性如图3.4所示。观图可知，FIR滤波器符合设计要求。3.4 信号滤波处理用滤波器对信号进行滤波y_fil,ny=conv_m(y1,t,h,n);%用设计好的滤波器处理干扰后的音乐信号 y_fil=y_fil(1:400

20、00);%限定滤波后信号的长度 Y_fil=abs(fft(y_fil);% 对原始信号和加噪信号进行fft变换,取幅度谱 Y_fil=Y_fil(1:N/2);% 截取前半部分3.5 结果分析(1) 在波形及频谱上比较滤波前后音乐信号subplot(3,2,1);plot(t,y);axis(0,5,-1.3,1.3);grid onxlabel('时间(单位:s)');ylabel('幅度');title('原始音乐信号'); subplot(3,2,2);plot(f,Y);grid onxlabel('频率(单位:Hz)'

21、);ylabel('幅度谱');title('音乐信号幅度谱图'); subplot(3,2,3);plot(t,y1);axis(0,5,-1.3,1.3);grid on xlabel('时间(单位:s)');ylabel('幅度');title('加入单频干扰后的音乐信号'); subplot(3,2,4);plot(f,Y1);axis(0,4000,0,1500);grid on xlabel('频率(单位:Hz)');ylabel('幅度谱');title('加入

22、干扰后的音乐信号幅度谱图'); subplot(3,2,5);plot(t,y_fil);axis(0,5,-1.3,1.3);grid on xlabel('时间(单位:s)');ylabel('幅度');title('滤波后的音乐信号');subplot(3,2,6);plot(f,Y_fil);grid on xlabel('频率(单位:Hz)');ylabel('幅度谱');title('滤波后音乐信号幅度谱图');图.滤波后信号幅度和频谱与原始信号和加噪后信号的比较分析图3.5，通

23、过FIR滤波器处理后的音乐信号和原始音乐信号在时域上基本相同，只是在频域1500Hz左右衰减较严重。(2) 在回放中比较滤波前后音乐信号wavplay(y,Fs);%播放原始音乐信号wavplay(y1,Fs);%播放加入单频噪声的音乐信号wavplay(y_fil,Fs);%回放滤波后的音乐信号，与原始信号进行对比通过在Matlab上回放原始音乐信号、加噪后的音乐信号和滤波后的音乐信号进行比较，滤波后的音乐信号在音质上和原始音乐基本相同，说明设计的FIR滤波器滤波效果较明显，设计成功。4 出现的问题及解决方法1、加入的噪声如果不限定其振幅，对原始信号的影响较大。为了便于观察和分析，限定其振幅

24、为0.3的单频噪声（0.3*sin(fn*2*pi*t)）。2、在开始录制音乐信号并将其导入MATLAB中的时候出现过错误，原因是我录制的语音信号是双音频信号，不符合要求，在老师的指导下我将语音信号变成了单声道信号，再次导入的时候MATLAB不再报错。3、在用最优等波纹法设计FIR滤波器的时候，因为本滤波器的设计指标相对较高，所以最后得到的滤波器的阶数会比较大，如果迭代运算时N的值每次只加1，会很耗时。为了解决这个问题，在通过firpmord函数得到近似的滤波器的阶后，开始可以将N的值每次加大点的数值，到发现Asd的值和As很接近了的时候，再减小N每次加的数值，直到Asd刚大于或等于As。4、

25、在用设计好的FIR滤波器处理加噪信号时，不太会用fftfilt函数，所以就使用conv_m代替。这两个函数的内部算法应该差不多，也达到了同样的效果。5 结束语通过本次课程设计，我对数字信号处理有了更具体的认识，更深地了解了滤波器的设计一般设计步骤，尤其是对利用最优等波纹法设计FIR滤波器。在本次课程设计中，我的题目是音乐信号的滤波去噪，相对平时作业题更贴近现实。不像作业中设计滤波器最多是处理余弦信号等，处理音乐信号让我更形象地认识到滤波器在实际生活中的应用，也增强了学习兴趣。还有就是在写课程设计报告时，你需要对原先课堂上学的知识有更透彻的认识，能对其原理和来龙去脉一清二楚。所以写报告在一定程度

26、上也是整理思路，加深对课程的理解。课程设计是我们运用所学知识，动手实践的一个很好的机会。它既可以帮助我们加深对所学知识的理解，又能提高我们运用知识，联系实际，动手实践的能力。而且在设计过程中可能用到我们没学过的知识，需要我们去查阅资料获取相关信息，这又提高了我们查找信息和学习新知识的能力。在实物的调试与检测过程中，又会遇到许多意想不到的问题，需要我们去分析原因和解决问题。也体会到真正的去独立地完成一件事情是很困难的，同学以及老师的帮助和提醒是必须的。通过这次课程设计，我拓宽了知识面，锻炼了实际操作能力，综合素质也得到了提高，进一步加深了了我们对专业的认识和激发了我们对专业的兴趣。虽然课程设计结

27、束了，但是我们的学习还没结束，对知识的进一步学习还需要继续，很开心成功地完成了这次设计。6 参考文献1美维纳·K·英格尔，约翰·G·普罗克斯.刘叔棠.数字信号处理（MATLAB版）Digital Signal Processing Using MATLAB.第2版.西安：西安交通大学出版社，2008年2 张志涌 等.精通MATLAB 6.5版.北京：北京航空航天大学出版社，2003年3邵玉斌.Matlab/Simulink通信系统建模与仿真实例分析. 北京：清华大学出版社，20084 邓华MATLAB通信仿真及应用实例详解北京：国防工业出版社，2003附

28、录1：音乐信号滤波去噪使用最优等波纹法设计的FIR滤波器 程序一：%放出原音频x,fs,bits=wavread('e:wangjun1.wav'); % 输入参数为文件的全路径和文件名，%输出的第一个参数是每个样本的值，%fs是生成该波形文件时的采样率，%bits是波形文件每样本的编码位数。sound(x,fs,bits); % 按指定的采样率和每样本编码位数回放程序二：%原音乐信号和加噪后的音乐信号的时域和频域的波形图绘制N=length(x); % 计算信号x的长度fn=2100; % 单频噪声频率，此参数可改t=0:1/fs:(N-1)/fs; % 计算时间范围，样本数

29、除以采样频率x=x(:,1)' y=x+0.1*sin(fn*2*pi*t); sound(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; % 计算频谱频率范围subplot(2,2,1);plot(t,x);xlabel('时间(t)');ylabel('幅度'); title(&#

30、39;原始音乐信号');axis(0,2.5,-1.5,1.5);grid on;subplot(2,2,2);plot(f,X);xlabel('频率(f)');ylabel('幅度谱');title('原始音乐信号幅度谱');axis(0，800,0,800);grid on;subplot(2,2,3);plot(t,y);xlabel('时间(t)');ylabel('幅度');title('加干扰后的音乐信号');axis(0,2.5,-1.5,1.5);grid on;subpl

31、ot(2,2,4);plot(f,Y);xlabel('频率(f)');ylabel('幅度谱');title('加干扰后的音乐信号幅度谱');axis(0,800,0,800);grid on;程序三：fp1=1400;fs1=1450;fs2=1550;fp2=1600; Rp=1;As=50;% 带阻滤波器设计指标wp1=fp1/Fs*2*pi;ws1=fs1/Fs*2*pi;ws2=fs2/Fs*2*pi;wp2=fp2/Fs*2*pi;%将Hz为单位的模拟频率换算为rad为单位的数字频率delta1,delta2=db2delta(Rp

32、,As);%把以dB计的通带波纹和阻带衰减转换成相应的容度（或波纹）f1=wp1,ws1,ws2,wp2/pi;m=1,0,1;delta=delta1,delta2,delta1;%设置firpmord的输入参数N,f1,m,weights=firpmord(f1,m,delta);N%估猜窗的阶NN = 288 h=firpm(N,f1,m,weights);%通过Parks-McClellan算法计算N阶FIR数字滤波器的频率响应db,mag,pha,grd,w=freqz_m(h,1); % 调用自编函数计算滤波器的频率特性 delta_w=2*pi/1000; ws1i=floor(

33、ws1/delta_w)+1;wp1i=floor(wp1/delta_w)+1;ws2i=floor(ws2/delta_w)+1;wp2i=floor(wp2/delta_w)+1; Asd=-max(db(ws1i:1:ws2i)Asd = 43.7257%校核最小阻带衰减并与As比较，然后增大（或降低）这个滤波器系数%- N=N+1;h=firpm(N,f1,m,weights);db,mag,pha,grd,w=freqz_m(h,1); % 调用自编函数计算滤波器的频率特性Asd=-max(db(ws1i:1:ws2i)Asd = 42.3696%校核最小阻带衰减并与As比较，然后

34、增大（或降低）这个滤波器系数%-中间省略了一些重复的步骤N=N+1;h=firpm(N,f1,m,weights); db,mag,pha,grd,w=freqz_m(h,1); % 调用自编函数计算滤波器的频率特性Asd=-max(db(ws1i:1:ws2i)Asd = 50.1366%一直比较判断Asd和As的大小，当Asd第一次大于等于As时完成%此时Asd>As，完成NN = 413M=N+1M = 414 %滤波器长度M=N+1（3） 画图分析滤波器的各参数n=0:M;subplot(2,2,1);stem(n,h);axis(0,M,-0.1,1)xlabel('n

35、');ylabel('h(n)');title('滤波器脉冲响应图');grid onsubplot(2,2,2);plot(w/pi,db); axis(0,1,-60,5);set(gca,'Xtick',0，(wp1+ws1)/(2*pi),(wp2+ws2)/(2*pi)，1,'Ytick',-50,-1);grid onxlabel('w/pi');ylabel('dB');title('滤波器幅度响应图');subplot(2,2,3);plot(w/pi,pha

36、); axis(0,1,-4,4);grid onxlabel('w/pi');ylabel('相位');title('滤波器相位响应图');subplot(2,2,4);plot(w/pi,mag); axis(0,1,0,1.1);grid onxlabel('w/pi');ylabel('幅度mag');title('滤波器幅度响应图');程序四：在波形及频谱上比较滤波前后音乐信号subplot(3,2,1);plot(t,y);axis(0,5,-1.3,1.3);grid onxlabel

37、('时间(单位:s)');ylabel('幅度');title('原始音乐信号'); subplot(3,2,2);plot(f,Y);grid onxlabel('频率(单位:Hz)');ylabel('幅度谱');title('音乐信号幅度谱图'); subplot(3,2,3);plot(t,y1);axis(0,5,-1.3,1.3);grid on xlabel('时间(单位:s)');ylabel('幅度');title('加入单频干扰后的音乐信号

38、'); subplot(3,2,4);plot(f,Y1);axis(0,4000,0,1500);grid on xlabel('频率(单位:Hz)');ylabel('幅度谱');title('加入干扰后的音乐信号幅度谱图'); subplot(3,2,5);plot(t,y_fil);axis(0,5,-1.3,1.3);grid on xlabel('时间(单位:s)');ylabel('幅度');title('滤波后的音乐信号');subplot(3,2,6);plot(f,Y_f

39、il);grid on xlabel('频率(单位:Hz)');ylabel('幅度谱');title('滤波后音乐信号幅度谱图');在回放中比较滤波前后音乐信号wavplay(y,Fs);%播放原始音乐信号wavplay(y1,Fs);%播放加入单频噪声的音乐信号wavplay(y_fil,Fs);%回放滤波后的音乐信号，与原始信号进行对比附录（程序用到的两个函数）： ：freqz_m的定义源程序清单% 程序名称：freqz_m的M文件内容% 程序功能：定义Matlab中freqz_m函数语句% 程序修改日期：2013-3-6function db,mag,pha,grd,w = freqz_m(b,a);% freqz 子程序的改进版本% -%