吹管音乐滤波去噪使用双线性变换法设计切比雪夫II型滤波器.doc

邹磊吹管音乐滤波去噪使用双线性变换设计切比雪夫II型滤波器第23页 共23页 吹管音乐滤波去噪使用双线性变换法设计切比雪夫II型滤波器学生姓名：邹 磊 指导老师：胡双红摘 要 本课程设计的主要目的是设计一个切比雪夫II型滤波器，并对吹管音乐进行滤波去噪。开发平台为MATLAB7.1,设计方法为双线性变换法。用酷狗音乐剪辑一段吹管音乐，并转换为WAV格式，绘制波形并观察其频谱，给定相应技术指标，用双线性变换法设计的一个满足指标的切比雪夫II型IIR滤波器，对该语音信号进行滤波去噪处理，比较滤波前后的波形和频谱并进行分析。程序通过调试运行，初步实现了设计目标，并且经过适当完善后，已达到设计要求。关键词 课程设计；滤波去噪；IIR滤波器；双线性变换法；MATLAB7.11 引 言用酷狗音乐采集一段8000Hz，8bit的单声道吹管音乐，绘制波形并观察其频谱，给定通带截止频率为2000Hz，阻带截止频率为2150Hz，通带波纹为1dB，阻带波纹为35dB，用双线性变换法设计的一个满足上述指标的切比雪夫II型IIR滤波器，对该吹管音乐进行滤波去噪处理。1.1 课程设计目的综合运用数字信号处理的理论知识进行频谱分析和滤波器设计，通过理论推导得出相应结论，再利用 MATLAB 作为编程工具进行计算机实现，从而加深对所学知识的理解，建立概念。本次使学生加深对理论知识理解的同时增强其逻辑思维的能力，另一方面对课堂教学中的理论知识做一个总结和补充。熟悉和巩固模拟滤波器的设计方法和原理 ，掌握双线性变换方法切比雪夫滤波器设计方法 ，实现滤波器设计的有关经典算法 ，熟练掌握使用MATLAB语言设计各种要求的数字滤波器 。可综合运用这些知识解决一定的实际问题，使学生在所学知识的综合运用能力上以及分析问题、解决问题能力上得到一定的提高。同时通过课程设计培养学生严谨的科学态度，认真的工作作风和团队协作精神。1.2课程设计的要求（1）学会 MATLAB 的使用，掌握 MATLAB 的程序设计方法；（2）滤波器指标必须符合工程实际，根据模拟滤波器的性能指标，确定数字滤波器指标；（3）采用双线性变换法，设计满足上述性能指标要求的ChebyshevII型数字低通滤波器；（4）设计完后应检查其频率响应曲线是否满足指标；（5）处理结果和分析结论应该一致，而且应符合理论；（6）独立完成课程设计并按要求编写课程设计报告书；1.3设计平台MATLAB是一种既可交互使用又能解释执行的计算机变成语言，可以用直观的数学表达式来描述问题，从而避开繁琐的底层编程，并且是解决工程技术问题的计算平台。它拥有丰富的函数资源和工具箱资源，MTLAB最大的特点就是简洁开放的程序代码和直观实用的开发环境。具备库函数资源丰富；语言精练，代码灵活；面向对象，控制功能优良；程序设计优良，图形功能强大；源代码开放，形形色色工具箱等优点，但是也有运行速度慢，占内存较大的缺点。MATLAB已经成为国际上最流行的科学与工程计算软件工具，现在的MATLAB已经不仅仅是一个“矩阵实验室”了，它已经成为一种具有广泛应用前景的全新的计算机高级语言，在国内外高校和研究部门中正扮演重要角色。2 设计原理 在网站上下载一段吹管音乐，绘制波形并观察其频谱，给定相应技术指标，用双线性变换法设计的切比雪夫II 型IIR滤波器，对该吹管音乐进行滤波去噪处理，比较滤波前后的波形和频谱并进行分析。2.1 IIR滤波器IIR滤波器具有无限长脉冲响应，因此能够与模拟滤波器相匹敌；因此，IIR滤波器设计的基本方法是利用复值映射将大家熟知的模拟滤波器变换为数字滤波器。这一方法的优势在于各种模拟滤波器设计（AFD）表格和映射在文献中普遍都能获得。这个基本方法称为A/D（模拟-数字）滤波器变换。然而，AFD表格仅对低通滤波器适用，而同时要得到设计其他频率选择性滤波器（高通、带通、带阻等等）。为此，需要对低通滤波器实行频带变换。这种IIR滤波器设计的基本方法存在两种途径：设计模拟低通滤波器途径1：实行滤波器变换SZ实行频带变换 S S期望的IIR滤波器途径2：设计模拟低通滤波器实行频带变换ZZ实行滤波器变换SZ期望的IIR滤波器我们将研究途径二的设计基本方法，在这种IIR滤波器设计方法将按下列步骤进行：（1）设计模拟低通滤波器。（2）研究并实行滤波器变换以得到数字低通滤波器。（3）研究并实行频带变换以便从数字低通滤波器得到其他数字滤波器。设IIR滤波器的输入序列为x(n)，则IIR滤波器的输入序列x(n)和输出序列y(n)之间的关系可用下列方程式表示：y(n)=i=0Mbixn-i+j=1Najy(n-j) (2-1)其中bi和aj是滤波器的系数，且aj中至少有一个非零。与上述差分方程对应的传递函数为：Hz=Y(z)X(z)=b0+b1z-1+bMz-M1-a1z-1-aNz-N (2-2)由传递函数可以发现无限长单位冲激响应滤波器有如下特点：（1） 单位冲激响应h(n)是无限长的。（2）系统传递函数H(z)在有限z平面上有极点存在。（3）结构上存在着输出到输入的反馈，也就是结构上是递归型的。设Ha（j）是某个模拟滤波器的频率响应，那么低通滤波器在幅度平方响应上的技术指标给出为：11+2|Ha（j）|21，|p （2-3）0|Ha（j）|21A2，s| （2-4）式中是通带波纹参数，p是通带截止频率以rad/s（弧度每秒）计，A是阻带衰减参数，以及s是阻带截止频率以rad/s计。图2.1 模拟低通滤波器技术指标要求由图可见，|Ha（j）|2必须满足|Haj|2=11+2，=p （2-5）|Haj|2=1A2，=p （2-6）参数和A是分别与以dB计的参数Rp和As有关的，这些关系是Rp=-10log1011+2=10Rp10-1 （2-7）As=-10log101A2A=10As20 （2-8）波纹1和2的绝对标尺是通过下式与和A有关的：1-11+1=11+2=211-1 21+1=1AA=1+12 利用幅度平方响应给出的模拟滤波器要求（2-1）和（2-2）式不包含任何相位信息。现在，为了求s域的系统函数Ha（s），考虑Haj=Ha（s）|s=j 那么有|Haj|2=HajHa*j=HajHa-j=HasHa-s|s=j （2-12）或者有HasHa-s=|Haj|2|=sj （2-9）因此，幅度平方函数的零点和极点相对于j轴是以镜像对称方式分布的【2】。2.2 切比雪夫II型滤波器切比雪夫滤波器是在通带或阻带上频率响应幅度等波纹波动的滤波器。在通带波动的为“I型切比雪夫滤波器”，在阻带波动的为“II型切比雪夫滤波器”。切比雪夫滤波器在过渡带比巴特沃斯滤波器的衰减快，但频率响应的幅频特性不如后者平坦。切比雪夫滤波器和理想滤波器的频率响应曲线之间的误差最小，但是在通频带内存在幅度波动。2.3 双线性变换法双线性变换由下式定义：s=2Tz-1z+1z=1+sT21-sT2 （2-11）这里T是一个参数。也叫线性分式变换，因为当乘开之后得到T2sz+T2s-z+1=0这在每个变量上都是线性的，如果另一个固定的话；或者说在s和z上是双线性的。在（2-11）式的制约下的复平面映射如下图所示。图2.2 双线性变换中的复平面映射从上图可知：（1）在（2-11）中利用s=+j得到 z=（1+T2+jT2）（1-T2-jT2） （2-12）（2）整个左半平面映射到单位圆内，所以这是一个稳定的变换。（3）虚轴映射到单位圆是以一种一对一的方式进行的，因此，在频域不存在混叠现象。在（2-12）中代入=0，得到（幅度为1）z=1+jT21-jT2=ej作为的函数解出得到=2tan-1（T2）或=2Ttan（2） （2-13）这表明是与非线性关联的，但是不存在混叠。因此，在（2-13）式中我们说是被预畴变到的。3 设计步骤3.1 设计流程图开始结束回放吹管音乐比较滤波前后吹管音乐的波形及频谱设计好的滤波器进行滤波处理验证并进行频谱分析双线性变换法设计切比雪夫II型滤波器快速傅里叶变换，并且画频谱图吹管音乐的采集（wavread函数），画时域图设定滤波器性能指标，通带截止频率fb=2000，阻带截止频率fc=2150，通带波纹Rp=1，阻带波纹As=35图3.1 双线性变换法切比雪夫II型滤波器对吹管音乐去噪流程图3.2 吹管音乐的采集在网上下载音乐“卡农”，利用酷狗音乐播放器，剪辑吹管音乐“卡农”，让其时间在2-3 s左右。转化格式为WAV格式，并调整为单声道，采样速率8kHz，8位码的音频格式。然后在Matlab软件平台下，利用函数wavread对语音信号进行采样，函数为y,Fs,bits = wavread(lgy.wav)，记住采样频率和采样点数。图3.2 语音信号设置3.3吹管音乐的频谱分析在MATLAB中编辑m函数，使用wavread函数读取采集的声音文件（.wav)将它赋值给某一向量，再对其进行采样,然后使用plot语句画出相关的频谱图形。（1）Wavread函数调用格式：x,fs,nbits=wavread(file)功能说明：采样值放在向量x中，fs表示采样频率（Hz），bits表示采样位数。（2）快速傅里叶变换算法FFT计算DFT的函数fft,其调用格式如下：Xk=fft(x,n)参数x为被变换的时域序列向量，N是DFT变换区间长度，当n大于x的长度时，fft函数自动在x后面补零。，当n小于xn的长度时，fft函数计算x的前n个元素，忽略其后面的元素。在本次课程设计中，我们利用fft函数对拉弦音乐进行快速傅里叶变换，就可以得到信号的频谱特性。（3）声音采样文件读取的程序（文件名：yuyin.wav)；吹管音乐的提取：x,fs,bits=wavread(22.wav); % 输入参数为文件的全路径和文件名，输出的第一个参数是每个样本的值，fs是生成该波形文件时的采样率，bits是波形文件每样本的编码位数。sound(x,fs,bits); % 按指定的采样率和每样本编码位数回放N=length(x); % 计算信号x的长度fn=2200; % 单频噪声频率，此参数可改t=0:1/fs:(N-1)/fs; % 计算时间范围，样本数除以采样频率x=x(:,1); y=x+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.4所示。图3.4 原始吹管音乐和加噪后信号的时域和频谱图3.4滤波器设计 设计指标：通带截止频率为2000Hz，阻带截止频率为2150Hz，通带波纹为1dB，阻带波纹为35dB，用双线性变换法设计的一个满足上述指标的切比雪夫II型滤波器双线性变换法设计切比雪夫II型滤波器fp=fn-200;fc=fn-50; %定义通带和阻带截止频率Rp=1;As=35; % 定义通带波纹和阻带衰减wp=fp/fs*2*pi; ws=fc/fs*2*pi; %计算对应的数字频率T=1;Fs = 1/T; %定义采样间隔Omegap=(2/T)*tan(wp/2);Omegas=(2/T)*tan(ws/2); %截止频率预畸变cs,ds = afd_chb1(Omegap,Omegas,Rp,As); %选择滤波器最小阶数b,a = bilinear(cs,ds,Fs);C,B,A = dir2cas(b,a) %双线性变换法实现模拟滤波器到数字滤波器的转换db,mag,pha,grd,w=freqz_m(b,a); %绘制数字滤波器频率响应幅度图delta=1,zeros(1,99);ha=filter(b,a,delta); ha=filter(b,a,delta);得到切比雪夫滤波器的幅度和相位谱如图3.5所示。图3.5 切比雪夫滤波器的幅度和相位谱3.5 信号滤波处理滤波程序：y_fil=filter(b,a,y); % 用设计好的滤波器对y进行滤波Y_fil=abs(fft(y_fil);Y_fil=Y_fil(1:N/2); % 计算频谱取前一半画出滤波前后拉弦音乐的时域和频谱图如图3.6所示。并比较加噪与去噪的吹管音乐频谱图如图3.7所示。图3.6 滤波前后吹管音乐音乐的时域和频谱图图3.7 加噪与去噪的吹管音乐频谱图比较上面是加语音信号幅度谱图，由图可知：吹管音乐的噪声是加在高频上，而我用切比雪夫低通滤波器除去噪声，得到滤波的信号与原始的吹管音乐相差很小。3.6结果分析在MATLAB中，经sound函数，对经过切比雪夫II型滤波器之后的信号进行回放，可以听出滤波之后的信号比原始信号更清晰一些，清除了环境噪音。通过以下语句来进行语音信号回放比较：sound（x，fs，bits） %播放原始吹管音乐sound（y，fs，bits）%播放加噪后的吹管音乐sound（y_fil，fs，bits） %播放经过滤波处理后的吹管音乐所得结果证明了切比雪夫I型滤波器去噪设计成功。4出现的问题及解决方法（1）采集的语音信号需要设置成8000Hz，8位，单声道，在得到理想的波形时，通过多次尝试，和查找书籍及同学讨论，最后猜得到理想的吹管音乐的时域图和频谱图；（2）最好将角频率的频谱图转化为频率的频谱图，这样在观察通带截止频率和阻带截止频率是更加精确，误差更小；（3）由于使用的MATLAB软件版本不同，导致在自己机器上始终不能画出正确的波形，最后只能在其他机器上完成此次的波形图；（4）我们设计的滤波器的阻带截止频率，在画出的频谱图上应该更改坐标，由此观察设计的滤波器是否满足要求，若不满足则可以及时更改。（5）所有的时间波形横坐标都要化为时间，滤波前后频谱的横坐标应是频率，这样在观察通带截止频率和阻带截止频率时更加精确，误差较小。5 结束语本次数字信号处理课程设计历时两个星期，此次课程设计采用Matlab进行数字信号处理课程设计，明显可以看出，通过这次设计，让我们加深了对课堂抽象概念的理解，巩固了课堂上所学的理论知识，并能很好地理解与掌握数字信号处理中的基本概念、基本原理、基本分析方法。在本次课程设计中，让学生自己剪辑音乐，设计滤波器对声音进行处理，大大激发了我们的学习兴趣，使我们很快地掌握编程方法和解决实际问题的技巧，取得了良好的教学效果。我们在学习这门课程时,普遍感到数字信号处理的概念抽象,对其中的分析方法与基本理论不能很好地理解与掌握。因此,如何帮助学生理解与掌握课程中的基本概念、基本原理、基本分析方法以及综合应用所学知识解决实际问题的能力,是本课程教学中所要解决的关键问题。通过这次课程设计，让我对这些抽象的概念的理解转化为实实在在的运用与扩展，从中感受到不同与课本的知识运用，加强了我对本科目的理解和升华。在此我要特别感谢我们的指导老师们，还有各位同学的帮助，你们的教诲让我受益非浅，我也将更加努力进步学习和完善自己。最后向老师和同学们表示最真诚的谢意。 6 参考文献1 程佩青数字信号处理教程M北京：清华大学出版社，20022 薛年喜主编MATLAB 在数字信号处理中的应用M北京：清华大学出版社，20023 维纳K恩格尔，约翰G普罗克斯. 数字信号处理 M.西安交通大学出版社，20024 董长虹等. MATLAB信号处理与应用M.北京：国防工业出版社，20055 美 M.H.海因斯 著，张建华等 译.数字信号处理M.北京：科学出版社，20026 张葛祥，李 娜. MATLAB仿真技术与应用M.北京：清华大学出版社，2007附录1：语音信号滤波去噪设计源程序清单% 程序名称：吹管音乐滤波去噪使用双线性变换设计切比雪夫II型滤波器% 程序功能：采用双线性变换设计切比雪夫II型滤波器对含噪吹管音乐进行滤波去噪处理。% 程序作者：邹磊/% 最后修改日期：2012-3-16x,fs,bits=wavread(22.wav); % 输入参数为文件的全路径和文件名，输出的第一个参数是每个样本的值，fs是生成该波形文件时的采样率，bits是波形文件每样本的编码位数。sound(x,fs,bits); % 按指定的采样率和每样本编码位数回放N=length(x); % 计算信号x的长度fn=2200; % 单频噪声频率，此参数可改t=0:1/fs:(N-1)/fs; % 计算时间范围，样本数除以采样频率x=x(:,1); y=x+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); % 截取前半部分Warning: Integer operands are required for colon operator when used as index.Warning: Integer operands are required for colon operator when used as index.deltaf=fs/N; % 计算频谱的谱线间隔f=0:deltaf:fs/2-deltaf; % 计算频谱频率范围subplot(2,2,1);plot(t,x);xlabel(时间(单位:s);ylabel(幅度);title(原始吹管音乐);gridsubplot(2,2,2);plot(f,X);xlabel(频率(单位:Hz);ylabel(幅度谱);title(原始吹管音乐幅度谱图); axis(0,4*103,0,6000);gridsubplot(2,2,3);plot(t,y);xlabel(时间(单位:s);ylabel(幅度);title(加噪后的吹管音乐);gridsubplot(2,2,4);plot(f,Y);xlabel(频率(单位:Hz);ylabel(幅度谱);title(加噪后的吹管音乐幅度谱图); axis(0,4*103,0,6000);grid使用双线性变换法设计切比雪夫II型滤波器：fp=fn-200;fc=fn-50; %定义通带和阻带截止频率Rp=1;As=35; % 定义通带波纹和阻带衰减wp=fp/fs*2*pi; ws=fc/fs*2*pi; %计算对应的数字频率T=1;Fs = 1/T; %定义采样间隔Omegap=(2/T)*tan(wp/2);Omegas=(2/T)*tan(ws/2); %截止频率预畸变cs,ds = afd_chb2(Omegap,Omegas,Rp,As); %选择滤波器最小阶数* Chebyshev-2 Filter Order = 11b,a = bilinear(cs,ds,Fs);C,B,A = dir2cas(b,a) %双线性变换法实现模拟滤波器到数字滤波器的转换C = 0.1071B = 1.0000 1.7641 1.0000 1.0000 1.2499 1.0000 1.0000 0.7567 1.0000 1.0000 0.4193 1.0000 1.0000 0.2552 1.0000 1.0000 1.0000 0A = 1.0000 0.7936 0.2346 1.0000 0.6003 0.3535 1.0000 0.3718 0.5142 1.0000 0.1798 0.6939 1.0000 0.0702 0.89031.0000 0.4354 0db,mag,pha,grd,w=freqz_m(b,a); %绘制数字滤波器频率响应幅度图delta=1,zeros(1,99);ha=filter(b,a,delta);figure(2);Subplot(2,2,1);plot(w/pi,db);Subplot(2,2,2);plot(w/pi,mag);Subplot(2,2,3);plot(w/pi,pha);Subplot(2,2,4);plot(ha);figure;%对加噪信号进行滤波处理并画出加噪前后的信号时域和频谱：y_fil=filter(b,a,y); % 用设计好的滤波器对y进行滤波Y_fil=abs(fft(y_fil);Y_fil=Y_fil(1:N/2); % 计算频谱取前一半subplot(3,2,1);plot(t,x);xlabel(时间(单位:s);ylabel(幅度);title(原始吹管音乐);gridsubplot(3,2,2);plot(f,X);xlabel(频率(单位:Hz);ylabel(幅度谱);title(原始吹管音乐幅度谱图);gridaxis(0,4*103,0,2000)subplot(3,2,3);plot(t,y);xlabel(时间(单位:s);ylabel(幅度);title(加噪后的吹管音乐);gridsubplot(3,2,4);plot(f,Y);xlabel(频率(单位:Hz);ylabel(幅度谱);title(加噪后的吹管音乐幅度谱图);gridaxis(0,4*103,0,2000)subplot(3,2,5);plot(t,y_fil);xlabel(时间(单位:s);ylabel(幅度);title(滤波后的吹管音乐);gridaxis(0,3,-1,1)subplot(3,2,6);plot(f,Y_fil);xlabel(频率(单位:Hz);ylabel(幅度谱);title(滤波后的吹管音乐幅度谱图);gridaxis(0,4*103,0,2000)figuresubplot(211);plot(f,20*log10(abs(Y)/max(abs(X);axis tight;grid on;% 绘制加噪吹管音乐频谱图subplot(212);plot(f,20*log10(abs(Y_fil)/max(abs(X);axis tight;grid on;% 绘制去噪吹管音乐频谱图sound(y_fil,fs,bits);附录二：afd_chb2函数function b,a = afd_chb2(Wp,Ws,Rp,As);% Analog Lowpass Filter Design: Chebyshev-2% -% b,a = afd_chb2(Wp,Ws,Rp,As);% b = Numerator coefficients of Ha(s)% a = Denominator coefficients of Ha(s)% Wp = Passband edge frequency in rad/sec; Wp 0% Ws = Stopband edge frequency in rad/sec; Ws Wp 0% Rp = Pass