数字信号处理课程设计

1、西安石油大学课 程 设 计 任 务 书题 目 数字信号处理课程设计 学院(部) 信 息 工 程 学 院 专 业 电 子 信 息 工 程 班 级 电子信息工程0801 学生姓名 任斌 学 号 25 2010.12.27 至 2010.12.10 共 1 周指导教师(签字) 吴银川 数字信号处理综合设计 一、实验目的1学会matlab的使用，掌握matlab的程序设计方法；2掌握在windows环境下语音信号采集的方法；3掌握数字信号处理的基本概念、基本理论和基本方法；4掌握matlab设计fir和iir数字滤波器的方法；5学会用matlab对信号进行分析和处理。二、实验原理参考数字信号处理教材。

2、三、主要实验仪器及材料微型计算机、matlab6.5教学版、tc编程环境。四、实验内容1语音信号的采集要求利用windows下的录音机（开始程序附件娱乐录音机，文件属性立即转换8000khz，8位，单声道）或其他软件，录制一段自己的话音，时间控制在1秒左右。然后在matlab软件平台下，利用函数wavread对语音信号进行采样，记住采样频率和采样点数。通过wavread函数的使用，要求理解采样频率、采样位数等概念。wavread函数调用格式：y=wavread(file)，读取file所规定的wav文件，返回采样值放在向量y中。y,fs,nbits=wavread(file)，采样值放在向量y

3、中，fs表示采样频率（hz），nbits表示采样位数。y=wavread(file,n)，读取前n点的采样值放在向量y中。y=wavread(file,n1,n2)，读取从n1点到n2点的采样值放在向量y中。2语音信号的频谱分析要求首先画出语音信号的时域波形；然后对语音信号进行频谱分析，在matlab中，可以利用函数fft对信号进行快速付立叶变换，得到信号的频谱特性；从而加深对频谱特性的理解。3设计数字滤波器和画出频率响应根据语音信号的特点给出有关滤波器的性能指标：1）低通滤波器性能指标，fp=1000hz，fc=1200 hz， as=100db，ap=1db；2）高通滤波器性能指标，fc=

4、2800 hz，fp=3000 hz as=100db，ap=1db；3）带通滤波器性能指标，fp1=1200 hz，fp2=3000 hz，fc1=1000 hz，fc2=3200 hz，as=100db，ap=1db。要求学生首先用窗函数法设计上面要求的三种滤波器，在matlab中，可以利用函数fir1设计fir滤波器，然后在用双线性变换法设计上面要求的三种滤波器；之后再利用函数butter和cheby1设计上面要求的三种iir滤波器。最后，利用matlab中的函数freqz画出各滤波器的频率响应。4用滤波器对信号进行滤波比较fir和iir两种滤波器的性能，然后用性能好的各滤波器分别对采集

5、的信号进行滤波，在matlab中，fir滤波器利用函数fftfilt对信号进行滤波，iir滤波器利用函数filter对信号进行滤波。5比较滤波前后语音信号的波形及频谱要求在一个窗口同时画出滤波前后的波形及频谱。6回放语音信号在matlab中，函数sound可以对声音进行回放。其调用格式：sound(x,fs,bits)；可以感觉滤波前后的声音有变化。五、实验思考1双线性变换法中和之间的关系是非线性的，在实验中你注意到这种非线性关系了吗？从哪几种数字滤波器的幅频特性曲线中可以观察到这种非线性关系？ 2能否利用公式完成脉冲响应不变法的数字滤波器设计？为什么？六、实验报告要求 1.简述实验原理及目的

6、。2.按照实验步骤及要求，比较各种情况下的滤波性能。3.总结实验所得主要结论。4.简要回答思考题。数字信号处理课程设计报告学院(部) 电子信息工程学院 专 业 电子信息工程 班 级 电子信息工程0801 学生姓名 任斌 学 号 25 语音信号的数字滤波 利用双线性变换法实现iir数字滤波器的设计一课程设计的目的通过对常用数字滤波器的设计和实现，掌握数字信号处理的工作原理及设计方法；熟悉用双线性变换法设计 iir 数字滤波器的原理与方法，掌握利用数字滤波器对信号进行滤波的方法，掌握数字滤波器的计算机仿真方法，并能够对设计结果加以分析。二课程设计原理1用窗函数法设计fir滤波器根据过渡带宽及阻带衰

7、减要求，选择窗函数的类型并估计窗口长度n（或阶数m=n-1），窗函数类型可根据最小阻带衰减as独立选择，因为窗口长度n对最小阻带衰减as没有影响，在确定窗函数类型以后，可根据过渡带宽小于给定指标确定所拟用的窗函数的窗口长度n，设待求滤波器的过渡带宽为w，它与窗口长度n近似成反比，窗函数类型确定后，其计算公式也确定了，不过这些公式是近似的，得出的窗口长度还要在计算中逐步修正，原则是在保证阻带衰减满足要求的情况下，尽量选择较小的n，在n和窗函数类型确定后，即可调用matlab中的窗函数求出窗函数wd（n）。根据待求滤波器的理想频率响应求出理想单位脉冲响应hd(n)，如果给出待求滤波器频率应为hd，

8、则理想的单位脉冲响应可以用下面的傅里叶反变换式求出：在一般情况下，hd(n)是不能用封闭公式表示的，需要采用数值方法表示；从w=0到w=2采样n点，采用离散傅里叶反变换（idft）即可求出。用窗函数wd(n)将hd(n)截断，并进行加权处理，得到如果要求线性相位特性， 则h(n)还必须满足：根据上式中的正、 负号和长度n的奇偶性又将线性相位fir滤波器分成四类。 要根据所设计的滤波特性正确选择其中一类。 例如， 要设计线性相位低通特性可选择h(n)=h(n-1-n)一类，而不能选h(n)=-h(n-1-n)一类。 验算技术指标是否满足要求，为了计算数字滤波器在频域中的特性，可调用freqz子程

9、序，如果不满足要求，可根据具体情况，调整窗函数类型或长度，直到满足要求为止。2.用双线性变换法设计iir数字滤波器脉冲响应不变法的主要缺点是产生频率响应的混叠失真。这是因为从s平面到平面是多值的映射关系所造成的。为了克服这一缺点，可以采用非线性频率压缩方法，将整个频率轴上的频率范围压缩到-/t/t之间，再用z=est转换到z平面上。也就是说，第一步先将整个s平面压缩映射到s1平面的-/t/t一条横带里；第二步再通过标准变换关系z=es1t将此横带变换到整个z平面上去。这样就使s平面与z平面建立了一一对应的单值关系，消除了多值变换性，也就消除了频谱混叠现象，映射关系如图1所示。图1双线性变换的映

10、射关系为了将s平面的整个虚轴j压缩到s1平面j1轴上的-/t到/t段上，可以通过以下的正切变换实现 （1）式中,t仍是采样间隔。当1由-/t经过0变化到/t时，由-经过0变化到+，也即映射了整个j轴。将式（1）写成将此关系解析延拓到整个s平面和s1平面，令j=s，j1=s1，则得再将s1平面通过以下标准变换关系映射到z平面z=es1t从而得到s平面和z平面的单值映射关系为： (2) (3)式（2）与式（3）是s平面与z平面之间的单值映射关系，这种变换都是两个线性函数之比，因此称为双线性变换式（1）与式（2）的双线性变换符合映射变换应满足的两点要求。首先,把z=ej，可得 (4)即s平面的虚轴映

11、射到z平面的单位圆。其次，将s=+j代入式（4），得因此由此看出，当0时，|z|0时，|z|1。也就是说，s平面的左半平面映射到z平面的单位圆内，s平面的右半平面映射到z平面的单位圆外，s平面的虚轴映射到z平面的单位圆上。因此，稳定的模拟滤波器经双线性变换后所得的数字滤波器也一定是稳定的。双线性变换法优缺点双线性变换法与脉冲响应不变法相比，其主要的优点是避免了频率响应的混叠现象。这是因为s平面与z平面是单值的一一对应关系。s平面整个j轴单值地对应于z平面单位圆一周，即频率轴是单值变换关系。这个关系如式（4）所示，重写如下：上式表明，s平面上与z平面的成非线性的正切关系，如图2所示。由图2看出，

12、在零频率附近，模拟角频率与数字频率之间的变换关系接近于线性关系；但当进一步增加时，增长得越来越慢，最后当时，终止在折叠频率=处，因而双线性变换就不会出现由于高频部分超过折叠频率而混淆到低频部分去的现象，从而消除了频率混叠现象。图2双线性变换法的频率变换关系但是双线性变换的这个特点是靠频率的严重非线性关系而得到的，如式（4）及图2所示。由于这种频率之间的非线性变换关系，就产生了新的问题。首先，一个线性相位的模拟滤波器经双线性变换后得到非线性相位的数字滤波器，不再保持原有的线性相位了；其次，这种非线性关系要求模拟滤波器的幅频响应必须是分段常数型的，即某一频率段的幅频响应近似等于某一常数（这正是一般

13、典型的低通、高通、带通、带阻型滤波器的响应特性），不然变换所产生的数字滤波器幅频响应相对于原模拟滤波器的幅频响应会有畸变，如图3所示。图3双线性变换法幅度和相位特性的非线性映射对于分段常数的滤波器，双线性变换后，仍得到幅频特性为分段常数的滤波器，但是各个分段边缘的临界频率点产生了畸变，这种频率的畸变，可以通过频率的预畸来加以校正。也就是将临界模拟频率事先加以畸变，然后经变换后正好映射到所需要的数字频率上。三：课程设计的步骤1.语音信号的采集：利用windows下的录音机（开始程序附件娱乐录音机，文件属性立即转换8000khz，8位，单声道），录制一段自己的话音“信号”， 时间控制在1秒左右，然

14、后将音频文件保存“xh.wav”（2）在matlab软件平台下，利用函数wavread对语音信号进行采样，记住采样频率和采样点数。2.语音信号的频谱分析首先画出语音信号的时域波形z1=wavread(d:laoa.wav);plot(z1);图像输出如图1 图1 语音信号时域波形 对语音信号进行频谱分析，在matlab中，可以利用函数fft对信号进行快速付立叶变换，得到信号的频谱特性z1=wavread(d:laoa.wav);y1=z1(1:8192);y1=fft(y1);n=0:8191;plot(n,y1);图像输出如图2：图2 语音信号频谱分析图3. 设计数字滤波器和对信号滤波（1）

15、窗函数设计低通滤波器 程序设计如下：clear;close allz1,fs,bits=wavread(d:laoa.wav)y1=z1(1:8192);y1=fft(y1);fp=1000;fc=1200;as=100;ap=1;fs=8000;wc=2*pi*fc/fs; wp=2*pi*fp/fs;wdel=wc-wp;beta=0.112*(as-8.7);n=ceil(as-8)/2.285/wdel);wn= kaiser(n+1,beta); ws=(wp+wc)/2/pi;b=fir1(n,ws,wn);figure(1);freqz(b,1);x=fftfilt(b,z1);

16、x=fft(x,8192);figure(2);subplot(2,2,1);plot(abs(y1);axis(0,1000,0,1.0);title(滤波前信号频谱);subplot(2,2,2);plot(abs(x);axis(0,1000,0,1.0);title(滤波后信号频谱);subplot(2,2,3);plot(z1);title(滤波前信号波形);subplot(2,2,4);plot(x);title(滤波前信号波形);sound(x,fs,bits);图形分析如下：（2）窗函数设计高通滤波器程序设计如下：clear;close allz1,fs,bits=wavrea

17、d(d:laoa.wav)y1=z1(1:8192);y1=fft(y1);fp=2800;fc=3000;as=100;ap=1;fs=8000;wc=2*pi*fc/fs; wp=2*pi*fp/fs;wdel=wc-wp;beta=0.112*(as-8.7);n=ceil(as-8)/2.285/wdel);wn= kaiser(n,beta); ws=(wp+wc)/2/pi;b=fir1(n-1,ws,high,wn);figure(1);freqz(b,1);x=fftfilt(b,z1);x=fft(x,8192);figure(2);subplot(2,2,1);plot(a

18、bs(y1);axis(0,1000,0,1.0);title(滤波前信号频谱);subplot(2,2,2);plot(abs(x);axis(0,1000,0,1.0);title(滤波后信号频谱);subplot(2,2,3);plot(z1);title(滤波前信号波形);subplot(2,2,4);plot(x);title(滤波前信号波形);sound(x,fs,bits);图形分析如下：（3）窗函数设计带通滤波器程序设计如下：clear;close allz1,fs,bits=wavread(d:laoa.wav)y1=z1(1:8192);y1=fft(y1);fp1=120

19、0 ;fp2=3000 ;fc1=1000 ;fc2=3200 ;as=100 ;ap=1 ;fs=8000 ;wp1=2*pi*fp1/fs; wc1=2*pi*fc1/fs; wp2=2*pi*fp2/fs; wc2=2*pi*fc2/fs;wdel=wp1-wc1;beta=0.112*(as-8.7);n=ceil(as-8)/2.285/wdel);ws =(wp1+wc1)/2/pi,(wp2+wc2)/2/pi;wn= kaiser(n+1,beta); b=fir1(n,ws,wn);figure(1);freqz(b,1)x=fftfilt(b,z1);x=fft(x,819

20、2);figure(2);subplot(2,2,1);plot(abs(y1);axis(0,1000,0,1.0);title(滤波前信号频谱);subplot(2,2,2);plot(abs(x);axis(0,2000,0,0.0003);title(滤波后信号频谱)subplot(2,2,3);plot(z1);title(滤波前信号波形);subplot(2,2,4);plot(x);title(滤波前信号波形);sound(x,fs,bits);图形分析如下：（4）双线性变换法设计低通滤波器选用butter程序设计如下：clear;close allz1,fs,bits=wavr

21、ead(d:laoa.wav)y1=z1(1:8192);y1=fft(y1);fp=1000;fc=1200;as=100;ap=1;fs=8000;wc=2*fc/fs; wp=2*fp/fs; n,ws=buttord(wc,wp,ap,as);b,a=butter(n,ws);figure(1);freqz(b,a,512,fs);x=filter(b,a,z1);x=fft(x,8192);figure(2);subplot(2,2,1);plot(abs(y1);axis(0,1000,0,1.0);title(滤波前信号频谱);subplot(2,2,2);plot(abs(x)

22、;axis(0,1000,0,1.0);title(滤波后信号频谱);subplot(2,2,3);plot(z1);title(滤波前信号波形);subplot(2,2,4);plot(x);title(滤波前信号波形);sound(x,fs,bits);图形分析选用cheby1程序设计如下：clear;close allz1,fs,bits=wavread(d:laoa.wav)y1=z1(1:8192);y1=fft(y1);fp=1000;fc=1200;as=100;ap=1; ;fs=8000;wc=2*fc/fs;wb=2*fp/fs;n,wp=cheb1ord(wc,wb,ap

23、,as);b,a=cheby1(n,ap,wp);figure(1);freqz(b,a);x=filter(b,a,z1);x=fft(x,8192);figure(2);subplot(2,2,1);plot(abs(y1);axis(0,1000,0,1.0);title(滤波前信号频谱);subplot(2,2,2);plot(abs(x);axis(0,4000,0,0.03);title(滤波后信号频谱);subplot(2,2,3);plot(z1);title(滤波前信号波形);subplot(2,2,4);plot(x);title(滤波前信号波形);sound(x,fs,b

24、its);图形分析如下：（5），双线性变换法设计高通滤波器选用butter程序设计如下：clear;close allz1,fs,bits=wavread(d:laoa.wav)y1=z1(1:8192);y1=fft(y1);fc=2800 ;fp=3000 ;as=100;ap=1;fs=8000;wc=2*fc/fs; wp=2*fp/fs; n,ws=buttord(wc,wp,ap,as);b,a=butter(n,ws,high);figure(1);freqz(b,a,512,fs);x=filter(b,a,z1);x=fft(x,8192);figure(2);subplot

25、(2,2,1);plot(abs(y1);axis(0,1000,0,1.0);title(滤波前信号频谱);subplot(2,2,2);plot(abs(x);axis(0,1000,0,1.0);title(滤波后信号频谱);subplot(2,2,3);plot(z1);title(滤波前信号波形);subplot(2,2,4);plot(x);title(滤波前信号波形);sound(x,fs,bits);图形分析如下：选用cheby1程序设计如下：clear;close allz1,fs,bits=wavread(d:laoa.wav)y1=z1(1:8192);y1=fft(y1

26、);fc=2800 ;fp=3000 ;as=100;ap=1; fs=8000;wc=2*fc/fs;wb=2*fp/fs;n,wp=cheb1ord(wc,wb,ap,as);b,a=cheby1(n,ap,wp,high);figure(1);freqz(b,a);x=filter(b,a,z1);x=fft(x,8192);figure(2);subplot(2,2,1);plot(abs(y1);axis(0,1000,0,1.0);title(滤波前信号频谱);subplot(2,2,2);plot(abs(x);axis(0,4000,0,0.03);title(滤波后信号频谱)

27、;subplot(2,2,3);plot(z1);title(滤波前信号波形);subplot(2,2,4);plot(x);title(滤波前信号波形);sound(x,fs,bits);图形分析如下：（6），双线性变换法设计带通滤波器选用butter程序设计如下clear;close allz1,fs,bits=wavread(d:laoa.wav)y1=z1(1:8192);y1=fft(y1);fp1=1200 ;fp2=3000; fc1=1000 ;fc2=3200 ;as=100;ap=1; fs=8000; wc=2*fc1/fs,2* fc2/fs;wp=2*fp1/fs,2

28、*fp2/fs; n,ws=buttord(wc,wp,ap,as);b,a=butter(n,ws,stop);figure(1);freqz(b,a,512,fs);x=filter(b,a,z1);x=fft(x,8192);figure(2);subplot(2,2,1);plot(abs(y1);axis(0,1000,0,1.0);title(滤波前信号频谱);subplot(2,2,2);plot(abs(x);axis(0,1000,0,1.0);title(滤波后信号频谱);subplot(2,2,3);plot(z1);title(滤波前信号波形);subplot(2,2,

29、4);plot(x);title(滤波前信号波形);sound(x,fs,bits);图形分析如下：选用cheby1程序设计如下：clear;close allz1,fs,bits=wavread(d:laoa.wav)y1=z1(1:8192);y1=fft(y1);fp1=1200 ;fp2=3000; fc1=1000 ;fc2=3200 ;as=100;ap=1; fs=8000; wc=2*fc1/fs,2* fc2/fs;wb=2*fp1/fs,2*fp2/fs; n,wp=cheb1ord(wc,wb,ap,as);b,a=cheby1(n,ap,wp,stop);figure(1);freqz(b,a);x=filter(b,a,z1);x=fft(x,8192);figure(2);subplot(2,2,1);plot(abs(y1);axis(0,1000,0,1.0);title(滤波前信号频谱);subplot(2,2,2);plot(abs(x);axis(0,4000,0,0.03);title(滤波后信号频谱);subplot(2,2,3);plot(z1);title(滤波前信号波形);subplot(2,2,4);plot(x);title(滤波前信号波形);sound(x,fs,bits);图形分析如下：四、结果分析：双线性变换法分析：双