音乐信号滤波去噪 ——使用脉冲响应不变法设计的切比雪夫I型滤波器要点

音乐信号滤波去噪

——使用脉冲响应不变法设计的切比雪夫I

型滤波器

学生姓名：李柳指导老师：黄红兵

摘要本课程设计主要内容是使用脉冲响应不变法设计的切比雪夫I型滤波器，对一段音乐信号进行滤波去噪处理并根据滤波前后的波形和频谱分析滤波性能。本课程设计仿真平台为MATLAB7.0，开发工具是M语言编程。首先下载一段音乐，并人为加入一单频噪声，然后对信号进行频谱分析以确定所加噪声频率，并设计滤波器进行滤波去噪处理，最后比较滤波前后的波形和频谱并进行分析。由分析结果可知，滤波器后的音乐信号与原始信号基本一致，即设计的滤波器能够去除信号中所加单频噪声，达到了设计目的。

关键词滤波去噪；脉冲响应不变法；切比雪夫I型滤波器；MATLAB

1 引言

此次课程设计主要是在网上采集一段8000Hz，8位的单声道PCM格式音乐信号，并绘制波形观察其时域和频域的波形图，再在MATLAB平台上，将该音乐信号进行滤波去噪处理，对比滤波前和滤波后的时域和频域的波形图，根据结果和学过的理论得出合理的结论。

1.1 课程设计目的

课程设计有利于我们对基础知识的理解，并将所学的知识应用起来，此次课程设计用到Matlab，数字信号处理，以及办公软件Visio等知识，平时总是在分析滤波器，其实并不太理解滤波器跟我们的生活有什么联系，而课程设计要求我们自己动手操作，从原始信号的采集到加入噪声之后的信号到使用我们自己设计的滤波器对加噪信号进行滤波处理之后恢复出原始信号，这个过程让我们真正理解我们所学的知识在我们生活中的用处，从而让增强我们对这门学科以及我们专业的兴趣。另外此次课程设计也有利于逻辑

思维的锻炼，《数字信号处理》课程设计是在学生完成数字信号处理和MATLAB的结合后的基本实验以后开设的。本课程设计的目的是为了让学生综合数字信号处理和MATLAB并实现一个较为完整的小型滤波系统。这一点与验证性的基本实验有本质性的区别。平时的学习都是分模块进行，并没有系统的自己一个人独立完成设计到操作的过程，这样系统的设计正好锻炼了我们这方面的能力。开设课程设计课程的主要目的是通过系统设计、软件仿真、程序安排与调试、写课程设计报告等步骤，使学生初步掌握工程设计的具体步骤和方法，提高分析问题和解决问题的能力，提高实际应用水平。

1.2课程设计的要求

（1）滤波器指标必须符合工程实际。

（2）设计完后应检查其频率响应曲线是否满足指标。

（3）处理结果和分析结论应该一致，而且应符合理论。

（4）独立完成课程设计并按要求编写课程设计报告书。

1.3设计平台

此次课程设计在MATLAB平台下进行设计，MATLAB（矩阵实验室）是MATrix LABoratory的缩写，是一款由美国The MathWorks公司出品的商业数学软件。MATLAB 是一种用于算法开发、数据可视化、数据分析以及数值计算的高级技术计算语言和交互式环境。除了矩阵运算、绘制函数/数据图像等常用功能外，MATLAB还可以用来创建用户界面及与调用其它语言（包括C，C++和FORTRAN）编写的程序。MATLAB可以进行矩阵运算、绘制函数和数据、实现算法、创建用户界面、连接其他编程语言的程序等，主要应用于工程计算、控制设计、信号处理与通讯、图像处理、信号检测、金融建模设计与分析等领域。高效的数值计算及符号计算功能，能使用户从繁杂的数学运算分析中解脱出来，且具有完备的图形处理功能，实现计算结果和编程的可视化，友好的用户界面及接近数学表达式的自然化语言，使学者易于学习和掌握，功能丰富的应用工具箱(如信号处理工具箱、通信工具箱等) ，为用户提供了大量方便实用的处理工具。

2 设计原理

在网上采集一段音乐信号，绘制波形并观察其频谱，给定相应技术指标，用脉冲响应不变法设计的一个满足指标的切比雪夫I型IIR滤波器，并对该信号进行滤波去噪处理，

之后比较滤波前后的波形和频谱并进行分析。

2.1 IIR 滤波器

IIR 滤波器具有无限长脉冲响应，因此能够与模拟滤波器相匹敌；一般来说，所有的模拟滤波器都有无限长脉冲响应。因此，IIR 滤波器设计的基本方法是利用复值映射将大家熟知的模拟滤波器变换为数字滤波器。这个基本方法称为A/D （模拟-数字）滤波器变换。IIR 数字滤波器采用递归型结构，即结构上带有反馈环路。IIR 滤波器运算结构通常由延时、乘以系数和相加等基本运算组成，可以组合成直接型、级联型、并联型三种结构形式，都具有反馈回路。通常有IIR 数字滤波器的直接和间接设计法，所谓模拟滤波器设计数字滤波器，第一步将设计的归一化样本模拟低通滤波器经模拟—模拟频带变换法转换成模拟低通滤波器，第二步，然后数字化采用脉冲响应不变法成各相应频带的数字滤波器[2]。

2.2切比雪夫I 型滤波器

切比雪夫滤波器，又名“车比雪夫滤波器”，是在通带或阻带上频率响应幅度等波纹波动的滤波器。切比雪夫滤波器来自切比雪夫分布，以“切比雪夫”命名，是用以纪念俄罗斯数学家巴夫尼提·列波维其·切比雪夫。切比雪夫滤波器在过渡带比巴特沃斯滤波器的衰减快，但频率响应的幅频特性不如后者平坦。切比雪夫滤波器和理想滤波器的频率响应曲线之间的误差最小，但是在通频带内存在幅度波动。I 型切比雪夫滤波器在通带（或称“通频带”）上频率响应幅度等波纹波动的滤波器称为“I 型切比雪夫滤波器；n 阶第一类切比雪夫滤波器的幅度与频率的关系可用下列公式表示：

201

()|()|n n n G w H jw w w ==?? ??? （2.1） 是滤波器在截止频率：的放大率arccos );0n Ω)arccos ?Ω

2.3脉冲响应不变法

脉冲响应不变法是从滤波器的脉冲响应出发，使数字滤波器的单位脉冲响应序列h(n)逼近模拟滤波器的冲击响应ha(t), 使h(n)正好等于ha(t)的抽样值，即满足

()()a h n h nT = （2.6）

其中，T 为采样周期[2]。

如以Ha(s)及H （z ）分别表示ha(t)的拉氏变换及h(n)的z 变换，

设 ()()a a h s L h t =????

()()h z Z h n =????

且设()a H s 只有单级阶点，且假定()a H s 的分母阶次N>M ，则根据采样序列z 变换与模拟信号拉氏变换的关系，得:

()12|st a z e m H z H s j m T T π∞=→∞??=+ ???

∑ （2.7） 上式表明，采用脉冲响应不变法将模拟滤波器变换为数字滤波器时，它所完成的S 平面到Z 平面的变换，正是以前讨论的拉氏变换到Z 变换的标准变换关系，即首先对Ha(s)作周期延拓，然后再经过z=e st 的映射关系映射到Z 平面上[2]。

z=e st 的映射关系表明，S 平面上每一条宽为2π/T 的横带部分，都将重叠地映射到Z 平面的整个全部平面上。每一横带的左半部分映射到Z 平面单位圆以内，每一横带的右半部分映射到Z 平面单位圆以外，j Ω轴映射在单位圆上，但j Ω轴上的每一段2π/T 都对应于绕单位圆一周，如下图所示：

(a)s 平面 (b) z 平面

图2-2 脉冲响应不变法的映射关系

应当指出，Z=e st 的映射关系反映的是H a (s)的周期延拓与H （z ）的关系，而不是H a (s)本身与H （z ）的关系，因此，使用脉冲响应不变法时，从H a (s)到H(z)并没有一个由S 平

面到Z 平面的简单代数映射关系，即没有一个s=f(z)的代数关系式。另外，数字滤波器的频响也不是简单的重现模拟滤波器的频响，而是模拟滤波器频响的周期延拓，周期为ΩS =2π/T=2πf s ，即

()1212jw

a a m m n w n H e H j j H j T T T T ππ∞∞→∞→∞+????=Ω+= ? ?????∑∑ （2.8） 脉冲响应不变法特别适用于用部分分式表达的传递函数，模拟滤波器的传递函数若只有单阶极点，且分母的阶数高于分子阶数N ＞M ，则可表达为部分分式形式：

1()N i a i i

A H s s s ==-∑ 其拉氏反变换为：

()1()i N

s t a i

i H t Ae u t ==∑ 其中u(t)为单位阶跃函数。 对ha(t)采样就得到数字滤波器的单位脉冲响应序列，再对h(n)取Z 变换，得到数字滤波器的传递函数：

()()()11110

i i n N N s T s T n i i n i i n H z A e

Z A e z ∞∞--======∑∑∑∑ （2.9）

第二个求和为等比级数之和，

()111|1i i k s T k s T e z e z -→∞---，要收敛的话，必有()1|0i k s T k e z -→∞=，所以

有: ()111i N i s T i A H z e z -==-∑ （2.10）

比较部分分式形式的Ha(s)和上式H(z)可以看到，S 平面上的极点s=s i ,变换到Z 平面

上是极点i s T i z e =,而Ha(s)与H(z)中部分分式所对应的系数不变。如果模拟滤波器是稳定的，则所有极点s i 都在S 左半平面，即Re[s i ]＜0,那么变换后H(z)的极点i s T i z e =也都在单位圆以内，即[]Re ||1i i s T s T e e =<，因此数字滤波器保持稳定。

虽然脉冲响应不变法能保证S 平面与Z 平面的极点位置有一一对应的代数关系，但这并不是说整个S 平面与Z 平面就存在这种一一对应的关系，特别是数字滤波器的零点位,

与S 平面上的零点就没有一一对应关系，而是随着()a H s 的极点i s 与系数i A 的不同而不同。

()jw H e 是()a H j Ω的周期延拓（周期为fs ），因()a H j Ω并不是带限，即在超过fs 频率部分并不为0，所以就产生了混迭。当为低通或带通滤波器时， fs 越大，则()a H j Ω的下一周期相隔越远，混迭也就越小。当为带阻或高通滤波器时，()a H j Ω在超过/2fs 频率部分全为通带，这样就不满足抽样定理，发生了完全的混迭，所以脉冲响应不变法不能设计带阻或高通滤波器。

3设计步骤

3.1 设计流程图

认真阅读老师下发的任务书，上网搜集查找相关材料，大致的画出设计流程图，即音乐信号滤波去噪—使用脉冲不变响应法设计的切比雪夫I 型滤波器的设计流程如图3-1所示：

图3-1 设计流程图

3.2 录制音乐信号

在网上下载一段音乐信号，利用windows下的录音机工具通过点击录音机界面上的文件—属性，在弹出的对话框中选择马上转换，设置录音格式为PCM编码，属性为8kHz，单声道，这样的设置在相同的时间下数据最少，处理时间也就最短。点击文件，保存，存为全英文的文件名，后缀为“.wav”。如“liliu.wav”。并将其放入MATLAB安装盘下的work 文件夹中就好了。如下图所示：

图3-2 音乐信号格式转换图

8000Hz的采样率是指：声音信号在“模→数”转换过程中单位时间内采样8000次。采样值是指每一次采样周期内声音模拟信号的积分值。对于单声道声音文件，采样数据为八位的短整数，wave文件是计算机领域最常用的数字化声音文件格式之一，它是微软专门为Windows系统定义的波形文件格式（Waveform Audio），由于其扩展名为"*.wav"。PCM也被称为脉码编码调制。PCM中的声音数据没有被压缩，如果是单声道的文件，采样数据按时间的先后顺序依次存入,它的基本组织单位是BYTE(8bit)。

采集信号完成以后，利用MATLAB画出原始信号在时域和频域的波形图，编写M文

件如下：

[x,fs,bits]=wavread('liliu.wav'); % 输入参数为文件的全路径和文件名，输出的第一个参数是每个样本的值，fs是生成该波形文件时的采样率，bits是波形文件每样本的编码位数。

sound(x,fs,bits); % 按指定的采样率和每样本编码位数回放

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

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

X=abs(fft(x)); % 对原始信号进行fft变换,取幅度谱

X=X(1:N/2); % 截取前半部分

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

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

得到的原始信号的时域与频域的波形图如下图所示：

图3-3 原始音乐信号时域与频域的波形

画出音乐信号的时域波形；然后对音乐号进行快速傅里叶变换，得到信号的频谱特性，从原始音乐信号的幅度谱中可以看出，音乐信号的频率主要集中在0～1000Hz ，2000Hz 以后的频率段几乎没有能量分布。

采集到合适格式的音乐信号之后，对原始音乐信号进行加噪处理，然后对加入单频率的噪声的音乐信号进行fft变换,取幅度谱，画出加噪的音乐信号在时域和频域的波形，并

对比进行分析，主要程序如下：

x=x'; y=x+0.1*sin(fn*2*pi*t); %加噪信号

sound(y,fs,bits); %加噪信号声音

Y=abs(fft(y)); %对原始信号和加噪信号进行fft变换,取幅度谱

Y=Y(1:N/2); % 截取前半部分

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

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

运行程序的得到的加噪信号的时域与频域的波形图如下所示：

图3-4 加噪信号的时域与频域的波形图

由于此次课程设计加入的噪声为fn=2000Hz的单频信号，运行程序的时候可以明显的听到加噪之后的音乐信号对比原始信号有明显的单频呼啸声，由运行结果图3-4右下角的频谱图也可明显的看到在频率为2000Hz时存在单频噪声，时域图也有明显的变化。

3.3 滤波器设计

我们这次设计的主要任务是对单频干扰进行滤波处理，因此设计一个低通滤波器，设计指标如下：

取fp=fn-200=1800 fc=fn-50=1950，通带波纹1p R dB =，阻带衰减40s A dB =，用脉冲响应不变法设计的切比雪夫I 型滤波器对单频干扰进行滤波处理。

设计滤波器的重要程序如下：

fp=fn-200;fc=fn-50; %定义通带和阻带截止频率

Rp=1;As=40; % 定义通带波纹和阻带衰减

wp=fp/fs*2*pi;ws=fc/fs*2*pi; %计算对应的数字频率

T=1; %定义采样间隔

Omegap=wp/T; Omegas=ws/T; %模数指标变换

[N,Omegac]=cheb1ord(Omegap,Omegas,Rp,As,'s') %计算切比雪夫滤波器阶数和截止频率

[cs,ds]=cheby1(N,Rp,Omegac,'s');% 计算切比雪夫滤波器系统函数分子分母系数

% 第三步 模数滤波器变换

[b,a]=impinvar(cs,ds,1/T) % 脉冲响应不变法得到数字滤波器系统函数分子分母系数

[db,mag,pha,grd,w]=freqz_m(b,a); %验证滤波器是否达到指定性能

delta=[1,zeros(1,99)];ha=filter(b,a,delta);% 计算脉冲响应

在MATLAB 中运行滤波器的程序，得到滤波器的特性图，如下图所示：

图3-5 滤波器特性图

由程序运行结果可知wp/pi=0.450 ，ws/pi =0.4875，N =15，Omegac =1.4137,对比图3-4中滤波器的幅度特性图可知当w/pi>0.4875时，dB都大于40，所以利用脉冲响应不变法设计的切比雪夫滤波器符合要求。

3.4 信号滤波处理

在3.3节中已经用使用脉冲响应不变法设计的切比雪夫I型滤波器，现在只需用自己设计的切比雪夫I型滤波器对采集的信号进行滤波，在Matlab中，IIR滤波器利用函数filter 对含噪信号y进行滤波去噪处理，得到滤波后的信号存储为y1，并计算y1的频谱。部分程序如下：

y1=filter(b,a，y); % 用设计好的滤波器对y进行滤波处理,得到时域图

Y1=fft(y1);Y1=Y1(1:length(Y1)/2); % 计算频谱取前一半

对比滤波前后的时域与频域的波形图3-5，具体分析滤波是否达到要求。

图3-6 滤波前后的时域与频域的波形对比图

我们观察到图3-5滤波前后音乐信号的波形对比图，发现在时域波形中滤波后的波形与加干扰噪声后有明显的变化，与原始音乐信号几乎相同，说明设计的滤波器达到要求。而且在频域波形中，很明显地反应出设计的滤波器滤去了我们采集的音乐信号中的噪声，即滤去了频率高于2000Hz的信号。

3.5 结果分析

综合上述的实验步骤及其实验得出的结果图，在图3－4中我们可以发现加入的

2000Hz单频噪声是直接累加在原始音乐信号上的，属于加性噪声，在图3-5 滤波器特性图中可以得出，我们利用脉冲响应不变法设计的切比雪夫I滤波器属于低通滤波器，能将2000Hz以上的噪声滤掉，分析图3－6滤波前后的时域与频域的波形对比图，我们就能明显看出我们利用脉冲响应不变法设计的切比雪夫I滤波器滤去了我们采集的音乐信号中的噪声，达到了设计要求。

3.6 滤波器结构设计

IIR滤波器可以组合成直接型、级联型、并联型三种基本结构形式，都具有反馈回路。根据设计的滤波器的原始代码的运行结果，我们可以得知滤波器的阶数N＝15，相对来说滤波器的阶数不是特别大，这里我们采用直接型的结构。滤波器的结构图如下图所示：

图3－6 滤波器结构图

由滤波器的源代码的运行结果可以知道系统函数分子分母系数b 和a 的值，如下图所示：

图3－7数字滤波器系统函数分子分母系数

根据系统函数的分子分母系数以及IIR 滤波器直接型结构的特点很容易可以画图3－6滤波器结构图，N 阶滤波器需要N 个延时单元，左侧从上到下是分母系数的负数0~N a a --,右侧是分子系数0~N b b --。

4出现的问题及解决方法

这次的课程设计中我们设计的滤波器对音乐信号的滤噪处理或多或少会因为我们不仔细和实践能力不足会遇到各种问题，当遇到问题时我们运用自己学到的知识以及同学和老师的帮助将其解决，不断完善和修改。设计中出现的问题及解决方法有以下几个方面：

（1）在设计的过程中，刚开始从网上下载的音乐信号是立体声，不符合课程设计的要求，经过和同学的讨论，利用WINDOWS 下的录音机工具通过点击录音机界面上的文件—属性，在弹出的对话框中选择马上转换，设置录音格式为PCM 编码，属性为8kHz ，单声道这样就能符合原始信号的要求了。

（2）加噪后的信号的幅度谱图，单频噪声幅度远远高于其他频谱分量，导致其他频谱分量幅度过低而看不清楚加噪前后的频谱关系，通过使用axis 调整纵坐标的范围即可解决。

（3）其次在设滤波器的参数时取fp=fn-200=1800 fc=fn-50=1950，通带波纹Rp=1dB 阻带衰减As=16dB ，滤波后发现不能完全将噪声滤掉，仍然存在噪声，在设置的通带截止频率和阻带截止频率的差值太大就是太小，都不能完全将噪声滤掉，经过一次次调整与修改最终将参数设置为fp=fn-200=1800 fc=fn-50==1950，通带波纹Rp=1dB ，阻带衰减

As=40dB，这样就能更好的将噪声滤除，基本达到指标要求。

（4）刚开始设计好滤波器是，不知道怎样才是符合要求，将设计好的滤波器拿给老师看，老师说在Wp/pi,Ws/pi以及Rp,As，出各画一条直线能更好的看出滤波器是否满足要求，即，当w/pi大与ws/pi时都大于As=40dB就说明满足要求了。

（5）最后就是在使用Visio软件时的问题，由于是第一次使用，比较陌生，最开始使用的线是尾端带箭头的，这样每次画一段就需要两次修改，最后仔细的在工具箱中寻找，发现并使用了中间带箭头的线，这样画起图来方便很多。

5 结束语

通过这段时间的亲身经历，我感觉自己学到了：收集、整理资料、共同协作、分析及处理问题等许多方面的知识。在做设计的过程中，我有许多不懂得地方，在老师的指导下我一步步的解决问题完成论文，在完成过程中老师指导我去怎么选择资料，如何去利用网络资源，在这个学习的过程中，我了解到MATLAB的实用价值，同时利用MATLAB简单编程的特点，可以实现较为复杂的系统，因此MATLAB在通信系统仿真方面具有强大的功能和优越性，在课程设计的过程中明显觉得自己对MATLAB软件的不熟练，很多函数不会使用，而这次的课程设计让我深刻认识到MATLAB的实用性，同时也更加熟练的使用MATLAB，掌握了MATLAB的更多用法。同时这次课程设计是《数字信号处理》的，在设计的过程中明显发现对课本知识的不熟练，用到的很多课本知识仍要去翻阅课本，这让我认识到在以后的学习中，应该更加认真更加深入，不能浅尝辄止。在以后的学习或者其他设计中，应该提前做好充足的准备，对要做的东西有足够的了解之后再开始进行，其次以后在工作中要脚踏实地，在学术上要严谨，在思维上要活跃，在学业上要勤奋刻苦，不能懵懂的就以为自己已经会了。

最后我真诚感谢这期间老师给予我的全力帮助，细心指导以及对我的严格要求，是老师在我遇到问题时，不辞辛苦帮我解决，感谢老师在设计和任务安排上长时间的指导。

参考文献

[1] 张圣勤.MATLAB7.0实用教程[M].北京：机械工程出版社,2006.P140-180

[2] 程佩青.数字信号处理教程[M].北京：清华大学出版社,2014,P269-463

[3] 张志美. MATLAB完全自学手册[M].电子工业出版社,2013,P142-405

[4] 邵毅,张明玉,温艳,基于MATLAB的IIR数字滤波器的设计[J].宿州学院学报,2009年02期

[5] 王艳文,杨楠, 基于MATLAB/脉冲响应不变法的切比雪夫I型数字滤波器设计[J]. 新乡医学院生命科学技术学院.2013年18期

附录：音乐信号滤波去噪设计源程序清单

% 程序名称：liliu.m

% 程序功能：音乐信号滤波去噪——使用脉冲响应不变法设计的切比雪夫I型滤波器

% 程序作者：李柳

% 最后修改日期：2015-3-13

[x,fs,bits]=wavread('liliu.wav'); % 输入参数为文件的全路径和文件名，输出的第一个参数是每个样本的值，fs是生成该波形文件时的采样率，bits是波形文件每样本的编码位数。sound(x,fs,bits); % 按指定的采样率和每样本编码位数回放

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

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

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

x=x'; 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;% 计算频谱频率范围

fn=2000; % 单频噪声频率

fs=8000;

fp=fn-200;fc=fn-50; %定义通带和阻带截止频率

Rp=1;As=40; % 定义通带波纹和阻带衰减

wp=fp/fs*2*pi;ws=fc/fs*2*pi; %计算对应的数字频率

T=1; %定义采样间隔

% 第一步模数指标变换

Omegap=wp/T; Omegas=ws/T; %截止频率线性变换

% 第二步根据指标设计模拟滤波器

[N,Omegac]=cheb1ord(Omegap,Omegas,Rp,As,'s') %计算滤波器阶数和截止频率

[cs,ds]=cheby1(N,Rp,Omegac,'s');% 计算滤波器系统函数分子分母系数

% 第三步模数滤波器变换

[b,a]=impinvar(cs,ds,1/T) % 脉冲响应不变法得到数字滤波器系统函数分子分母系数

[db,mag,pha,grd,w]=freqz_m(b,a); %验证滤波器是否达到指定性能

delta=[1,zeros(1,99)];ha=filter(b,a,delta);% 计算脉冲响应

Rp = -min(db(1:fix(wp/(2*pi/1000))+1)) % 实际的通带波动

As = -round(max(db(fix(ws/(2*pi/1000))+1:501))) % 最小阻带衰减

y1=filter(b,a,y) % IIR滤波器对信号signal进行滤波处理，其中b,a为上面设计好的滤波器参数。

Y1=abs(fft(y1));Y1=Y1(1:length(Y1)/2); % 计算频谱取前一半

figure(1)

subplot(221);plot(t,x) % 绘制截取后的语音信号

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

subplot(223);plot(t,y) % 绘制截取后的语音信号

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

subplot(222);plot(f,abs(X)) % 绘制频谱图

title('原始音乐信号幅度谱图');xlabel('频率（单位：Hz）');ylabel('幅度谱') ;axis tight subplot(224);plot(f,abs(Y)) % 绘制频谱图

title('加入干扰后的音乐信号幅度谱图');xlabel('频率（单位：Hz）');ylabel('幅度谱') ;axis ([0 4000 0 400])

figure(2)

Subplot(221);plot(w/pi,db);axis([0 1 -50 20]);title(‘滤波器幅度响应');xlabel('w/pi');ylabel('dB');grid on;

line([wp/pi,wp/pi],[-50,20],'Color','r','LineWidth',2,'LineStyle','--');

line([Omegac/pi,Omegac/pi],[-50,20],'Color','r','LineWidth',2,'LineStyle','--');

line([ws/pi,ws/pi],[-50,20],'Color','r','LineWidth',2,'LineStyle','--');

line([0,1],[-Rp,-Rp],'Color','r','LineWidth',2,'LineStyle','--');

line([0,1],[-As,-As],'Color','r','LineWidth',2,'LineStyle','--');

Subplot(222);plot(w/pi,mag);axis([0 1 0 1.5]);title('滤波器幅度响应图');xlabel('w/pi');ylabel('幅度mag');grid on;

Subplot(223);plot(w/pi,pha);axis([0 1 -4 4]);title('滤波器相位响应图');xlabel('w/pi');ylabel('相位pha');grid on;

Subplot(224);plot(ha);axis([ 0 100 -0.2 0.6 ]);title('滤波器脉冲响应图');xlabel('n');ylabel('h（n）');grid on;

figure(3);

subplot(321);plot(t,x);axis tight;grid on;

title('原始音乐信号时间x');xlabel('时间t');ylabel('幅度');

% 绘制原始录制的音乐信号时间波形图

subplot(323);plot(t,y);axis tight;grid on;

title('加干扰音乐信号时间y');xlabel('时间t');ylabel('幅度');

% 绘制加噪音乐信号时间波形图

subplot(325);plot(t,y1);axis tight;grid on;

title('滤波后音乐信号时间y1');xlabel('时间t');ylabel('幅度');

% 绘制去噪音乐信号时间波形图

subplot(322);plot(f,abs(X));axis tight;grid on;

title('原始音乐信号幅度谱X');xlabel('频率f');ylabel('幅度');

% 绘制原始录制的音乐信号频谱图

subplot(324);plot(f,abs(Y));axis([0 4000 0 400]);grid on;

title('加干扰音乐信号幅度谱Y');xlabel('频率f');ylabel('幅度');

% 绘制加噪音乐信号频谱图

subplot(326);plot(f,abs(Y1));axis tight;grid on;

title('滤波后音乐信号幅度谱Y1');xlabel('频率f');ylabel('幅度')

% 绘制去噪音乐信号频谱图

sound(y1,fs,bits); %回放滤波去噪后的音乐

wavwrite(y1,fs,' liliuquzao.wav') %将去噪后音乐信号写成wav文件