数字信号处理课程设计

课程设计说明书NO.1语音信号的数字滤波——IIR数字滤波器的双线性变换法设计1、课程设计的目的通过对常用数字滤波器的设计和实现，掌握数字信号处理的工作原理及设计方法；掌握利用数字滤波器对信号进行滤波的方法。并能够对设计结果加以分析。2、设计步骤2.1．设计方案论证2.1.1.滤波器简介：滤波器从广义上来说对特定的频点或频点以外的频率进行有效滤波的电路，这种电路保留输入信号中的有用信息，滤除不需要的信息，从而达到信号的检测、提取、识别等不同的目的。如果处理的信号是时域离散信号，那么相应的处理系统就称为数字滤波器，由于在实际工作中被处理的信号都是幅度量化的数字信号，因此，数字滤波器实际上是用有限精度的算法实现一个线性时不变的时域离散系统。数字滤波器是一种用来过滤时间离散信号的数字系统，通过对抽样数据进行数学处理来达到频域滤波的目的。数字滤波器也是具有一定传输选择特性的数字信号处理装置，其输入、输出均为数字信号，实质上是一个由有限精度算法实现的线性时不变离散系统。数字滤波器的种类很多，分类方法也不同，滤波器在功能上总的可分为四类，即低通(LP)、高通(HP)、带通(BP)、带阻(BS)滤波器等，每一种又有模拟滤波器和数字滤波器两种形式。如果滤波器的输人和输出都是离散时间信号，则该滤波器的冲击响应也必然是离散的，这种滤波器称之为数字滤波器。从实现方法上，由有限长冲激响应的数字滤波器被称为FIR滤波器，具有无限长冲激响应的数字滤波器称作IIR滤波器。处理和设计。因此，选择用Matlab进行课程设计。2.1.2.MATLAB简介：MATLAB是一种用于科学工程数值计算和可视化的人机交互式和基于矩阵的体系。它是一种广泛应用于工程计算及数值分析领域的新型高级语言，其功能强大、简单易学、编程效率高，通常被用来解决一些实际课题和数学模型问题.。典型的应用包括数沈阳大学课程设计说明书NO.2值计算、算法预设计与验证，以及一些特殊的矩阵计算应用,如数理统计、数字信号处理、时间序列分析、动态系统仿真等。MATLAB软件的诞生,使数字信号处理系统的分析与设计问题变得简单了,它为数字信号处理系统的设计与仿真提供了一个强有力的工具。在设计数字滤波器时,通常采用MATLAB来进行辅助设计和仿真.。Matlab同时还具有信号分析工具箱，不需具备很强的编程能力，就可以很方便地进行信号分析、2.1.3.双线性变换法简介：在设计中为了克服冲激响应不变法所存在的缺点，从而提出采用双线性变换法。采用MATLAB用来解决一些实际课题和数学模型问题，使数字信号处理系统的分析与设计问题变得简单了，它为数字信号处理系统的设计与仿真提供了一个强有力的工具。2.2.课程设计的基本原理：2.2.1.双线性变换法的基本原理：本次设计用的是双线性变换法，因此这里简单介绍一下。双线性变换法的S域与Z域间的变换关系为：Re(s)Re(z)jIm(z)jIm(s)s平面z平面单位圆图1S域与Z域间的变换关系沈阳大学课程设计说明书NO.3双线性变换法实际上是将整个频率轴上的频率范围压缩到-π/T～π/T之间，再用转换到Z平面上。也就是说，第一步先将整个S平面压缩映射到S1平面的-π/T～π/T一条横带里；第二步再通过标准变换关系z=Ts1e将此横带变换到整个Z平面上去。这样就使S平面与Z平面建立了一一对应的单值关系，消除了多值变换性，也就消除了频谱混叠现象。模拟滤波器的系统函数)(asH一般具有无限带宽，为使映射后不发生频谱混叠，就必须将)(asH的频谱从无限频率范围分布压缩至有限范围。若设新的频率变量为1，根据频谱不混叠的条件，以1为变量的频谱必须分布在||1T的频率范围内。实现频谱压缩的数学关系式为)/arctan(21CT（1）式中C为一正常数。当从经0至变化时，T1则从经0至变化，式（1）为非线性单值映射，其逆关系则为)2tan(1TC（2）为了求出满足上式的复频域映射关系，设与和1对应的复频率变量分别为s和1s，给式（2）两边同乘以j，用s置换j，用1s置换1j，有2/2/2/2/11111eeee)2th(TsTsTsTsCTsCs或1e1e11TsTsCs（3）式（3）显示，)(asH映射至1s平面后，以1s为变量的系统函数是Ts1e的有理函数，因而正好能使用标准映射关系Tsz1e映射到z平面，从而有11zzCs（4）沈阳大学课程设计说明书NO.4上式称之为双线性变换公式，它建立了s平面与z平面之间的变换关系，且能够将有双线性变换的复频率平面映射若用表示数字角频率，与式（7.6-4）对应的频率关系式为)2tan(C（5）式（5）对每一有惟一的与之对应，因而数字滤波器的响应在形状上与模拟滤波器的一致，二者极值点的数目一定相同。理)(asH转换为有理)(zH。2.3.语音的信号采集：（1）按“开始－＞程序－＞附件－＞娱乐－＞录音机”的顺序操作打开Windows系统中的录音机软件，具体操作如图：图2运行录音操作步骤沈阳大学课程设计说明书NO.5（2）录入自己的声音信号并保存成文件（语音信号的长度不得少于1秒）（3）记录以下内容：①语音信号文件保存的文件名、格式（扩展名，一般应为“.wav”）②语音信号的采样速率（建议保存语音文件时选择采样频率等于8000Hz）③语音信号的时间长度图3录音及音频保存过程沈阳大学课程设计说明书NO.62.4.语音信号的分析：（1）将上一步骤中保存下来的语音信号文件“冯男.wav”复制到计算机装有Matlab软件的磁盘中相应Matlab目录中的“work”文件夹中（2）双击桌面上Matlab软件的快捷图标，打开Matlab软件（3）在菜单栏中选择“File－＞new－＞M-File”，打开m文件编辑器。（4）将代码复制到空白处，将自己的语音信号导入Matlab工作台并保存到“work”文件夹中。2.5.语音信号的时域及频谱分析：利用MATLAB中的“wavread”命令来读入（采集）语音信号，将它赋值给某一向量。再对其进行采样，记住采样频率和采样点数。在本次设计中，利用fft对语音信号进行快速傅里叶变换，就可以得到信号的频谱特性。其程序如下：[s,fs]=wavread('dota.wav');plot(s)s1=s(6300:6300+fs-1,1);wavwrite(s1,fs,'s1.wav')plot(s1)S1=fft(s1);figure(2)plot(s1)plot(abs(S1))沈阳大学课程设计说明书NO.7程序结果如下图：图4声音信号的时域波形图5声音信号的频谱图6单边谱沈阳大学课程设计说明书NO.8如图可知，读音信号的能量主要频段分布在（HLff,）fL=150Hz，fH=1800Hz。因为语音起伏大的波峰为主要能量（即红色框框住的部分），这一部分需要保留其余部分都要滤掉，由于该频谱图是以4000HZ为中心点左右对称的图形，3dB截止频率为4000HZ的左侧第一个波峰所在的点，通带截止频率为3dB截止频率左侧第一个波峰所在的点，阻带截止频率为3dB截止频率加上3db截止频率和通带截止频率之间的距离。上图是程序程序执行结果，其中第二个图是信号的FFT结果，其横坐标的具体值是X(k)中的序号k；第三个图是确定滤波频率范围的参考图，其横坐标的具体值应当是遵循DFT定义式和频率分辨率求得的：10)()]([)(NnknNWnxnxDFTkX，当k等于0时，020jknNjkknNeeW，从数字角频率上看，对应的正好是0即直流的位置，也就是说，大家在取滤波频段时，应当将主要能量（即红色框框住的部分）保留，其余频段部分的信号滤除。至于对于主要能量具体值的确定由大家自己设定。频率分辨率：X(k)=DFT[x(n)]相当于是信号x(n)的实际频谱X(ejω)=DTFT[x(n)]的采样，而x(n)又是连续时间语音信号x(t)的采样。X(k)的每两个相邻取值之间的频率间隔大小对应到语音信号x(t)的频谱中去，其频率间隔大小正好是采样结果的长度采样速率Nffs，称频率分辨率。2..6语音信号频谱分析由语音采样可知采样频率Fs=8000Hz,N=52000由下面程序段可知语音信号的长度：[s,fs,nbits]=wavread('fn');sound(s,fs,nbits);N=length(s)沈阳大学课程设计说明书NO.9运行出后得到N=52000Nffs8000Hz/52000=0.15Hz/采样2.7.滤波器的设计2.7.1确定滤波器的参数由频谱图可知：通带截止频率fp=2600Hz×0.15Hz/采样=390Hz阻带截止频率fst=3000Hz×0.15Hz/采样=450Hz3dB截止频率fc=3400Hz×0.15Hz/采样=510Hz同时设通带最大衰减δ1=2dB阻带最小衰减δ2=40dB在本次设计中将要设计的滤波器的类型为的巴特沃思低通滤波器（LPF）。2.7.2.设计滤波器：用双线性方法设计IIR滤波器由上面参数可以得出Ωp=2πfp=780πrad/sΩst=2πfs=900πrad/sΩc=2πfc=1020πrad/s则ωp=ΩpT=780π/8000=0.10πωst=ΩsT=900π/8000=0.12πωc=ΩcT=1020π/8000=0.14π沈阳大学课程设计说明书NO.10确定滤波器的阶数为:N≥)/lg(2)]110/()110lg[(11.021.0cs=2.4所以取N=3则系统函数为Ha(s)=)3)(2)(1(3SSSSSSc经查表可得归一化的6阶巴特沃思模拟低通滤波器的系统函数为：Ha(s)=10000000.20000000.2233SSSc=10000000.20000000.2101.12393SSS用双线性设计数字滤波器可得出：由)2/(cCotcC得出C=1020π×Cot(0.14π/2)=4590π=14420所以数字滤波器的系统函数为：H（z）=Ha(s)|s=c1111ZZ=4590π1111ZZH（z）=1)114590(00000.2)114590(00000.2)114590(101.11121131193ZZZZZZ由上可得系统函数的频率响应为：H(jwe)=1)114590(00000.2)114590(00000.2)114590(101.12393jwjwjwjwjwjweeeeee沈阳大学课程设计说明书NO.11IIR数字滤波器的主要程序:[s,fs]=wavread('冯男.wav');wp=0.1*pi;ws=0.12*pi;fp=390*pi;fs=450*pi;[n11,wn11]=buttord(wp,ws,1,50,'s');%求低通滤波器的阶数和截止频率[b11,a11]=butter(n11,wn11,'s');%求S域的频率响应的参数[num11,den11]=bilinear(b11,a11,0.5);%双线性变换实现S域到Z域的变换[h,w]=freqz(num11,den11);%根据参数求出频率响应plot(w*8000*0.5/pi,abs(h));legend('用butter设计');运行结果如下：图7IIR低通滤波器沈阳大学课程设计说明书NO.122.8.对信号进行加噪：(1)产生一个长度为N=8000的噪声信号。噪声信噪比SNR1取值范围15——25dB，在本次课设取值为20dB。（2）将该噪声信号迭加到语音信号“冯男.wav”上，然后导出，保存成“冯男1.wav”文件。（3）将语音与噪声合成的信号“冯男1.wav”在Matlab中按照相应的指令编写程序，得到迭加噪声后的语音信号频谱图，加噪程序如下：S1=fft(s1);figure(2)subplot(21