音乐信号滤波去噪—使用椭圆滤波器 课程设计内容

音乐信号滤波去噪——使用脉冲响应不变法设计的椭圆滤波器

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

摘要本设计是用麦克风采集一段音乐信号，绘制波形并观察其频谱，给定相应技术指标，用脉冲响应不变法设计的一个满足指标的椭圆IIR滤波器，对该音乐信号进行滤波去噪处理，比较滤波前后的波形和频谱并进行分析，根据结果和学过的理论得出合理的结论。与不同信源相同滤波方法的同学比较各种信源的特点，与相同信源不同滤波方法的同学比较各种滤波方法性能的优劣。

关键词MATLAB7.0；数字信号处理；滤波去噪；频谱；相位；IIR滤波器噪声

1 引言

本课程设计主要是设计一个音乐信号滤波去噪——使用脉冲响应不变法设计的椭圆IIR滤波器。对该音乐信号进行滤波去噪处理，比较滤波前后的波形和频谱并进行分析。

1.1 课程设计目的

数字信号处理是一门以算法为核心，理论和实践性较强的学科。是电子信息工程、通信工程专业、电子信息科学与技术专业的一门重要的专业技术基础课。数字信号处理课程是在学习完数字信号处理的相关理论后，进行的综合性训练课程，其目的是：

1．使学生进一步巩固数字信号处理的基本概念、理论、分析方法和实现方法；2．增强学生应用Matlab语言编写数字信号处理的应用程序及分析、解决实际问题的能力；

1.2 课程设计要求

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

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

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

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

1.3 设计平台

MATLAB的名称源自Matrix Laboratory，1984年由美国Mathworks公司推向市场。它是一种科学计算软件，专门以矩阵的形式处理数据。MATLAB将高性能的数值计算和可视化集成在一起，并提供了大量的内置函数，从而被广泛地应用于科学计算、控制系统、信息处理等领域的分析、仿真和设计工作。

MATLAB在信号与系统中的应用主要包括符号运算和数值计算仿真分析。由于信号与系统课程的许多内容都是基于公式演算，而MATLAB借助符号数学工具箱提供的符号运算功能能基本满足信号与系统课程的需求。例如，解微分方程、傅里叶正反变换、拉普拉斯正反变换、z正反变换等。MATLAB在信号与系统中的另一主要应用是数值计算与仿真分析，主要包括函数波形绘制、函数运算、冲激响应与阶跃响应仿真分析、信号的时域分析、信号的频谱分析、系统的S域分析、零极点图绘制等内容。数值计算仿真分析可以帮助学生更深入理解信号与系统的理论知识，并为将来使用MATLAB进行信号处理领域的各种分析和实际应用打下基础。

语言简洁紧凑，使用方便灵活，库函数极其丰富。MATLAB程序书写形式自由，利用其丰富的库函数避开繁杂的子程序编程任务，压缩了一切不必要的编程工作。由于库函数都由本领域的专家编写，用户不必担心函数的可靠性。可以说，用MATLAB进行科技开发是站在专家的肩膀上。

2 设计原理

2.1脉冲响应不变法

脉冲响应不变法是一种将模拟滤波器转化为数字滤波器的基本方法。它通过对模拟滤波器的单位冲激响应h（t）等间隔抽样来获取对应数字滤波器的单位脉冲响应h[k],即h[k]=h(t)|t=kT =h(kT)，其中T是等间隔抽样。

优点：数字滤波器和模拟滤波器的频率关系为线性（ω=ΩT）。如果Ha(s)是稳定的，即其极点在S左半平面，映射到H(Z)也是稳定的。

缺点：有频谱周期延拓效应，因此只能用于带限的频响特性，如衰减特性很好的低通或带通。而高频衰减越大，频响的混淆效应越小，至于高通和带限滤波器,由于它们在高频部分不衰减，因此将完全混淆在低频响应中。所以用脉冲响应不变法实现高通和带限滤波器时，应增加一保护滤波器，滤掉高于折叠频率以上的频带，然后再用脉冲响应不变法转换为数字滤波器，这会增加设计的复杂性和滤波器的阶数，只有在一定需要频率线性关系或保持网络瞬态响应时才采用。

变换原理：

数字滤波器的单位脉冲响应h(n)模仿模拟滤波器的冲激响应h a(t)，使h(n)等于h a(t)的采样值，即

h(n)＝ha(nT) （2-1）

设模拟滤波器的传输函数H a(s)只有单阶极点，且分母的阶数高于分子阶数，则

（2-2）其拉普拉斯反变换

（2-3）对h a(t)采样，得到数字滤波器的单位脉冲响应为

（2-4）对h(n)取z变换，即得到数字滤波器的系统函数

（2-5）可见，s平面的极点si映射到z平面的极点，而对应系数不变。

z平面与s平面的映射关系：

（2-6）

令，，代入上式，得

（2-7）

因此，

z平面与s平面的映射关系如图2.1所示。

图2.1 z平面与s平面的映射关系

讨论：

（1）σ=0，r=1，表明s平面虚轴映射为z平面的单位圆。

（2）σ<0，r<1；σ>0，r>1。表明s左半平面映射为z平面的单位圆内部，而s右半平面则映射为z平面单位圆外部。

（3）由于ω=ΩT，当ω自0至±π变化时，Ω自0至±π/T变化。

s平面上每一条宽为2π/T的横条，都将重叠的映射到整个z平面上。反映了是和的周期延拓之间有变化关系，脉冲响应不变法不是从s 平面到z平面的简单代数映射关系。

虽然脉冲响应不变法能保证S平面与Z平面的极点位置有一一对应的代数关系，但这并不是说整个S平面与Z平面就存在这种一一对应的关系，特别是数字滤波器的零点位置与S平面上的零点就没有一一对应关系，而是随着Ha(s)的极点si与系数Ai的不同而不同。

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

数字滤波器的频率响应和模拟滤波器频率响应间的关系：

由，可得

数字滤波器的频率响应和模拟滤波器频率响应间的关系如图2.2所示。

图2.2数字滤波器的频率响应和模拟滤波器频率响应间的关系

数字滤波器的频响是模拟滤波器频响的周期延拓。只有当模拟滤波器的频响是带限的，且带限于折叠频率以内，即

，（2-8）

有

，（2-9）

此时，数字滤波器的频率响应重现模拟滤波器的频率响应而不产生混叠失真。

2.2 IIR数字滤波器

IIR(Infinite Impulse Response)数字滤波器，又名“无限脉冲响应数字滤波器”，或“递归滤波器”。递归滤波器，也就是IIR数字滤波器，顾名思义，具有反馈，一般认为具有无限的脉冲响应。

IIR滤波器有以下几个特点：

(1)封闭函数 IIR数字滤波器的系统函数可以写成封闭函数的形式。

(2)IIR数字滤波器采用递归型结构。IIR数字滤波器采用递归型结构，即结构上带有反馈环路。IIR滤波器运算结构通常由延时、乘以系数和相加等基本运算组成，可以组合成直接型、正准型、级联型、并联型四种结构形式，都具有反馈回路。由于运算中的舍入处理，使误差不断累积，有时会产生微弱的寄生振荡。

(3)借助成熟的模拟滤波器的成果。IIR数字滤波器在设计上可以借助成熟的模拟滤波器的成果，如巴特沃斯、契比雪夫和椭圆滤波器等，有现成的设计数据或图表可查，其设计工作量比较小，对计算工具的要求不高。在设计一个IIR 数字滤波器时，我们根据指标先写出模拟滤波器的公式，然后通过一定的变换，将模拟滤波器的公式转换成数字滤波器的公式。

(4)需加相位校准网络。IIR数字滤波器的相位特性不好控制，对相位要求较高时，需加相位校准网络。

IIR数字滤波器的设计：

利用MATLAB信号处理工具箱中的滤波器设计和分析工具(FDATool)可以很方便地设计出符合应用要求的未经量化的IIR数字滤波器。需要将MATLAB设计出的IIR数字滤波器进一步分解和量化，从而获得可用FPGA实现的滤波器系数。IIR数字滤波器的设计步骤:

由于采用了级联结构，因此如何将滤波器的每一个极点和零点相组合，从而使得数字滤波器输出所含的噪声最小是个十分关键的问题。为了产生最优的量化后的IIR数字滤波器，采用如下步骤进行设计。

(1)首先计算整体传递函数的零极点；

(2)选取具有最大幅度的极点以及距离它最近的零点，使用它们组成一个二阶基本节的传递函数；

(3)、对于剩下的极点和零点采用与3.2相类似的步骤，直至形成所有的二阶基本节。

通过上面三步法进行的设计可以保证IIR数字滤波器中N位乘法器产生的量化舍入误差最小。

2.3 椭圆滤波器

椭圆滤波器（Elliptic filter ）又称考尔滤波器（Cauer filter ），是在通带和阻带等波纹的一种滤波器。椭圆滤波器相比其他类型的滤波器，在阶数相同的条件下有着最小的通带和阻带波动。它在通带和阻带的波动相同，这一点区别于在通带和阻带都平坦的巴特沃斯滤波器，以及通带平坦、阻带等波纹或是阻带平坦、通带等波纹的切比雪夫滤波器。

椭圆低通滤波器：

定义：幅度平方响应满足下式的滤波器

（2-10）

N: 滤波器阶数;

ε : 通带波纹；

UN(x)： N 阶 Jacobian 椭圆函数 特性：通带和阻带都呈现等波纹特性 在给定指标下能实现最小阶次

（2-11）

（2-12）

（2-13）

第一类完全椭圆积分

可由MATLAB 提供的ellipke 函数计算 [z,p,k]=ellipap(N,Rp,As);

设计阶次为N ，通带波纹为Rp ，阻带衰减As 的归一化椭圆模拟原型滤波器

???

? ??ΩΩ+=

Ωc N

a U j H 22211

|)(|ε()(

)

2

12

1

12/2

2

2

()1,()1,,

1()1sin p s

K k K

k N K k K k k k A d K x x πε

θ

θ

-=

-Ω=

=

Ω-=-?

[b,a] = u_ellipap(N,Rp,As,Omegac)

设计非归一化椭圆滤波器

[b,a] = afd_elip(Wp,Ws,Rp,As)

设计椭圆IIR滤波器满足技术指标：fb=3400;

fc=3500;

As=75;

Rp=0.25;

T=1;

3设计步骤

3.1 设计流程图

语音信号滤波去噪——使用脉冲响应不变法设计的椭圆滤波器的设计流程

图3.1 设计流程图

开始 语音信号的采集（wavread 函数），画时域图,同时作快速傅里叶变换，并且画频设定滤波器性能指标，通带截止频率fb=3400，阻带截止频率fc=3500，通带波纹Rp=0.25，阻带波纹As=75 用脉冲响应不变法设计椭圆滤波器，验证并进行频谱分析 比较滤波前后语音信号的波形及频谱 回放语音信号 结束 设计好的滤波器进行滤波处理

3.2 语音信号的采集

用麦克风采集一段音乐信号，在录音的过程中人为的加入噪声。时间在30s 左右。然后在Matlab软件平台下，利用函数wavread对语音信号进行采样，记住采样频率和采样点数。录音机使用如图3.2所示。

图3.2 录音机使用及其属性图

3.3.产生原始语音信号

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

[x,Fs,bits]=wavread('d:\binbin.wav'); % x：语音数据； fs：采样频率；

bits：采样点数

sound(x,Fs,bits); %话音回放

N=length(x);

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

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

x=x';y1=x+0.15*sin(fn*2*pi*t);

sound(y1,Fs,bits);

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

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

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

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

其中，[y,Fs,bits]=wavread(‘'D:\binbin.wav'’); % 输入参数为文件的全路径和文件名，输出的第一个参数是每个样本的值。

sound(y,Fs,bits);用于对声音的回放。向量y则就代表了一个信号（也即

一个复杂的“函数表达式”）也就是说可以像处理一个信号表达式一样处理这

个声音信号。

程序运行后，可以听到清楚的音一段乐，同时可以得到音乐的信号波形图。

图3.3原始音乐信号及加入噪声后音乐信号的比较图

3.4 椭圆IIR滤波器的设计

fb=3400;fc=3500;

As=75;Rp=0.25; %音乐信号处理时采用的滤波器性能指标：

T=1;

wp=(fb/Fs)*2*pi;

ws=(fc/Fs)*2*pi;

OmegaP=wp/T;OmegaS=ws/T; %按频率转换公式进行转换[cs,ds]=afd_elip(OmegaP,OmegaS,Rp,As); %计算模拟滤波器的最小阶数[b,a]=imp_invr(cs,ds,T); %低通转换为低通滤波器的频率转化

[db,mag,pha, grd,w]=freqz_m(b,a);

*** 椭圆滤波器阶次 = 14

图3.4 利用脉冲响应不变法设计的椭圆滤波器（w（单位：错误！未找到引用源。）（a）幅度响应；（b）相位响应；（c）幅度（dB）；（d）群延时

3.5 音乐信号的滤波和去噪

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

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

f=0:deltaf:Fs/2-deltaf;

t=0:1/Fs:(N-1)/Fs;

图3.5加入噪声后音乐信号和滤波后音乐信号的比较图

图3.6原始音乐信号、加入噪声后音乐信号、滤波后音乐信号的比较图

回放音乐信号

sound (y_fil,Fs,bits);

3.6实验结论：

给一段音乐信号加入噪声，然后用脉冲响应不变法设计的椭圆IIR滤波器滤波，将加在音乐信号中的噪声滤除，用函数sound可以对声音进行回放。其调用格式：sound (y_fil,Fs,bits);滤波前可以明显听出有尖锐的单频啸叫声，经过滤波器滤波后噪声被滤除了。可以感觉滤波前后的声音有明显变化。可听到清楚且没有噪声的一段音乐。

4 出现的问题及解决方法

在这次课程设计中我设计的滤波器对语音信号的滤波处理或多或少会因为我的疏忽或实践能力的不够还是会遇到很多问题的，当遇到问题时我不但运用自己所学到的知识还得到同学和老师的帮助将其解决，不断修善和修改。设计中出现的问题及其解决方法如下。

（1）我的是进行音乐信号滤波去噪——使用脉冲响应不变法设计的椭圆滤波器，我开始错误的人为加入噪声，频率为2100Hz。结果滤波前后的频谱一致，后来，用高频余弦函数加入噪声。才得出正确的结果。

（2）在运行matlab设计的滤波器时，当编辑完程序时无法放出声音，后来发现我应当把采集的音乐信号wav文件放到matlab的work文件夹中。

（3）要区分清楚ft,fs,fp所代表的意思。

5结束语

经过两个星期的课程设计实习，我有了很大的感受。不同于平时学习，做作业，课程设计每个人都有一个自己的题目，完全得靠自己，不能有依赖心理。把平时书本的知识应用在了实践中。

通过本次课程设计，真正加强了我对课本知识的理解，有机会去把课本中的知识运用到实践中，我了解很多关于专业的知识，以前每次学这些时，总是不知道这些东西具体拿来有什么用，我这才发现，原来身边有个这么好的工具，而且通过这么些天的研究我知道了关于simulink的一些基本知识与用途。巩固了所学理论，培养了初步的实际工作能力和设计技术能力。

总的来说，这次课程设计过程还是比较愉快轻松的，虽然中间有过一些困难，但是在老师与同学的指点下我还是渡过了，而经过人生中的第一次课程设计，我相信我以后会越做越好的。

Maltab无疑是一款功能超强的数学软件，基于它可实现多领域的科学研究和开发。在本次课程设计中,收获最多的是如何看待问题以及如何解决问题。一

个很现实的问题就是搞科研开发，遇到问题是必不可免，关键是如何看待这样的问题，持怎样的心态，持怎样的目的。再就是问题是遇到了，要如何解决以及用怎样最优的方案去处理，都得下心思去思考和摸索。再一个很现实的问题是，要学会学习，学会自学。毕竟社会日新月异，科技蓬勃发展，新领域新知识不断涌现，怎么样才能在社会保持有效的立足，怎样才能翻弄潮流，面对社会变迁做到运筹帷幄之中呢？在我看来，关键要保持良好心态，保持一种积极永不妥协永不放弃的信念，不断进取不断开拓。

在学好基本课内的专业知识之外，应尽量多补充课外的和本专业相协调的专业知识。强势未来就业利剑。

参考文献

（1）维纳.K.英格尔、约翰.G.普罗克斯主编《数字信号处理》（MATLAB版）——西安交通大学出版社；

（2）孙志涌编著《精通MATLAB6.5版》——北京航空航天工业出版社

（3）孔屹主编《MATLAB通信仿真开发手册》——国防工业出版社；

(4) 李文军，李师贤，周晓聪. C++作为计算机专业程序设计入门语言的实践与探讨. 计算机科学，1999，26（4）：80～83

(5) F.Brokken and K.Kubat. C++ Annotations. Version 4.4.0m，ICCE，University of Groningen，Netherlands，1990. 250～280

附件：

产生原始语音信号程序：

[x,Fs,bits]=wavread('d:\binbin.wav'); % x：语音数据； fs：采样频率；bits：采样点数

sound(x,Fs,bits); %话音回放

N=length(x);

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

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

x=x';y1=x+0.15*sin(fn*2*pi*t);

sound(y1,Fs,bits);

Y=abs(fft(x));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,x); %画出原始音乐信号的波形xlabel('时间t');

ylabel('幅度x(t)');

title('原始音乐信号时域波形');

axis([0,8,-1,1])

grid on

subplot(2,2,2);plot( f,Y);

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

axis([0,4000,0,1000]);grid on

subplot(2,2,3);plot(t,y1); axis([0,8,-1,1]);grid on

xlabel('时间(单位:s)');ylabel('幅度');title('加入噪声后的音乐信号时域波形'); subplot(2,2,4);plot(f,Y1);axis([0,4000,0,1000]);grid on

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

椭圆IIR滤波器的设计

音乐信号处理时采用的滤波器性能指标：

fb=3400;fc=3500; As=75;Rp=0.25;

T=1;

wp=(fb/Fs)*2*pi;

ws=(fc/Fs)*2*pi;

OmegaP=wp/T;OmegaS=ws/T; %按频率转换公式进行转换

[cs,ds]=afd_elip(OmegaP,OmegaS,Rp,As); %计算模拟滤波器的最小阶数[b,a]=imp_invr(cs,ds,T); %低通转换为低通滤波器的频率转化

[db,mag,pha, grd,w]=freqz_m(b,a);

subplot(2,2,1);plot(w,mag);title('椭圆滤波器幅度响应')

xlabel('w');ylabel('幅值|H|'); axis([2,3.5,0,1.3])

grid on

subplot(2,2,2);plot(w,pha);title('椭圆滤波器相位响应')

xlabel('w');ylabel('相位（单位：π ）');

axis([1,3.5,-4,4])

grid on

subplot(2,2,3);plot(w,db);title('椭圆滤波器幅度(dB)')

xlabel('w');ylabel('dB');

axis([2,3.5,-100,50])

grid on

subplot(2,2,4);plot(w,grd);title('群时延')

xlabel('w');ylabel('样本');

axis([2.3,3,-200,200])

grid on

音乐信号的滤波和去噪

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