毕业论文设计--语音信号的数字滤波处理十一--椭圆函数(hanning窗)滤波器.doc

CHANGSHA UNIVERSITY OF SCIENCE & TECHNOLOGY数字信号处理课程设计 语音信号的数字滤波处理(十一) 题目： 椭圆函数（hanning窗）滤波器 学生姓名： 学 号： 班 级: 专 业： 指导教师： 实习起止时间: 2015年6月29日至2015年7月3日 语音信号的数字滤波处理(九) 题目 数字信号处理语音信号的数字滤波处理(十一)学生姓名： 学 号： 班 级： 所在院(系): 指导教师： 语音信号的数字滤波处理(十一) 摘要 本次设计的内容为切比雪夫及hanning低通、高通、带通滤波器，并利用MATLAB平台进行设计。首先通过声音处理语句得到声音信号的时域数据，利用FFT变换可得到频域数据，以此进行频率分析。然后对原语音信号进行加噪处理，得到被污染信号。最后将被污染信号通过设计的滤波器，实现滤波功能，得到滤波后的语音信号。滤波器分别用切比雪夫II型和hanning窗设计，间接法设计IIR滤波器采用双线性变换法，滤波器设计指标由频谱分析得到。通过声音播放语句进行语音播放，可观察声音的变化；通过图形处理语句和FFT得到时域图和频谱图，可分析得到滤波器对频率的滤波功能。关键词：切比雪夫；声音处理；hanning；MATLAB；FFT目 录1 绪论.21.1 课题背景.21.2 课题目的.22 课程设计预习与原理.32.1 课程设计预习.32.1.1 卷积运算的演示.32.2.2 采样定理的演示.92.2 课程设计原理.112.2.1 频谱分析原理.11 2.2.2 IIR设计原理.12 2.2.3 FIR设计原理.123 课程设计步骤和过程.153.1 IIR设计步骤和过程.153.2 FIR设计步骤和过程.154 设计程序的调试和运行结果.164.1切比雪夫低通滤波器程序的调试和运行结果.164.2切比雪夫高通滤波器程序的调试和运行结果.184.3切比雪夫带通滤波器程序的调试和运行结果.204.4 hamming滤波器程序的调试和运行结果.235 总结.25参考文献.27附录.28附录 A.28附录 B.33第 38 页 共 56 页1 绪论1.1 课题背景随着软硬件技术的发展，仪器的智能化与虚拟化已成为未来实验室及研究机构的发展方向。虚拟仪器技术的优势在于可由用户定义自己的专用仪器系统，且功能灵活，很容易构建，所以应用面极为广泛。基于计算机软硬件平台的虚拟仪器可代替传统的测量仪器，如示波器、逻辑分析仪、信号发生器、频谱分析仪等。从发展史看，电子测量仪器经历了由模拟仪器、智能仪器到虚拟仪器，由于计算机性能的飞速发展，已把传统仪器远远抛到后面，并给虚拟仪器生产厂家不断带来连锅端的技术更新速率。目前已经有许多较成熟的频谱分析软件，如SpectraLAB、RSAVu、dBFA等。声卡是多媒体计算机最基本的配置硬件之一，价格便宜，使用方便。MATLAB是一个数据分析和处理功能十分强大的工程实用软件，他的数据采集工具箱为实现数据的输入和输出提供了十分方便的函数和命令。MATLAB是MathWorks公司于1982年推出的一款功能强大、易于使用的高效数值计算和可视化软件，它为进行算法开发、数据计算、信号分析与可视化提供了交互式应用开发环境，主要包括基本数学计算、编程环境(M语言)、数据可视化、GUIDE等。并附加了大量支持建模、分析、计算应用的工具箱来扩展MATLAB基本环境，用于解决特定领域的工程问题，如今MATLAB已广泛应用于通信、信号处理、生物医学、自动控制等领域。数字滤波器(Digital Filter, DF)在信号的过滤、检测和参数估计等方面起着重要的作用。信号往往夹杂着噪声及无用信号成分，必须将这些干扰成分滤除。数字滤波器对信号进行筛选，可通过特定频段的信号。一般来说，噪声信号往往是高频信号，而经典滤波器正是假定有用信号与噪声信号具有不同的频段，所以利用经典滤波器可以去除噪声。但如果有用信号和无用信号，或有用信号和噪声的频谱相互重叠，那么经典滤波器则不能实现理想的滤波性能。现代滤波器的作用是从含有噪声的信号中估计出信号的某些特征或信号本身，一旦信号被估计出，那么估计出来的信号与原信号相比会有更高的信噪比。这类滤波器主要有维纳滤波器、卡尔曼滤波器、线性预测滤波器以及自适应滤波器等。1.2 课题目的综合运用数字信号处理的理论知识进行频谱分析和滤波器设计，通过理论推导得出相应结论，再利用 MATLAB 作为编程工具进行计算机实现，从而加深对所学知识的理解，建立概念。2.课程设计预习与原理2.1 课程设计预习2.1.1 卷积运算的演示(1) 线性卷积序列x1(n)=2 0 1 2 5 7 0 5 0 1 1 2，序列x2(n)=2 0 1 2 5 7 0 5。动态演示两个序列进行线性卷积x1(n)x2(n)的翻转、移位、乘积、求和的过程。其中翻转采用fliplr1，程序如下：n=-7:18; M=17;yn=zeros(1,19);figure(1)stem(yn);xlabel(n)ylabel(y(n)xn1=2 0 1 2 5 7 0 5 0 1 1 2;xm1=zeros(1,7) xn1 zeros(1,7);figure(2)stem(n,xm1)xlabel(m)ylabel(x1(m)xn2=2 0 1 2 5 7 0 5;xm2=fliplr(xn2) zeros(1,18);figure(3)stem(n,xm2)xlabel(m)ylabel(x2(N-m)title(n=0)yn(1)=sum(xm1.*xm2);figure(4)stem(yn)xlabel(n)ylabel(y(n)title(n=N)for N=1:17xm3=zeros(1,N) fliplr(xn2) zeros(1,M);figure(5)stem(n,xm3)xlabel(m)ylabel(x2(N-m)title(n=N)M=M-1;yn(N+1)=sum(xm1.*xm3);figure(6)stem(yn)xlabel(n)ylabel(y(n)title(n=N)endxm3=zeros(1,18) fliplr(xn2)figure()stem(xm3);xlabel(m)ylabel(x2(N-m)title(n=N)yn(19)=sum(xm1.*xm3);figure()stem(yn)xlabel(n)ylabel(y(n)title(n=N)线性卷积动态演示过程见：附录A。(2) 循环卷积序列x1(n)=2 0 0 7 5 7 1 7 0 1 1 9，序列x2(n)=2 0 0 7 4 7 0 3，N=12。动态演示两个序列进行圆周卷积2, 3x1(n)x2(n)的翻转、移位、乘积、求和的过程。程序如下：n=0:11; yn=zeros(1,12);figure(1)stem(yn);xlabel(n)ylabel(y(n)xn1=2 0 0 7 5 7 1 7 0 1 1 9;figure(2)stem(n,xn1)xlabel(m)ylabel(x1(m)xn2=2 0 0 7 4 7 0 3;xm2=xn2 zeros(1,length(xn1)-length(xn2);figure(3)stem(n,xm2)xlabel(m)ylabel(x2(m)title(n=0);yn(1)=sum(xn1.*xm2);figure(4)stem(yn)xlabel(n)ylabel(y(n)title(n=N)for N=1:11xm1=fliplr(xn1(1:N) fliplr(xn1(N+1:12);figure()stem(n,xm1)xlabel(m)ylabel(x1(N-m)title(n=N)yn(N)=sum(xm1.*xm2);figure()stem(n,yn)xlabel(n)ylabel(y(n)title(n=N)endfigure()xm1=fliplr(xn1);stem(n,xm1)xlabel(m)ylabel(x1(N-m)title(n=N)yn(12)=sum(xm1.*xm2);figure()stem(n,yn)xlabel(n)ylabel(y(n)title(n=N)当N12时，圆周卷积发生混叠；当N20时，圆周卷积等于线性卷积。圆周卷积动态演示过程见：附录A。(3) 声音文件线性卷积序列x1(n)=2 0 1 2 5 7 0 5,读取一段声音数据（Windows XP 菜单命令），当循环卷积长度大于或等于两序列长度之和时，循环卷积等于线性卷积。因为直接用FFT进行1024点卷积大于两序列长度，所以可用线性卷积替代圆周卷积，其程序如下：n=-7:686; M=685;yn=zeros(1,687);xn1,fs,nbits=wavread(G:Windows XP 菜单命令);xn1(:,2)=xm1=zeros(1,7) xn1 zeros(1,7);xn2=2 0 1 2 5 7 0 5;xm2=fliplr(xn2) zeros(1,686);xm1=xm1(:,1:694)yn(1)=sum(xm1.*xm2);for N=1:685xm3=zeros(1,N) fliplr(xn2) zeros(1,M);M=M-1;yn(N+1)=sum(xm1.*xm3)endxm3=zeros(1,686) fliplr(xn2);yn(687)=sum(xm1.*xm3);figure(1)stem(yn)xlabel(m)ylabel(y(n)title(n=N)线性卷积结果如图2.1所示。 图2.1 线性卷积结果(4) FFT 由循环卷积定理可知：对于时域序列循环卷积，可先进行FFT变换，然后频率相乘，最后对结果进行IFFT变换，即可得到时域循环卷积结果。其程序如下：程序：y,fs,nbits=wavread(G:Windows XP 菜单命令);Y=fft(y,1024);xn2=2 0 1 2 5 7 0 5;X=fft(xn2,1024);X1=rot90(X,3);Z=X1.*Y;z=ifft(Z,1024)figure(2)stem(z)axis(0,700,-1,0.6);FFT求得结果如图2.2所示。 图2.2 1024点FFT求得结果由图2.1和图2.2可得：直接用卷积求得结果与间接用FFT求得结果相同。同时可观察到间接用FFT求得速度要比直接用卷积求得速率要快，效率要高。2.1.2采样定理的演示 编写程序演示采样定理（时域采样、频谱周期延拓），同时演示采样频率小于2fc时，产生的混叠效应： 对下面连续信号进行采样：，A为幅度因子，a为衰减因子，为模拟角频率，其中n为学号(例如，王墨同学n=23)程序：fs=25;n=12;t=0:1/fs:1;t1=0:0.0005:1;A=n;a=sqrt(2)*n*pi;w=a;x=A*exp(-a*t).*sin(w*t);x1=A*exp(-a*t1).*sin(w*t1);subplot(5,1,1)plot(x1);xn=ones(1,length(t);subplot(512)stem(xn)yn=xn.*x;subplot(5,1,3)stem(yn);y=fft(yn,fs);f=fs*(0:fs/2)/fs;subplot(5,1,4)plot(f,y(1:fs/2+1);yn=ifft(y);subplot(5,1,5)plot(yn)其结果如果2.3.1所示；图2.3.1 采样定理演示1 要求输入采样频率fs（根据程序处理需要指定范围）后，在时域演示信号波形、采样脉冲及采样后信号；在频域演示不同采样频率下对应信号的频谱。程序：fs=60n=12t=0:1/fs:1;t1=0:0.0005:1;A=n;a=sqrt(2)*n*pi;w=a;x=A*exp(-a*t).*sin(w*t);x1=A*exp(-a*t1).*sin(w*t1);subplot(5,1,1)plot(x1);xn=ones(1,length(t);subplot(512)stem(xn)yn=xn.*x;stem(yn);y=fft(yn,fs);f=fs*(0:fs/2)/fs;subplot(5,1,4)plot(f,y(1:fs/2+1);yn=ifft(y);subplot(5,1,5)plot(yn)其结果如果2.3.2所示；. 图2.3.2 采样定理演示22.2 课程设计原理2.2.1 频谱分析原理时域分析只能反映信号的幅值随时间的变化情况，除单频率分量的简单波形外，很难明确提示信号的频率组成和各频率分量大小，而频谱分析能很好的解决此问题。(1) DFT与FFT对于给定的时域信号y，可以通过Fourier变换得到频域信息Y。Y可按下式计算 （2-1）式中，N为样本容量，为采样间隔。采样信号的频谱是一个连续的频谱，不可能计算出所有的点的值，故采用离散Fourier变换(DFT)，即 （2-2）式中，。但上式的计算效率很低，因为有大量的指数(等价于三角函数)运算，故实际中多采用快速Fourier变换(FFT)。其原理即是将重复的三角函数算计的中间结果保存起来，以减少重复三角函数计算带来的时间浪费。由于三角函数计算的重复量相当大，故FFT能极大地提高运算效率。(2) 频谱图为了直观地表示信号的频率特性，工程上常常将Fourier变换的结果用图形的方式表示，即频谱图。以频率f为横坐标，为纵坐标，可以得到幅值谱；以频率f为横坐标，为纵坐标，可以得到相位谱；以频率f为横坐标，为纵坐标，可以得到实频谱；以频率f为横坐标，为纵坐标，可以得到虚频谱。根据采样定理，只有频率不超过的信号才能被正确采集，即Fourier变换的结果中频率大于的部分是不正确的部分，故不在频谱图中显示。即横坐标。2.2.2 IIR设计原理(1) 切比雪夫滤波器原理切比雪夫滤波器的幅频特性具有等波纹特性。它有两种形式：振幅特性在通带内是等波纹的，在阻带内是单调下降的切比雪夫I型；振幅特性在阻带内是等波纹的，在通带内是单调下降的切比雪夫II型。(2) 双线性变换法工作原理 双线性变换中数字域频率和模拟频率之间的非线性关系限制了它的应用范围，只有当非线性失真是允许的或能被忽略时，才能采用双线性变换法，通常低通、高通、带通和带阻等滤波器等具有分段恒定的频率特性，可以采用预畸变的方法来补偿频率畸变，因此可以采用双线性变换设计方法。(3) 脉冲响应不变法工作原理 冲激响应不变法遵循的准则是使数字滤波器的单位取样响应与参照的模拟滤波器的脉冲响应的取样值完全一样，即h(n)=ha(nT),其中T为取样周期。实际是由模拟滤波器转换成为数字滤波器，就是要建立模拟系统函数Ha（S）与数字系统函数H（z）之间的关系。脉冲响应不变法是从S平面映射到z平面，这种映射不是简单的代数映射，而是S平面的每一条宽为2/T的横带重复地映射到整个z平面。 2.2.3 FIR设计原理由于IIR数字滤波器能够保留一些模拟滤波器的优良特性，因此应用很广。但是这些特性是以牺牲线性相位频率特性为代价的，即用Butterworth、切比雪夫和椭圆法设计的数字滤波器逼近理想的滤波器的幅度频率特性，得到的滤波器往往是非线性的。在许多电子系统中，对幅度频率特性和线性相位特性都有较高的要求，所以IIR滤波器在这些系统中往往难以胜任。有限长单位冲激响应(FIR)数字滤波器具有以下优良的特点：(1) 可在设计任意幅度频率特性滤波器的同时，保证精确、严格的线性相位特性。(2) FIR数字滤波器的单位冲激响应h(n)是有限长的，可以用一个固定的系统来实现，因而FIR数字滤波器可以做成因果稳定系统。 (3) 允许设计多通带(多阻带)系统。 窗函数法就是设计FIR数字滤波器的最简单的方法。它在设计FIR数字滤波器中有很重要的作用，正确地选择窗函数可以提高设计数字滤波器的性能，或者在满足设计要求的情况下，减小FIR数字滤波器的阶次。常用的窗函数4有以下几种：矩形窗(Rectangular window)、三角窗(Triangular window)、汉宁窗(Hanning window)、海明窗(Hamming window)、布拉克曼窗(Blackman window)、切比雪夫窗(Chebyshev window)、巴特里特窗(Bartlett window)及凯塞窗(Kaiser window)。在MATLAB中，实现矩形窗的函数为boxcar和rectwin，其调用格式如下： w=boxcar(N)w=rectwin(N)其中N是窗函数的长度，返回值w是一个N阶的向量，它的元素由窗函数的值组成。实际上，w=boxcar(N)等价于w=ones(N,1)。在MATLAB中，实现三角窗的函数为triang，调用格式为： w=triang(N) 在MATLAB中，实现汉宁窗的函数为hann，调用格式如下： w=hann(N)w=hann(N,sflag)Hann函数中的参数sflag为采样方式，其值可取symmetric(默认值)或periodic。当sflagsymmetric时，为对称采样；当sflagperiodic时，为周期采样，此时hann函数计算N+1个点的窗，但是仅返回前N个点。在MATLAB中，实现海明窗的函数为hamming，调用格式分别如下： w=hamming (N)w=hamming (N,sflag)其中sflag的用法同上。在MATLAB中，实现布拉克曼窗的函数为blackman，调用格式如下： w=blackman (N)w=blackman (N,sflag)在MATLAB中，实现切比雪夫窗的函数为chebwin，调用格式为： w=chebwin (N,r)其中r 表示切比雪夫窗函数的傅里叶变换旁瓣幅度比主瓣低rdB(其默认值为100dB)，且旁瓣是等纹波的。在MATLAB中，实现巴特里特窗的函数为bartlett，调用格式为： w=bartlett (N)在MATLAB中，实现凯塞窗的函数为kaiser，调用格式为： w=kaiser (N,beta)其中beta为窗函数的参数。各种窗函数的性能比较可见表2-1，在设计FIR滤波器的过程中可以根据要求选择合适的窗函数：表2-1 各种窗函数的性能比较窗 函 数第一旁瓣相对于主瓣衰减/dB主 瓣 宽阻带最小衰减/dB矩形窗134/N21三角窗258/N25汉宁窗318/N44海明窗418/N53布拉克曼窗5712/N74凯塞窗可调可调可调切比雪夫窗可调可调可调3 课程设计步骤和过程3.1 IIR设计步骤和过程(1) 按一定的规则将给出的数字滤波器的技术指标转换为模拟低通滤波器的技术指标。(2) 根据转换后的技术指标设计模拟低通滤波器G(s)。(3) 再按一定的规则将G(s)转换成H(s)。若所设计的数字滤波器是低通的，那么上述设计工作可以结束，若所设计是高通、带通或带阻滤波器，那么还需进行以下步骤。(4) 将高通、带通或带阻数字滤波器的技术指标先转换为低通模拟滤波器的技术指标，然后按照上述步骤(2)设计出低通G(s)，再将G(s)转换为所需的H(z)。由于的频率映射关系是根据推导的，所以使j轴每隔2/Ts便映射到单位圆上一周，利用冲激响应不变法设计数字滤波器时可能会导致上述的频域混叠现象。为了克服这一问题，需要找到由s平面到z平面的另外的映射关系，这种关系应保证：(1) s平面的整个j轴仅映射为z平面单位圆上的一周；(2) 若G(s)是稳定的，由G(s)映射得到的H(z)也应该是稳定的；(3) 这种映射是可逆的，既能由G(s)得到H(z)，也能由H(z)得到G(s)；(4) 如果G(j0)=1，那么。双线性Z变换满足以上4个条件的映射关系，其变换公式为 (3-1)双线性Z变换的基本思路是：首先将整个s平面压缩到一条从/Ts到/Ts的带宽为2/Ts的横带里，然后通过标准的变换关系将横带变换成整个z平面上去，这样就得到s平面与z平面间的一一对应的单值关系。3.2 FIR课程设计步骤和过程(1)根据对阻带衰减及过渡带的指标要求，选择窗函数的类型，并估计窗口长度N；(2)构造希望逼近的频率响应函数Hd(ejw)，即Hd(ejw)=Hdg(w)e-jw（N-1）/2；(3)计算hd(n)；(4)加窗得到设计结果：h(n)=hd(n)w(n)。数字高通滤波器的设计分两种情况，即无相移和相移为/2两种情况。当无相移时，N为奇数时，所设计的FIR数字高通滤波器为型滤波器；当N为偶数时，为型滤波器。当相移为/2时，N为奇数时，所设计的FIR数字高通滤波器为型滤波器；当N为偶数时，为型滤波器。数字带通滤波器的设计也分两种情况讨论，即无相移和相移为/2两种情况。当无相移时，N为奇数时，所设计的FIR数字带通滤波器为型滤波器；当N为偶数时，为型滤波器。当相移为/2时，N为奇数时，所设计的FIR数字带通滤波器为型滤波器；当N为偶数时，为型滤波器。4 设计程序的调试和运行结果4.1 切比雪夫低通程序的调试和运行结果 如图4.1所示为原语音信号的时域图和频谱图。在MATLAB中通过wavplay(filename)读取语音信号数据。其程序如下：y,fs,nbits=wavread(G:Windows Background.wav); 图4.1 原语音信号时域图和频谱图如图4.2所示为被污染语音信号的时域图和频谱图。其加噪信号程序如下：t=0:0.000045:0.90855;x1=0.06*sin(2*pi*5500*t)+0.07*sin(2*pi*6000*t); 图4.2 被污染语音信号时域图和频谱图如图4.3所示为低通滤波器的幅频图。设计指标：wpz=4500Hz; wsz=4700Hz; rp=2; rs=100;其设计程序如下：wpz=2*4500/fs;wsz=2*4700/fs;wp=2*fs*tan(wpz*pi/2);ws=2*fs*tan(wsz*pi/2);rp=2;rs=100;N,wc=cheb2ord(wp,ws,rp,rs,s);b,a=cheby2(N,rs,wc,s);bz,az=bilinear(b,a,fs);Hk=freqz(bz,az,16384,fs); 图4.3 低通滤波器幅频图如图4.4所示为滤波后语音信号的时域图和频谱图。其滤波程序如下：Yk=Z.*Hk; 图4.4 滤波后语音信号的时域图和频谱图4.2 切比雪夫高通程序的调试和运行结果 如图4.5所示为原语音信号的时域图和频谱图。在MATLAB中通过wavplay(filename)读取语音信号数据。其程序如下：y,fs,nbits=wavread(G:Windows Background.wav); 图4.5 原语音信号的时域图和频谱图如图4.6所示为被污染语音信号的时域图和频谱图。其加噪信号程序如下：t=0:0.000045:0.90855;x1=0.1*sin(2*pi*100*t)+0.2*sin(2*pi*220*t);x=rot90(x1,3),rot90(x1,3);z=x+y;图4.6 被污染语音信号的时域图和频谱图如图4.7所示为高通滤波器的幅频图。设计指标：wpz=700Hz; wsz=500Hz; rp=10; rs=100;其设计程序如下：wpz=2*700/fs;wsz=2*500/fs;wp=2*fs*tan(wpz*pi/2);ws=2*fs*tan(wsz*pi/2);rp=10;rs=100;N,wc=cheb2ord(wp,ws,rp,rs,s);b,a=cheby2(N,rs,wc,high,s);bz,az=bilinear(b,a,fs);Hk=freqz(bz,az,16384,fs);图4.7 高通滤波器的幅频图如图4.8所示为滤波后语音信号的时域图和频谱图。其滤波程序如下：Yk=Z.*Hk;图4.8 滤波后语音信号的时域图和频谱图4.3 切比雪夫带通程序的调试和运行结果如图4.9所示为原语音信号的时域图和频谱图。在MATLAB中通过wavplay(filename)读取语音信号数据。其程序如下：y,fs,nbits=wavread(G:Windows Background.wav);图4.9 原语音信号的时域图和频谱图如图4.10所示为被污染语音信号的时域图和频谱图。其加噪信号程序如下：t=0:0.000045:0.90855;x1=0.1*sin(2*pi*100*t)+0.2*sin(2*pi*220*t)+ 0.07*sin(2*pi*6000*t);x=rot90(x1,3),rot90(x1,3);z=x+y; 图4.10 被污染语音信号的时域图和频谱图 如图4.11所示为带通滤波器的幅频图。设计指标：wpz=4500Hz; wsz=4700Hz; wpzl=700Hz; wszl=500Hz; rp=20; rs=50;其设计程序如下：wpz=2*700/fs,2*4500/fs;wsz=2*500/fs,2*4700/fs;wp=2*fs*tan(wpz*pi/2);ws=2*fs*tan(wsz*pi/2);rp=20;rs=50;N,wc=cheb2ord(wp,ws,rp,rs,s);b,a=cheby2(N,rs,wc,s);bz,az=bilinear(b,a,fs);Hk=freqz(bz,az,16384,fs);图4.11 带通滤波器的幅频图如图4.12所示为滤波后语音信号的时域图和频谱图。其滤波程序如下：Yk=Z.*Hk;图4.12 滤波后语音信号的时域图和频谱图4.4 hanning程序的调试和运行结果如图4.13所示为原语音信号的时域图和频谱图。在MAT