音乐信号滤波去噪——使用最优等波纹法设计的FIR滤波器

音乐信号滤波去噪——使用最优等波纹

法设计的FIR滤波器

学生姓名：*** 指导老师：***

摘要本课程设计是设计和使用FIR滤波器滤除信号中的噪声。本次设计是在MATLAB软件平台上进行的，采用最优等波纹法设计FIR滤波器。原始信号是录音采集的音乐信号，加入一个单频噪声后当作干扰后的信号，通过设计的FIR滤波器还原成原始信号。本次课程设计中设计的滤波器起到了很好的滤波作用，滤波后的信号和原始信号基本一样，达到了设计目的。

关键词音乐信号；FIR滤波器；最优等波纹法；MATLAB

1 引言

滤波器作为通信系统的基本模块，是指执行信号处理功能的电子系统，它专门用于去除信号中不想要的或者增强所需成分。根据性质，滤波器可以分为非线性的、线性的、时不变的、时变的（自适应的），连续的、离散的（数字的）、无限脉冲响应（IIR）的、有限脉冲响应（FIR）的等。

选择和过滤信号是滤波器的重要功能。从频率域上看，就是将有用的信号频率成分选择出来，而阻止其他频率成分的信号或干扰。根据信号过滤的频域特征，又可将滤波器分为低通、带通、高通、带通、带阻、全通以及梳状滤波器等类型。能够通过滤波器的信号频率部分称为通带，而被阻止的频率部分称为阻带。

1.1课程设计的目的

（1）巩固和扩大所学数字信号处理，信号与系统及通信原理等专业课程的知识。

（2）综合使用这些知识和技术完成通信系统中的数字信号处理算法的设计，激发学习兴趣和培养独立思考和动手能力。

1.2课程设计的要求

进一步熟悉数字信号处理算法和数字信号处理系统。

（1）将课程设计内容用MATLAB实现。

（2）设计过程中，要求独立思考，刻苦钻研，按时、全面的完成规定的设计任务。

（3）设计报告必须由设计者独立编写，不得抄袭。报告内容简明扼要，条例清楚。

（4）按时完成设计任务，准时上交设计报告的打印件。

1.3设计平台

本次课程设计的设计平台是MATLAB。

MATLAB是由美国mathworks公司发布的主要面对科学计算、可视化以及交互式程序设计的高科技计算环境。它将数值分析、矩阵计算、科学数据可视化以及非线性动态系统的建模和仿真等诸多强大功能集成在一个易于使用的视窗环境中，为科学研究、工程设计以及必须进行有效数值计算的众多科学领域提供了一种全面的解决方案，并在很大程度上摆脱了传统非交互式程序设计语言（如C、Fortran）的编辑模式，代表了当今国际科学计算软件的先进水平。

MATLAB 的应用范围非常广，包括信号和图像处理、通讯、控制系统设计、测试和测量、财务建模和分析以及计算生物学等众多应用领域。附加的工具箱（单独提供的专用MATLAB 函数集）扩展了MATLAB 环境，以解决这些应用领域内特定类型的问题。

特点

1) 高效的数值计算及符号计算功能,能使用户从繁杂的数学运算分析中解脱出来;

2) 具有完备的图形处理功能,实现计算结果和编程的可视化;

3) 友好的用户界面及接近数学表达式的自然化语言,使学者易于学习和掌握;

4) 功能丰富的应用工具箱(如信号处理工具箱、通信工具箱等) ,为用户提供了大量方

便实用的处理工具。

2 设计原理

2.1 FIR滤波器

FIR(Finite Impulse Response)滤波器：有限长单位冲激响应滤波器，是数字信号处理

系统中最基本的元件，它可以在保证任意幅频特性的同时具有严格的线性相频特性，同时其单位抽样响应是有限长的，因而滤波器是稳定的系统。因此，FIR 滤波器在通信、图像处理、模式识别等领域都有着广泛的应用。

有限长单位冲激响应（FIR ）滤波器有以下特点： (1) 系统的单位冲激响应h (n)在有限个n 值处不为零；

(2) 系统函数H(z)在|z|>0处收敛，极点全部在z = 0处（因果系统）；

(3) 结构上主要是非递归结构，没有输出到输入的反馈，但有些结构中（例如频率抽样结构）也包含有反馈的递归部分。

设FIR 滤波器的单位冲激响应h (n)为一个N 点序列，0 ≤ n ≤ N —1，则滤波器的系统函数为

H(z)=∑h(n)*z^-n

（2.1）

就是说，它有（N —1）阶极点在z = 0处，有（N —1）个零点位于有限z 平面的任何位置。FIR 滤波器的基本结构有：横截型（卷积型、直接型）、级联型、频率抽样型、快速卷积结构。

FIR 滤波器的基本设计方法有：窗口设计法、频率采样设计法、最优等波纹设计法。

2.2最优等波纹设计法

窗口法设计和频率采样设计都存在某些缺陷。首先，在设计中不能将边缘频率p ω和

s ω精确地给定。其次，不能够同时标定波纹因子δ1和δ2。最后，近似误差（即理想响

应和实际响应之间的差）在频率区间上不是均匀分布的。而最优等纹波设计法能解决上面三个问题。

对于线性相位FIR 滤波器来说，有可能导得一组条件，对这组条件能够证明，在最

大近似误差最小化的意义下这个设计解是最优的（有时就称最大值最小或切比雪夫（Chebyshev ）误差）。具有这种性质的滤波器就是通过最优等波纹设计法设计得到的等波纹滤波器。 切比雪夫近似问题现在能定义为：

确定这组系数)(~n a 或)(~n b 或)(~n c 或)(~n d [或等效为a(n)或b(n)或c(n)或d(n)]以使在

通带和阻带内E(w)的最大绝对值最小，即

)](max min[ωωE S

∈

（2..2）

交错点定理：

设S 是闭区间[0,π]内任意闭合子集，为使P （ω）是在S 上对)(ωdr H 的唯一最大值

最小近似，其必要与充分条件是E （ω）在S 内至少呈现出（L+2）个“交错点”或极值频率；这就是说，在S 内一定存在（L+2）个频率 i ω使之有

E(i ω)=-E(1-i ω)=|)(|max ωE S

±

（2.3）

S L ∈<<

将这个定理与前面的结论结合在一起，表明最优等波纹滤波器在S 内它的误差函数

不是有（L+2）个就是有（L+3）个交错点。 Parks-McClellan 算法：

假定滤波器长度M （或L ）和比值12/δδ已知，选取加权函数，正确的选定阶M ，

当这个解得到时就有2δδ=。M 和δ是互为关联的；M 愈大，δ就愈小。近似的M 可由下面

这个公式得到。

p

s M ωωωω

δδ-=?+?--=;1285.213log 20?2110 （2.4）

Parks-McClellan 算法从估算（L+2）个极值频率{i ω}开始并估计出在这些频率上的

最大误差δ。然后通过（2.3）式给出的点拟合一个L 阶多项式。在一个很细的密度上确定局部最大误差，并在这些新的极值上调整极值频率{i ω}。通过这些新的极值频率又拟合出一个新的L 阶多项式，这个过程一直重复下去。这一迭代过程一直持续到最优一组频率{i ω}和全局最大误差δ被找到为止。最后脉冲响应被计算出来。

由于是对M 的近似，最大误差δ可能不等于2δ。若是这样，那么必须增大M （若

δ>2δ），或者减小M （若δ<2δ），再求出一个新的δ。重复这个过程直到δ<2δ为止。这样，最优等波纹FIR 滤波器就确定了。

Parks-McClellan 算法在MATLAB 中作为一个称为firpm 的函数是可以得到的，这个函数一种的句法是 [h]=firpm(N,f,m,weights)；

设计一N 阶（注意，滤波器长度是M=N+1）FIR 数字滤波器，其频率响应由数组f 和m

给定。数组f 包含以π为单位的频带边缘频率，也即0.0

为了估计出滤波器的阶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 总体设计流程图

估猜滤波器的阶N 得到滤波器系数

判断最小阻带衰减

Asd是否大于等于As

M=N+1

完成，最后一次得到的滤波器系数即最终设计值N=N+1或N=N-1

设计指标

YES

ＮＯ

图3.2 最优等波纹法设计FIR滤波器流程图

本课程设计的总体设计流程如图3.1所示。

利用最优等波纹法设计FIR滤波器的流程图如图3.2所示。

3.2 录制语音信号

（1）音乐信号的采集

本课程设计的原始信号是在Matlab软件平台下通过wavrecord函数录制的一段音乐信号（歌曲“伤心太平洋”中的一段）。利用函数wavrecord对语音信号进行采样，记住采样频率和采样时间。通过wavrecord函数的使用，很快理解了采样频率、采样时间等概念。采集完成后在信号中加入一个频率1500Hz的单频噪声，设计的任务即为从含噪信号中滤除单频噪声，还原原始信号。

>> Fs=8000; %设置录音采样率为8000次/s

>> T=5; %录音时间长度为5s

>> y=wavrecord(T*Fs,Fs);

%以double格式记录样值，单声道，输出值y在-1~1之间

%回车之后开始录音，5s之后将自动执行下一句

>>wavplay(y,Fs); %播放刚刚录音的结果

>> N=length(y); % 计算信号y的长度

>> fn=1500; % 单频噪声频率，此参数可改

>> t=0:1/Fs:(N-1)/Fs; % 计算时间范围，样本数除以采样频率

>> y=y';y1=y+0.3*sin(fn*2*pi*t); %加入单频噪声

>>wavplay(y1,Fs); %播放加入单频噪声的音乐信号

（2）原始信号和加噪信号的幅度和频谱分析

首先画出原始信号和加噪信号的时域波形；然后对原始信号和加噪信号进行快速傅里叶变换，得到信号的频谱特性。

>> Y=abs(fft(y));Y1=abs(fft(y1)); % 对原始信号和加噪信号进行fft变换,取幅度谱>> Y=Y(1:N/2);Y1=Y1(1:N/2); % 截取前半部分

>> deltaf=Fs/N; % 计算频谱的谱线间隔

>> f=0:deltaf:Fs/2-deltaf; % 计算频谱频率范围

>> subplot(2,2,1);plot(t,y);axis([0,5,-1.3,1.3]);grid on

>> xlabel('时间(单位:s)');ylabel('幅度');title('原始音乐信号');

>> subplot(2,2,2);plot( f,Y);grid on

>>xlabel('频率(单位:Hz)');ylabel('幅度谱');title('音乐信号幅度谱图');

>> subplot(2,2,3);plot(t,y1);axis([0,5,-1.3,1.3]);grid on

>> xlabel('时间(单位:s)');ylabel('幅度');title('加入单频干扰后的音乐信号');

>> subplot(2,2,4);plot(f,Y1);axis([0,4000,0,1500]);grid on

>> xlabel('频率(单位:Hz)');ylabel('幅度谱');title('加入干扰后的音乐信号幅度谱图');

1

2345-1-0.500.5

1时间(单位:s)

幅度

原始音乐信号

010002000300040000200

400600

800

频率(单位:Hz)

幅度谱音乐信号幅度谱图

1

2345

-1-0.500.5

1时间(单位:s)

幅度

加入单频干扰后的音乐信号0

1000200030004000

0500

1000

1500频率(单位:Hz)

幅度谱

加入干扰后的音乐信号幅度谱图

图3.3 原始信号和加入噪声后信号的幅度及其幅频特性图

原始信号和加噪信号的幅度和频谱如图3.3所示。加入噪声后，原始信号幅度受到

了一定的干扰。加噪后的信号频谱图中也多了一个1500Hz 的冲击频率。

3.3 滤波器设计

（1） 设计数字滤波器

因为噪声的频率是1500Hz ，所以设计各滤波器的性能指标如下所示：

Hz

f Hz f Hz

f Hz f p s s p 16001550145014002211====，A s ＝50dB ，R p ＝1 dB

>>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('滤波器脉冲响应图');grid on

>> 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]);grid on >> xlabel('w/pi');ylabel('dB');title('滤波器幅度响应图');

>> subplot(2,2,3);plot(w/pi,pha); axis([0,1,-4,4]);grid on

>> xlabel('w/pi');ylabel('相位');title('滤波器相位响应图');

>> subplot(2,2,4);plot(w/pi,mag); axis([0,1,0,1.1]);grid on

>>xlabel('w/pi');ylabel('幅度mag');title('滤波器幅度响应图');

0.35630.3937

1

-50

-1

w/pi

d B

滤波器幅度响应图

100

200300400

0.5

1

n

h (n )

滤波器脉冲响应图

0.51

-4-2024w/pi

相位

滤波器相位响应图0

0.51

0.5

1

w/pi

幅度m a g

滤波器幅度响应图

图 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]);grid on >> xlabel('时间(单位:s)');ylabel('幅度');title('原始音乐信号'); >> subplot(3,2,2);plot(f,Y);grid on

>> xlabel('频率(单位:Hz)');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('加入干扰后的音乐信号幅度谱图'); >> 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('滤波后音乐信号幅度谱图');

1

2345-10

1时间(单位:s)

幅度

原始音乐信号

10002000300040000

500

1000

频率(单位:Hz)

幅度谱

音乐信号幅度谱图

1

2345

-10

1

时间(单位:s)幅度

加入单频干扰后的音乐信号

100020003000400005001000

1500频率(单位:Hz)

幅度谱

加入干扰后的音乐信号幅度谱图0

1

2345

-10

1时间(单位:s)

幅度

滤波后的音乐信号0

1000200030004000

500

1000

频率(单位:Hz)

幅度谱

滤波后音乐信号幅度谱图

图３.５ 滤波后信号幅度和频谱与原始信号和加噪后信号的比较

分析图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）在通过plot函数画加噪后信号的频谱图时，因为噪声的冲击值很大，如果不限定纵坐标显示的范围，原始信号的频谱图会显得很小。所以为了更明显地比较原始信号和加噪后信号在频谱图上的不同，适当地限定了纵坐标的显示范围。

（3）在用最优等波纹法设计FIR滤波器的时候，因为本滤波器的设计指标相对较高，所以最后得到的滤波器的阶数会比较大，如果迭代运算时N的值每次只加1，会很耗时。为了解决这个问题，在通过firpmord函数得到近似的滤波器的阶后，开始可以将N的值每次加大点的数值，到发现Asd的值和As很接近了的时候，再减小N每次加的数值，直到Asd刚大于或等于As。

（4）在用设计好的FIR滤波器处理加噪信号时，不太会用fftfilt函数，所以就使用conv_m代替。这两个函数的内部算法应该差不多，也达到了同样的效果。

5结束语

通过本次课程设计，我对数字信号处理有了更具体的认识，更深地了解了滤波器的设计一般设计步骤，尤其是对利用最优等波纹法设计FIR滤波器。

在本次课程设计中，我的题目是音乐信号的滤波去噪，相对平时作业题更贴近现实。不像作业中设计滤波器最多是处理余弦信号等，处理音乐信号让我更形象地认识到滤波器在实际生活中的应用，也增强了学习兴趣。还有就是在写课程设计报告时，你需要对原先课堂上学的知识有更透彻的认识，能对其原理和来龙去脉一清二楚。所以写报告在一定程度上也是整理思路，加深对课程的理解。

经过几天的设计和实践，已经比较顺利地完成了本次课程设计，在此感谢胡老师的指导和一些同学的帮助！

6 参考文献

[1] [美]维纳·K·英格尔，约翰·G·普罗克斯.刘叔棠.数字信号处理（MATLAB版）Digital Signal Processing Using MATLAB.第2版.西安：西安交通大学出版社，2008年

[2] 张志涌等.精通MA TLAB 6.5版.北京：北京航空航天大学出版社，2003年

[3] 邵玉斌.Matlab/Simulink通信系统建模与仿真实例分析. 北京：清华大学出版社，2008

[4] 邓华．MATLAB通信仿真及应用实例详解．北京：国防工业出版社，2003

附录（程序用到的两个自编函数）：

function [d1,d2] = db2delta(Rp,As)

% 从相对dB标度变换为绝对delta 标度.

% ----------------------------------------------------------

% [d1,d2] = db2delta(Rp,As)

% d1 = 通带容差

% d2 = 阻带容差

% Rp = 通带波动

% As = 阻带衰减

%

K = 10^(Rp/20);

d1 = (K-1)/(K+1);

d2 = (1+d1)*(10^(-As/20));

function [y,ny] = conv_m(x,nx,h,nh)

% 信号处理的改进卷积程序

% --------------------------------------------------

% [y,ny] = conv_m(x,nx,h,nh)

% y = 卷积结果

% ny = y 的基底(support)

% x =基底nx 上的第一个信号

% nx = x 的支架

% h =基底nh上的第二个信号

% nh = h 的基底

%

nyb = nx(1)+nh(1); nye = nx(length(x)) + nh(length(h)); ny = [nyb:nye];

y = conv(x,h);