数组m包含有在f所给定频率上的期望幅度响应。
数组weights给出了再每个频带内的加权函数。
为了估计出滤波器的阶N，SP工具箱提供了函数firpmord。
这个函数也能估计在firpm函数中用到的其他参数。
这个基本语法是
[N,f0,m0,weights]=firpmord（f,m,delta）；
这个函数计算窗的阶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软件平台下，利用函数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）/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）设计数字滤波器
因为噪声的频率是1500Hz，所以设计各滤波器的性能指标如下所示：
，As＝50dB，Rp＝1dB
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,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%估猜窗的阶N
N=
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（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]=freqz_m（h,[1]）;%调用自编函数计算滤波器的频率特性
Asd=-max（db（ws1i:
1:
ws2i））
Asd=
50.1366
%一直比较判断Asd和As的大小，当Asd第一次大于等于As时完成
%此时Asd>As，完成
N
N=
413
M=N+1
M=
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（'滤波器脉冲响应图'）;gridon
subplot（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]）;gridon
xlabel（'w/pi'）;ylabel（'dB'）;title（'滤波器幅度响应图'）;
subplot（2,2,3）;plot（w/pi,pha）;axis（[0,1,-4,4]）;gridon
xlabel（'w/pi'）;ylabel（'相位'）;title（'滤波器相位响应图'）;
subplot（2,2,4）;plot（w/pi,mag）;axis（[0,1,0,1.1]）;gridon
xlabel（'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:
40000）;%限定滤波后信号的长度
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]）;gridon
xlabel（'时间（单位:
s）'）;ylabel（'幅度'）;title（'原始音乐信号'）;
subplot（3,2,2）;plot（f,Y）;gridon
xlabel（'频率（单位:
Hz）'）;ylabel（'幅度谱'）;title（'音乐信号幅度谱图'）;
subplot（3,2,3）;plot（t,y1）;axis（[0,5,-1.3,1.3]）;gridon
xlabel（'时间（单位:
s）'）;ylabel（'幅度'）;title（'加入单频干扰后的音乐信号'）;
subplot（3,2,4）;plot（f,Y1）;axis（[0,4000,0,1500]）;gridon
xlabel（'频率（单位:
Hz）'）;ylabel（'幅度谱'）;title（'加入干扰后的音乐信号幅度谱图'）;
subplot（3,2,5）;plot（t,y_fil）;axis（[0,5,-1.3,1.3]）;gridon
xlabel（'时间（单位:
s）'）;ylabel（'幅度'）;title（'滤波后的音乐信号'）;
subplot（3,2,6）;plot（f,Y_fil）;gridon
xlabel（'频率（单位:
Hz）'）;ylabel（'幅度谱'）;title（'滤波后音乐信号幅度谱图'）;
图３.５滤波后信号幅度和频谱与原始信号和加噪后信号的比较
分析图3.5，通过FIR滤波器处理后的音乐信号和原始音乐信号在时域上基本相同，只是在频域1500Hz左右衰减较严重。
（2）在回放中比较滤波前后音乐信号
wavplay（y,Fs）;%播放原始音乐信号
wavplay（y1,Fs）;%播放加入单频噪声的音乐信号
wavplay（y_fil,Fs）;%回放滤波后的音乐信号，与原始信号进行对比
通过在Matlab上回放原始音乐信号、加噪后的音乐信号和滤波后的音乐信号进行比较，滤波后的音乐信号在音质上和原始音乐基本相同，说明设计的FIR滤波器滤波效果较明显，设计成功。
4出现的问题及解决方法
1、加入的噪声如果不限定其振幅，对原始信号的影响较大。
为了便于观察和分析，限定其振幅为0.3的单频噪声（0.3*sin（fn*2*pi*t））。
2、在开始录制音乐信号并将其导入MATLAB中的时候出现过错误，原因是我录制的语音信号是双音频信号，不符合要求，在老师的指导下我将语音信号变成了单声道信号，再次导入的时候MATLAB不再报错。
3、在用最优等波纹法设计FIR滤波器的时候，因为本滤波器的设计指标相对较高，所以最后得到的滤波器的阶数会比较大，如果迭代运算时N的值每次只加1，会很耗时。
为了解决这个问题，在通过firpmord函数得到近似的滤波器的阶后，开始可以将N的值每次加大点的数值，到发现Asd的值和As很接近了的时候，再减小N每次加的数值，直到Asd刚大于或等于As。
4、在用设计好的FIR滤波器处理加噪信号时，不太会用fftfilt函数，所以就使用conv_m代替。
这两个函数的内部算法应该差不多，也达到了同样的效果。
5结束语
通过本次课程设计，我对数字信号处理有了更具体的认识，更深地了解了滤波器的设计一般设计步骤，尤其是对利用最优等波纹法设计FIR滤波器。
在本次课程设计中，我的题目是音乐信号的滤波去噪，相对平时作业题更贴近现实。
不像作业中设计滤波器最多是处理余弦信号等，处理音乐信号让我更形象地认识到滤波器在实际生活中的应用，也增强了学习兴趣。
还有就是在写课程设计报告时，你需要对原先课堂上学的知识有更透彻的认识，能对其原理和来龙去脉一清二楚。
所以写报告在一定程度上也是整理思路，加深对课程的理解。
课程设计是我们运用所学知识，动手实践的一个很好的机会。
它既可以帮助我们加深对所学知识的理解，又能提高我们运用知识，联系实际，动手实践的能力。
而且在设计过程中可能用到我们没学过的知识，需要我们去查阅资料获取相关信息，这又提高了我们查找信息和学习新知识的能力。
在实物的调试与检测过程中，又会遇到许多意想不到的问题，需要我们去分析原因和解决问题。
也体会到真正的去独立地完成一件事情是很困难的，同学以及老师的帮助和提醒是必须的。
通过这次课程设计，我拓宽了知识面，锻炼了实际操作能力，综合素质也得到了提高，进一步加深了了我们对专业的认识和激发了我们对专业的兴趣。
虽然课程设计结束了，但是我们的学习还没结束，对知识的进一步学习还需要继续，很开心成功地完成了这次设计。
6参考文献
[1][美]维纳·K·英格尔，约翰·G·普罗克斯.刘叔棠.数字信号处理（MATLAB版）DigitalSignalProcessingUsingMATLAB.第2版.西安：
西安交通大学出版社，2008年
[2]张志涌等.精通MATLAB6.5版.北京：
北京航空航天大学出版社，2003年
[3]邵玉斌.Matlab/Simulink通信系统建模与仿真实例分析.北京：
清华大学出版社，2008
[4]邓华．MATLAB通信仿真及应用实例详解．北京：
国防工业出版社，2003
附录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;%计算时间范围，样本数除以采样频率
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（'原始音乐信号'）;
axis（[0,2.5,-1.5,1.5]）;gridon;
subplot（2,2,2）;plot（f,X）;
xlabel（'频率（f）'）;ylabel（'幅度谱'）;title（'原始音乐信号幅度谱'）;
axis（[0，800,0,800]）;gridon;
subplot（2,2,3）;plot（t,y）;
xlabel（'时间（t）'）;ylabel（'幅度'）;title（'加干扰后的音乐信号'）;
axis（[0,2.5,-1.5,1.5]）;gridon;
subplot（2,2,4）;plot（f,Y）;
xlabel（'频率（f）'）;ylabel（'幅度谱'）;title（'加干扰后的音乐信号幅度谱'）;
axis（[0,800,0,800]）;gridon;
程序三：
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,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%估猜窗的阶N
N=
288
h=firpm（N,f1,m,weights）;%通过Parks-McClellan算法计算N阶FIR数字滤波器的频率响应
[db,mag,pha,grd,w]=freqz_m（h,[1]）;%调用自编函数计算滤
升级会员 