专业知识341期

1、音乐信号滤波去噪使用双线性变换法设计的并联型椭圆滤波器第24页共24页音乐信号滤波去噪使用双线性变换法设计的并联型椭圆滤波器学生姓名 指导老师： 摘 要 本文简单、直观地介绍了椭圆低通滤波器的基本理论和设计思想，阐述了设计椭圆低通滤波器的具体步骤，并对仿真结果进行分析和处理。本课程设计主要内容是利用双线性变换法设计一个椭圆IIR滤波器，对一段音乐信号（含燥声）进行滤波去噪处理，对比滤波前后的波形和频谱图分析滤波器的性能。本课程设计仿真平台为MATLAB7.0，开发工具是M语言编程。首先在windows下用麦克风录制一段语音然后加入音乐背景，并人为加入一单频噪声，听加入噪声前后的音乐变化，然后对

2、信号进行频谱分析以确定所加噪声频率，然后按照设计滤波器进行滤波去噪处理，最后比较滤波前后的波形和频谱并进行分析，听滤波前后音乐的变化。由分析结果可知，滤波器后的语音信号与原始信号基本一致，即设计的IIR椭圆滤波器能够去除信号中所加单频噪声，达到了设计目的。关键词 滤波去噪；IIR滤波器；椭圆滤波器；MATLAB；双线性变换法1 引 言信号处理是科学研究和工程技术许多领域都需要进行的一个重要环节，传统上对信号的处理大都采用模拟系统实现。随着人们对信号处理要求的日益提高，以及模拟信号处理中一些不可克服的缺点，对信号的许多处理而采用数字的方法进行。近年来由于大规模集成电路和计算机技术的进步，信号的数

3、字处理技术得到了飞速发展。数字信号处理系统无论在性能、可靠性、体积、耗电量、成本等诸多方面都比模拟信号处理系统优越的多，使得许多以往采用模拟信号处理的系统越来越多地被数字处理系统所代替,数字信号处理技术在通信、语音、图像、自动控制、雷达、军事、航空航天、医疗和家用电器等众多领域得到了广泛的应用。在数字信号处理中,数字滤波器十分重要并已获得广泛应用,数字滤波器与模拟滤波器比较,具有精度高、稳定、体积小、重量轻、灵活、不要求阻抗匹配以及实现模拟滤波器无法实现的特殊滤波功能等优点。在各种滤波器中,椭圆滤波器具有其独特的优点1.1 课程设计目的本课程设计主要利用MATLAB，设计一个音乐信号滤波去噪使

4、用双线性变化法设计的并联型椭圆滤波器。数字信号处理课程设计是在学生完成数字信号处理和MATLAB的结合后的基本实验以后开设的。本课程设计的目的是为了让学生综合数字信号处理和MATLAB并实现一个较为完整的小型滤波系统。这一点与验证性的基本实验有本质性的区别。开设课程设计环节的主要目的是通过系统设计、软件仿真、程序安排与调试、写实习报告等步骤，使学生初步掌握工程设计的具体步骤和方法，提高分析问题和解决问题的能力，提高实际应用水平。1.2 课程设计的要求（1）滤波器指标必须符合工程实际。（2）设计完后应检查其频率响应曲线是否满足指标。（3）处理结果和分析结论应该一致，而且应符合理论。（4）独立完成

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

6、R滤波器从离散时间来看，若系统的单位抽样(冲激)响应延伸到无穷长，称之为“无限长单位冲激响应系统”，简称为IIR系统。无限长单位冲激响应（IIR）滤波器有以下几个特点:(1) 系统的单位冲激响应h(n)是无限长；(2) 系统函数H(z)在有限z平面（0<|z|<）；(3) 结构上存在着输出到输入的反馈，也就是结构上是递归型的。IIR滤波器采用递归型结构，即结构上带有反馈环路。同一种系统函数H（z）可以有多种不同的结构，基本网络结构有直接型、直接型、级联型、并联型四种，都具有反馈回路。同时，IIR数字滤波器在设计上可以借助成熟的模拟滤波器的成果，巴特沃斯(Butterworth)滤波

7、器、切比雪夫(Chebyshev)滤波器、椭圆(Cauer)滤波器、贝塞尔(Bessel)滤波器等，这些典型的滤波器各有特点。有现成的设计数据或图表可查，在设计一个IIR数字滤波器时，我们根据指标先写出模拟滤波器的公式，然后通过一定的变换，将模拟滤波器的公式转换成数字滤波器的公式。2.2 椭圆低通滤波器1.椭圆滤波器的概述常用数字滤波器的类型有巴特沃斯(Butterworth) ,切比雪夫(Chebyshev) 及椭圆型 滤波器，其中椭圆滤波器（Elliptic filter）又称考尔滤波器（Cauer filter），是一种性能优越的滤波器。从传递函数来看，巴特沃斯和切比雪夫滤波器的传输函数

8、都是一个常数除以一个多项式, 为全极点网络, 仅在无限大阻带处衰减为无限大, 而椭圆函数滤波器在有限频率上既有零点又有极点。极零点在通带内产生等波纹, 阻带内的有限传输零点减少了过渡区, 可获得极为陡峭的衰减曲线。也就是说在阶数相同的条件下,椭圆滤波器相比于其他类型的滤波器，能获得更窄的过渡带宽和较小的阻带波动, 就这点而言, 椭圆滤波器是最优的。它陡峭的过渡带特性是用通带和阻带的起伏为代价来换取的，并且在通带和阻带的波动相同，这一点区别于在通带和阻带都平坦的巴特沃斯滤波器，以及通带平坦、阻带等波纹或是阻带平坦、通带等波纹的切比雪夫滤波器。 总结起来，椭圆滤波器具有以下特点：1、 椭圆低通滤波

9、器是一种零、极点型滤波器，它在有限频率范围内存在传输零点和极点。2、 椭圆低通滤波器的通带和阻带都具有等波纹特性，因此通带，阻带逼近特性良好。 3、 对于同样的性能要求，它比前两种滤波器所需用的阶数都低，而且它的过渡带比较窄。但是椭圆滤波器传输函数是一种较复杂的逼近函数, 利用传统的设计方法进行电路网络综合要进行繁琐的计算, 还要根据计算结果进行查表, 整个设计, 调整都十分困难和繁琐。而用MATLAB设计椭圆滤波器可以大大简化设计过程。2.椭圆滤波器设计的数学推导椭圆滤波器的振幅平方函数为 ： （1）其中是雅可比(Jacobi) 椭圆函数，雅可比椭圆函数是阶数N的有理函数，N=5时的特性曲线

10、如图1所示。图1 N=5时雅可比椭圆函数的特性曲线 由图1 可见,在归一化通带内( - 1 1) ,() 在(0 ,1) 间振荡,而超过 后,() 在( , ) 间振荡。L 越大，也变大。这一特点使滤波器同时在通带和阻带具有任意衰减量。L 是一个表示波纹性质的参量。雅可比椭圆函数还具有以下性质： （2）阶数N等于通带和阻带内最大点和最小点的总和，为与通带衰减有关的参数。系统函数和阶数N是由系统下面的性能指标来确定的，主要有：截止频率,通带内最大衰减和阻带截止频率以及阻带内最小衰减。假定是频率归一化的基准频率，即 （3）定义频率的选择性因数k为 （4）则截止频率分别归一化为 （5）再次假定 （6

11、） （7） （8） （9）则得到椭圆滤波器的阶数N为 （10）这时，令归一化的基准频率为，则得到归一化后的椭圆低通滤波器的系统函数为 （11）式中， 所以，实际的椭圆低通滤波器就可以由归一化的系统函数来得到 （12） 图2 为典型N 为奇数的椭圆滤波器的幅度特性, 当，和A 确定后,阶次N 即可确定,进而可以设计出椭圆滤波器。 图2-1 椭圆滤波器的幅度特性3.关于归一化的讨论归一化是一种简化计算的方式，主要是为了数据处理方便提出来的，即将有量纲的表达式，经过变换，化为无量纲的表达式，成为纯量。比如，复数阻抗可以归一化书写：Z = R + jL = R(1 + jL/R) ，复数部分变成了纯数

12、量了，没有量纲。归一化方法（Normalization Method）把数据映射到01范围之内处理，更加便捷快速，应该归到数字信号处理范畴之内。其具体作用是归纳统一样本的统计分布性。归一化在01之间是统计的概率分布，归一化在-1+1之间是统计的坐标分布。归一化有同一、统一和合一的意思。无论是为了建模还是为了计算，其基本度量单位要同一。但是归一化处理并不总是合适的，根据输出值的分布情况，标准化等其它统计变换方法有时可能更好。具体情况还需具体分析。2.3 双线性变换法脉冲响应不变法的主要缺点是频谱交叠产生的混淆，这是从S平面到Z平面的标准变换zesT的多值对应关系导致的,为了克服这一缺点，设想变换

13、分为两步：第一步：将整个S平面压缩到S1平面的一条横带里； 第二步：通过标准变换关系将此横带变换到整个Z平面上去。由此建立S平面与Z平面一一对应的单值关系，如图2-1所示。图 2-2 双线性变换法中的复平面映射为了将s平面的j轴压缩到s1平面j轴上的- 一段上，可通过以下的正切变换实现： 这里C是待定常数，下面会讲到用不同的方法确定C，可使模拟滤波器的频率特性与数字滤波器的频率特性在不同频率点有对应关系。 经过这样的频率变换，当1由 时,由 ,即映射了整个j轴。 将这一关系解析延拓至整个s平面，则得到s平面 平面的映射关系：再将s1平面通过标准变换关系映射到z平面，即令 ，通常取C=2/T ，

14、最后得S平面与Z平面的单值映射关系：双线性换法的主要优点是S平面与Z平面一一单值对应，S平面的虚轴(整个j)对应于Z平面单位圆的一周，S平面的=0处对应于Z平面的=0处， 对应 即数字滤波器的频率响应终 止于折迭频率处，所以双线性变换不存在混迭效应。3设计步骤3.1 设计流程图音乐信号滤波去噪使用双线性变化法设计的并联型椭圆滤波器的设计流程如图3-1所示。在windows平台下用麦克风录制一段语音，并且加入音乐背景加入单频噪声对语音信号进行频谱分析，画出时域和频域波形图 使用双线性变换法设计并联型椭圆滤波器画出其频率响应用设计好的椭圆滤波器对语音信号进行滤波画出语音信号滤波前后波形并且进行比较

15、分析开始结束用wavread函数对该语音信号进行采集 图3-1 设计流程图3.2 拉弦音乐的采集先在windows平台下用麦克风录制一段语音，然后加入一段音乐背景。如图3-2所示。图3-2 录音机使用图然后在Matlab软件平台下，利用函数wavread对语音信号进行采样，记住采样频率和采样点数。Wavread函数调用格式：x,fs,nbits=wavread(file);功能说明：采样值放在向量x中，fs表示采样频率（Hz），bits表示采样位数。采集完成后在信号中加入一频率为2200Hz的单频噪声，设计的任务即为从含噪信号中滤除单频噪声，还原原始信号。x,fs,bits=wavread(&

16、#39;e:201054080226.wav'); % 输入参数为文件的全路径和文件名， %输出的第一个参数是每个样本的值， %fs是生成该波形文件时的采样率， %bits是波形文件每样本的编码位数。sound(x,fs,bits); % 按指定的采样率和每样本编码位数回放N=length(x); % 计算信号x的长度fn=2300; % 单频噪声频率，此参数可改t=0:1/fs:(N-1)/fs; % 计算时间范围，样本数除以采样频率xx=x(:,1)' y=x+0.1*sin(fn*2*pi*t); sound(y,fs,bits); %可以明显听出有尖锐的单频啸叫声3.3

17、语音信号的频谱分析在MATLAB中编辑m函数，使用wavread函数读取采集的声音文件（.wav)将它赋值给某一向量，再对其进行采样后使用plot语句画出相关的频谱图形。快速傅里叶变换算法FFT计算DFT的函数fft,其调用格式如下：Xk=fft(x,n)参数x为被变换的时域序列向量，N是DFT变换区间长度，当n大于x的长度时，fft函数自动在x后面补零。当n小于xn的长度时，fft函数计算x的前n个元素，忽略其后面的元素。在本次课程设计中，我们利用fft函数对单音频信号进行快速傅里叶变换，就可以得到信号的频谱特性。X=abs(fft(x);Y=abs(fft(y); % 对原始信号和加噪信号

18、进行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('时间(s)');ylabel('幅度'); title('原始音乐信号');axis(0,3,-0.5,0.5);subplot(2,2,2);plot(f,X);xlabel('频率(Hz)');ylabel('幅度谱');title('原

19、始音乐信号幅度谱');axis(0,3000,0,1000); subplot(2,2,3);plot(t,y);xlabel('时间(s)');ylabel('幅度');title('加干扰后的音乐信号');axis(0,3,-0.5,0.5); subplot(2,2,4);plot(f,Y);xlabel('频率(Hz)');ylabel('幅度谱');title('加干扰后的音乐信号幅度谱');axis(0,3000,0,1000);得到的原始单音频信号和加上噪音后的单音频信号的时域

20、波形和频谱图如图3-4所示。图3-4 原始单音频信号和加噪后信号的时域图和频谱图3.4滤波器设计设计指标：通带截止频率为2000Hz，阻带截止频率为2150Hz，通带波纹为1dB，阻带波纹为45dB，用双线性变换法设计的一个满足上述指标的椭圆IIR滤波器。fp=fn-200;fc=fn-50; %定义通带和阻带截止频率Rp=1;As=45; % 定义通带波纹和阻带衰减wp=fp/fs*2*pi;ws=fc/fs*2*pi; %计算对应的数字频率T=1; %定义采样间隔OmegaP=(2/T)*tan(wp/2);OmegaS=(2/T)*tan(ws/2); %截止频率预畸变c,d=afd_e

21、lip(OmegaP,OmegaS,Rp,As); % 计算滤波器系统函数分子分母系数b,a=bilinear(c,d,T); % 双线性变换得到数字滤波器系统函数分子分母系数db,mag,pha,grd,w=freqz_m(b,a); %计算椭圆低通滤波器的幅度，幅度响应，群延时和相位响应等数据图delta=1,zeros(1,99);ha=filter(b,a,delta); % 计算脉冲响应figure;subplot(2,2,1);plot(w/pi,db);grid;axis(0,1,-100,0);xlabel('w/pi');ylabel('dB'

22、); title('滤波器幅度响应图');subplot(2,2,2);plot(w/pi,mag);grid;axis(0,1,0,1);xlabel('w/pi');ylabel('幅度mag'); title('滤波器幅度响应图');subplot(2,2,3);plot(w/pi,pha);grid;axis(0,1,-4,4);xlabel('w/pi');ylabel('相位pha'); title('滤波器相位响应图');subplot(2,2,4);plot(ha);

23、grid;axis(0,100,-0.05,0.1);xlabel('n');ylabel('h(n)'); title('滤波器脉冲响应图');得到椭圆滤波器的幅度和相位谱如图3-5所示。图3-5 椭圆滤波器的幅度和相位谱3.5画出滤波器并联型结构图利用C,B,A=dir2par(b,a)将滤波器转换成级联型,得出C = -0.2948B = -0.0074 0.0280 0.1225 -0.0921 -0.6592 0.1720 0.9031 0A = 1.0000 0.1540 0.9715 1.0000 0.0077 0.8606 1.0

24、000 -0.4534 0.55091.0000 -0.4727 0根据以上参数画出的并联型结构图为图3-6 滤波器并联型结构图3.6滤波和结果分析滤波程序：y_fil=filter(b,a,y); % IIR滤波器对信号x1进行滤波处理，其中b,a为上面设计好的滤波器参数。Y_fil=abs(fft(y_fil); % 对fft变换,取幅度谱Y_fil=Y_fil(1:N/2); % 截取前半部分 deltaf=fs/N; % 计算频谱的谱线间隔 figure;subplot(3,2,1);plot(t,x);xlabel('时间t');ylabel('幅度'

25、); title('原始音乐信号');axis(0,3,-1,1);grid on;subplot(3,2,2);plot(f,X);xlabel('频率f');ylabel('幅度');title('原始音乐信号幅度谱');axis(0,3000,0,1000);grid on;subplot(3,2,3);plot(t,y);xlabel('时间t');ylabel('幅度');title('加干扰后的音乐信号');axis(0,3,-1,1);grid on;subplot(3

26、,2,4);plot(f,Y);xlabel('频率f');ylabel('幅度');title('加干扰后的音乐信号幅度谱');axis(0,3000,0,1000);grid on;subplot(3,2,5);plot(t,y_fil);xlabel('时间t');ylabel('幅度'); title('滤波后语音信号时间y');axis(0,3,-1,1);grid on;subplot(3,2,6);plot(f,Y_fil);xlabel('频率f');ylabel(&

27、#39;幅度');title('滤波后语音信号幅度谱Y');axis(0,3000,0,1000);grid on;比较滤波前后信号的波形及频谱图如图3-6所示。图3-7 滤波前后信号的时域和频谱图3.6滤波前后语音的对比在MATLAB中，经sound函数，对经过椭圆滤波器之后的信号进行回放，可以听出滤波之后的信号比原始信号更清晰一些，清除了环境噪音。通过以下语句来进行语音信号回放比较：sound（x，fs，bits） %播放原始双音频信号sound（y，fs，bits）%播放加噪后的双音频信号sound（y_fil，fs，bits） %播放经过滤波处理后的双音频信号所

28、得结果证明了椭圆滤波器去噪设计成功。4 出现的问题及解决方法 虽然已不是第一进行课程设计，但由于理论知识没学好还是遇到不少问题。问题总结如下：（1）在windows平台下用麦克风录制一段语音后不会加入音乐背景。解决办法：请教会的同学，在同学的指导下利用软件完成了背景音乐的加入。（2）问题：在用wavread函数加入语音信号时，提示错误。解决办法：输出语音文件的正确路径。（3）问题：加噪声前后幅度和幅度谱坐标不一样，不好做对比分析。解决办法：利用axis对坐标图做调整，使得前后坐标一样，便于对比分析。5 课程设计心得体会在这次数字信号课程设计中，虽然时间不长只有二周，但让我学到了很多东西。首先我

29、了解到了MATLAB功能的强大性与掌握编程各种函数和语句的重要性，再次体会到平时学习与实践之间的差距，让我意识到不仅要学好课本上的基本知识，而且要加强动手实践能力。这次我的课程设计内容是音乐滤波去噪使用双线性变换法设计的椭圆并联型滤波器，采用MATLAB设计滤波器，使原来非常繁琐复杂的程序设计变成了简单的函数调用，为滤波器的设和实现开辟了广阔的天地，尤其是MATLAB工具箱使各个领域的研究人员可以直观方便地进行科学研究与工程应用。其中的信号处理工具箱、图像处理工具箱、小波工具箱等更是为数字滤波研究的蓬勃发展提供了可能。两周的课程设计让我对MATLAB软件的使用更加的熟练，对椭圆滤波器的滤波原理

30、有了更深刻的认识。在二周的课程设计时间里，我每天都过的很充实，查资料、读程序，反复琢磨，学习的过程是艰辛的，但是同时也是快乐的。通过实实在在的课程设计，发现自己在课堂上所学的知识对于解决实际问题来说，是远远不够的，要想掌握技术，还需要加倍的努力。不过，在这短短的时间里，还是让我对数字信号处理有了更深了解，也学到了新的知识让我掌握了MATLAB等软件的简单使用，明白了软件仿真对设计的重要性。再次，这次课程设计让我充分认识到团队合作的重要性，只有分工协作才能保证整个项目的有条不紊。在整个课程设计过程中我再次认识到坚持、耐心、细心等品质的重要，这对今后的学习和工作是有很大帮助的。总之，这次课程设计让

31、我获益良多。最后在此我由衷地感谢各位帮助我的老师和同学！ 6 参考文献1 程佩青数字信号处理教程M北京：清华大学出版社，20022 薛年喜主编MATLAB 在数字信号处理中的应用M北京：清华大学出版社，20023 维纳K恩格尔，约翰G普罗克斯.数字信号处理M.西安交通大学出版社，20024 董长虹等. MATLAB信号处理与应用M.北京：国防工业出版社，20055 美 M.H.海因斯 著，张建华等 译.数字信号处理M.北京：科学出版社，20026 张葛祥，李 娜. MATLAB仿真技术与应用M.北京：清华大学出版社，20077 数字椭圆滤波器的Matlab设计与实现 王靖 李永全 现代电子技术

32、2007年第6 期总第245附录1：语音信号滤波去噪设计源程序清单% 程序名称：201054080226.m% 程序功能：使用双线性变换法设计的并联型椭圆滤波器对音乐信号滤波去噪。% 程序作者：马子森/% 最后修改日期：2013-3-7x,fs,bits=wavread('e:201054080226.wav'); % 输入参数为文件的全路径和文件名， %输出的第一个参数是每个样本的值， %fs是生成该波形文件时的采样率， %bits是波形文件每样本的编码位数。sound(x,fs,bits); % 按指定的采样率和每样本编码位数回放N=length(x); % 计算信号x的长

33、度fn=2300; % 单频噪声频率，此参数可改t=0:1/fs:(N-1)/fs; % 计算时间范围，样本数除以采样频率xx=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);p

34、lot(t,x);xlabel('时间(s)');ylabel('幅度'); title('原始音乐信号');axis(0,3,-0.5,0.5);subplot(2,2,2);plot(f,X);xlabel('频率(Hz)');ylabel('幅度谱');title('原始音乐信号幅度谱');axis(0,3000,0,1000); subplot(2,2,3);plot(t,y);xlabel('时间(s)');ylabel('幅度');title('加

35、干扰后的音乐信号');axis(0,3,-0.5,0.5); subplot(2,2,4);plot(f,Y);xlabel('频率(Hz)');ylabel('幅度谱');title('加干扰后的音乐信号幅度谱');axis(0,3000,0,1000);fp=fn-200;fc=fn-50; %定义通带和阻带截止频率Rp=1;As=45; % 定义通带波纹和阻带衰减wp=fp/fs*2*pi;ws=fc/fs*2*pi; %计算对应的数字频率T=1; %定义采样间隔OmegaP=(2/T)*tan(wp/2);OmegaS=(2/T)

36、*tan(ws/2); %截止频率预畸变c,d=afd_elip(OmegaP,OmegaS,Rp,As); % 计算滤波器系统函数分子分母系数b,a=bilinear(c,d,T); % 双线性变换得到数字滤波器系统函数分子分母系数db,mag,pha,grd,w=freqz_m(b,a); %计算椭圆低通滤波器的幅度，幅度响应，群延时和相位响应等数据图delta=1,zeros(1,99);ha=filter(b,a,delta); % 计算脉冲响应figure;subplot(2,2,1);plot(w/pi,db);grid;axis(0,1,-100,0);xlabel('w/pi');ylabel('dB'); title('滤波器幅度响应图');subplot(2,2,2);plot(w/pi,mag);grid;axis(0,1,0,1);xlabel('w/p