数字信号处理综合实验.docx

1、数字信号处理综合实验数字信号处理综合实验一、实验目的1掌握MATLAB的程序设计方法；2掌握数字信号处理的基本理论和基本方法；3掌握语音信号的采集与处理方法；4掌握用MATLAB设计FIR和IIR数字滤波器的方法；5掌握用MATLAB对信号进行分析和处理的方法二、实验原理1.采样定理：在进行模拟/数字信号的转换过程中，当采样频率fs.max大于信号中最高频率fmax的2倍时，即：fs.max=2fmax,则采样之后的数字信号完整地保留了原始信号中的信息，一般实际应用中保证采样频率为信号最高频率的510倍；采样定理又称奈奎斯特定理。2.数字滤波器设计的基本原理：数字滤波器是由常数乘法器、加法器和

2、延迟单元组成的一种算法或装置，是应用最广泛的线性时不变系统。它对输入序列进行处理，去除其频谱中不需要的成分，保留有效成分。（1）对于IIR滤波器，其逼近问题就是寻找滤波器的各项系数，使其系统函数逼近一个所要求的特性。先设计一个合适的模拟滤波器，然后变换成满足约定指标的数字滤波器。 用双线形变换法设计IIR数字滤波器的过程：a.将设计性能指标中的关键频率点进行“预畸” b.利用“预畸”得到的频率点设计一个模拟滤波器。c.双线形变换，确定系统函数（2）FIR滤波器的优点是有N个零点和N个极点，但其全部极点都在z平面的原点，FIR系统总是稳定的；在满足一定的对称条件下，可以实现严格的线性相位。缺点是

3、由于没有非零极点，要获得较好的过渡带特性，必须以较高的滤波器阶数作为代价。（3）IIR滤波器设计可以通过模拟滤波器设计和各种变换法实现，而FIR滤波器不能利用模拟滤波器的设计理论，因为FIR滤波器的系统函数是的多项式，并非有理分式，无法找到与其对应的原型模拟滤波器。3.频谱现象分析：（1）混叠：采样序列的频谱是被采样信号频谱的周期延拓，当采样频率不满足奈奎斯特采样定理的时候，就会发生混叠，使得刺痒后的序列信号的频谱不能真实的反映原采样信号的频谱。（2）泄露：根据理论分析，一个时间的信号其频带宽度为无限，一个时间无限的信号其频带宽度则为有限。因此对一个时间有限的信号，应用DFT进行分析，频谱混叠

4、难以避免。对一个时间无限的信号虽然频带有限，但在实际运算中，时间总是取有限值，在将信号截断的过程中，出现了分散的扩展谱线的现象，称之为频谱泄露或功率泄露。（3）栅栏效应：DFT是对单位圆上Z变换的均匀采样，所以它不可能将频谱视为一个连续函数，就在一定意义上看，用DFT来观察频谱就好象通过一个栅栏来观看一个景象一样，只能在离散点上看到真实的频谱，这样就有可能发生一些频谱的峰点和谷点被“尖桩的栅栏”所挡住，不能被我们观察到。三、实验内容1、连续信号的频率分析一个连续信号含两个频率分量：，用Hz的采样器对该信号进行采样，对采样后的序列分别用长度为16和128的FFT观察其频谱，你观察到什么现象？请解

5、释原因。2、数字滤波器的设计 分别用双线性变换法和窗口设计法设计IIR和FIR数字滤波器，包括低通、高通和带通三种类型。fs=8000HZ;三种滤波器的性能指标如下：（1）低通滤波器：通带和阻带边界频率分别为1000Hz和1200Hz，通带波动为1dB，阻带最小衰减为50dB；（2）高通滤波器：通带和阻带边界频率分别为3000Hz和2800Hz，通带波动为1dB，阻带最小衰减为50dB；（3）带通滤波器：通带边界频率分别为1200Hz和3000Hz，阻带边界频率分别为1000Hz和3200Hz，通带波动为1dB，阻带最小衰减为50dB。请给出每种滤波器的阶数、系统函数和对数幅频响应。3、语音信

6、号的采集与处理（1）用Windows自带的录音机或其他录音软件录制一段语音。录制格式如下：采样频率8000Hz，8位，单声道。录音内容为你的姓名和学号，时长控制在5秒左右，说话不要有过多的停顿。最后格式为.wav;（2）对语音信号进行频谱分析，请画出时域波形和对数幅频特性。（3）分别用“2、数字滤波器的设计”中的三种FIR滤波器对语音信号进行滤波，比较滤波前后语音信号的时域波形与幅频特性。并用MATLAB的sound函数回放语音，试听滤波前后语音的变化。请在实验结果与分析中解释实验现象。（4）分别用“2、数字滤波器的设计”中的IIR低通滤波器和线性相位FIR低通滤波器对语音信号进行滤波，比较两

7、种滤波器滤波后的时域波形与幅频特性。并用MATLAB的sound函数回放语音，试听两种滤波器滤波后语音的差异。请在实验结果与分析中解释实验现象。四、实验结果与分析1.实验结果： 实验结果分析：采样后序列：，N1=8; N2=64；所以N=64;两序列有效长度一样，分辨率一样，由于FFT点数不同，导致结果不同。N=1664时，无频谱泄漏。2.实验结果： (1)az=1.0000 -8.6797 35.5999 -90.8275 159.8643 -203.5543 191.1691 -132.3641 66.2101 -22.7930 4.8632 -0.4877bz=5.2944e-008 5

8、.8238e-007 2.9119e-006 8.7357e-006 1.7471e-005 2.446e-005 2.446e-005 1.7471e-005 8.7357e-006 2.9119e-006 5.8238e-007 5.2944e-008可得IIR滤波器的阶数为12阶。通带和阻带的数字边界频率:中心频率为： 选择窗函数和窗口长度。阻带衰减不小于50dB，因此选择海明窗(其阻带衰减为53dB，满足衰减要求)。滤波器的过渡带宽为，因此窗口长度为:线性相位延迟常数为:根据理想边界频率和线性相位延迟常数，求理想单位脉冲响应:窗函数与理想单位脉冲响应相乘，得到线性相位FIR低通滤波器的

9、单位脉冲响应:可得FIR滤波器的阶数为132阶(2)az=1.0000 8.6797 35.5999 90.8275 159.8643 203.5543 191.1691 132.3641 66.2101 22.7930 4.8632 0.4877bz=5.2944e-008 -5.8238e-007 2.9119e-006 -8.7357e-006 1.7471e-005 -2.446e-005 2.446e-005 -1.7471e-005 8.7357e-006 -2.9119e-006 5.8238e-007 -5.2944e-008可得IIR滤波器的阶数为12阶。通带和阻带的数字边界

10、频率:中心频率为： 选择窗函数和窗口长度。阻带衰减不小于50dB，因此选择海明窗(其阻带衰减为53dB，满足衰减要求)。滤波器的过渡带宽为，因此窗口长度为:线性相位延迟常数为:根据理想边界频率和线性相位延迟常数，求理想单位脉冲响应:窗函数与理想单位脉冲响应相乘，得到线性相位FIR高通滤波器的单位脉冲响应:可得FIR滤波器的阶数是132位。(3)az= 1.0000 0.8124 2.1656 1.6670 3.6730 2.2977 3.5842 1.8949 2.5651 0.9933 1.1179 0.2734 0.2754bz=0.034963 0.013788 -0.07315 -0.

11、010428 0.1368 0.009192 -0.15175 0.009192 0.1368 -0.010428 -0.07315 0.013788 0.034963可得IIR滤波器的阶数为13阶。通带和阻带的数字边界频率:选择窗函数和窗口长度。阻带衰减不小于50dB，因此选择海明窗(其阻带衰减为53dB，满足衰减要求)。滤波器的过渡带宽为，因此窗口长度为:线性相位延迟常数为:根据理想边界频率和线性相位延迟常数，求理想单位脉冲响应:窗函数与理想单位脉冲响应相乘，得到线性相位FIR带通滤波器的单位脉冲响应:3(2)(3)人说话语音在300到3400HZ之间，主要还是集中在低频部分，高频由极少数

12、的的高声调组成。用FIR低通滤波器进行滤波以后，人声基本都还保留着，但是高音调基本没有了，高频部分被滤除掉；用FIR高通滤波器进行滤波以后，基本上听不到人的声音，只剩下一点点模糊的声音和刺耳的噪声，用FIR带通滤波器进行滤波以后，人声比经过高通滤波器滤除要清晰，比经过低通滤波器滤除要模糊一些。(3)IIR滤波器和FIR滤波器的不同在于IIR滤波器阶数少，延迟小，但是有稳定性问题，非线性相位；FIR滤波器没有稳定性问题，线性相位，但阶数多，延迟大，对声音的处理不同于图像的处理，并不需要实现的严格的线性相位，FIR低通滤波器和IIR低通滤波器对语音信号进行滤波后的结果基本相同。五、实验小结六、程序

13、1. n=0:1:63;xa=sin(pi/4*n)+sin(2250/8000*pi*n);h1=abs(fft(xa,16);h2=abs(fft(xa,128);figure;n1=0:1:15;subplot(2,1,1);stem(n1,h1);xlabel(k);ylabel(xk);title(N=16,FFT幅度特性);n2=0:1:127;subplot(2,1,2);stem(n2,h2);xlabel(k);ylabel(xk);title(N=128,FFT幅度特性);2.(1) clear;fs=8000;fc=1000;fr=1200;rp=1;rs=50;%双线性

14、变换法wp=2*fs*tan(2*pi*fc/fs/2);ws=2*fs*tan(2*pi*fr/fs/2);N1, wn =cheb1ord(wp, ws, rp, rs, s);b a= cheby1(N1,rp,wn,low,s);bz,az=bilinear(b,a,fs);h,w=freqz(bz,az);f=w/(2*pi)*fs;figure; plot(f,20*log10(abs(h),-b);grid; xlabel(频率/Hz); ylabel(幅度);legend(双线性变换法);title(切比雪夫低通滤波器，对数幅度谱);%窗口设计法 wc=2*pi*fc/fs;w

15、r=2*pi*fr/fs;%N=rs*6.6*pi/(wc-wr); N=132;n=0:N-1;a=(N-1)/2; wn1=(wc+wr)/2; hd=(sin(n-a)*wn1)./(pi*(n-a); win=hamming(N); h1=win.*hd;H,w1=freqz(h1,1); H=20*log10(abs(H);figure; plot(w1/(pi), H); xlabel(omega/pi); ylabel(幅度/dB); grid;title(线性相位FIR低通滤波器，海明窗，N=132);(2) clear;fs=8000;fc=3000;fr=2800;rp=1

16、;rs=50;%双线性变换法wp=2*fs*tan(2*pi*fc/fs/2);ws=2*fs*tan(2*pi*fr/fs/2);N1, wn =cheb1ord(wp, ws, rp, rs, s);b a= cheby1(N1,rp,wn,high,s);bz,az=bilinear(b,a,fs);h,w=freqz(bz,az);f=w/(2*pi)*fs;figure; plot(f,20*log10(abs(h),-b);grid; xlabel(频率/Hz); ylabel(幅度);legend(双线性变换法);title(切比雪夫高通滤波器，对数幅度谱);%窗口设计法wc=2

17、*fc/fs;wr=2*fr/fs;N=132;n=0:N-1;M=N-1;a=(N-1)/2; wn1=(wc+wr)/2; h1=fir1(M-1,wc,high,hamming(M); H,w1=freqz(h1,1); H=20*log10(abs(H);figure; plot(w1./pi, H); xlabel(omega/pi); ylabel(幅度/dB); grid;title(线性相位FIR高通滤波器，海明窗，N=132);(3) clear;fs=8000;fc=1200 3000;fr=1000 3200;rp=1;rs=50;%双线性变换法wp=2*fs*tan(2

18、*pi*fc/fs/2);ws=2*fs*tan(2*pi*fr/fs/2);N1, wn = ellipord(wp, ws, rp, rs, s);b a= ellip(N1,rp,rs,wn,bandpass,s);bz,az=bilinear(b,a,fs);h,w=freqz(bz,az);f=w/(2*pi)*fs;figure; plot(f,20*log10(abs(h),-b);grid; xlabel(频率/Hz); ylabel(幅度);legend(双线性变换法);title(椭圆带通滤波器，对数幅度谱);%窗口设计法 wc=2*fc/fs; wr=2*fr/fs;N=

19、132;n=0:N-1;M=N-1;a=(N-1)/2; wn1=(wc+wr)/2; h1=fir1(M,wc,bandpass,hamming(N); H,w1=freqz(h1,1); H=20*log10(abs(H);figure; plot(w1/(pi), H); xlabel(omega/pi); ylabel(幅度/dB); grid;title(线性相位FIR带通滤波器，海明窗，N=132);3(2) clear;x=wavread(C:UserslenovoDesktop录音.wav);yk=abs(fft(x);yk1=20*log10(yk);subplot(2,1,

20、1);plot(x);title(原始声音信号时域波形);subplot(2,1,2);plot(yk1);title(对数幅频特性);x=wavread(C:UserslenovoDesktop录音.wav);N=132; n=0:N-1;a=(N-1)/2;wn=0.275*pi; hd=sin(0.275*pi*(n-a)./(pi*(n-a); win=0.54-0.46*cos(n*2*pi/(N-1); %win=hamming(N);h=win.*hd;subplot(2,1,1);y=conv(h,x);plot(y);xlabel(t);ylabel(y(t);title(

21、低通滤波后时域特性);yk=abs(fft(y);subplot(2,1,2);yk1=20*log10(yk);plot(yk1);xlabel(k);ylabel(y(k);axis(0,35000,-200,100);title( 低通滤波后幅频特性);grid;sound(y);x=wavread(C:UserslenovoDesktop录音.wav);N=132;n=0:N-1;a=(N-1)/2; wn=0.725*pi; fs=8000;hd=(sin(pi*(n-a)-sin(wn*(n-a)./(pi*(n-a); win=0.54-0.46*cos(n*2*pi/(N-1)

22、; %win=hamming(N);h=win.*hd;y=conv(h,x);subplot(2,1,1);plot(y);xlabel(t);ylabel(y(t);title( 高通滤波后时域特性);yk=abs(fft(y);yk1=20*log10(yk);subplot(2,1,2);stem(yk1);xlabel(k);ylabel(y(k);title( 高通滤波后幅频特性);grid;sound(y);x=wavread(C:UserslenovoDesktop录音.wav);N=132; n=0:N-1;a=(N-1)/2; w1=0.275*pi;w2=0.775*pi

23、; fs=8000;hd=(sin(w2*(n-a)-sin(w1*(n-a)./(pi*(n-a); win=0.54-0.46*cos(n*2*pi/(N-1); %win=hamming(N);h=win.*hd;y=conv(h,x);subplot(2,1,1);plot(y);xlabel(t);ylabel(y(t);title( 带通滤波后时域特性);yk=abs(fft(y);yk1=20*log10(yk);subplot(2,1,2);stem(yk1);xlabel(k);ylabel(y(k);title( 带通滤波后幅频特性);grid;sound(y);(3) c

24、lear;x=wavread(C:UserslenovoDesktop录音.wav);yk=abs(fft(x);yk1=20*log10(yk); fs=8000;fc=1000;fr=1200;rp=1;rs=50;%双线性变换法wp=2*fs*tan(2*pi*fc/fs/2);ws=2*fs*tan(2*pi*fr/fs/2);N1,wn =cheb1ord(wp, ws, rp, rs, s);b a= cheby1(N1,rp,wn,low,s);bz,az=bilinear(b,a,fs);m=fftfilt(bz,x);M=20*log10(abs(m);subplot(2,1

25、,1);plot(m);title(IIR低通滤波器之后的时域波形);subplot(2,1,2);plot(M);title(IIR低通滤波器之后对数幅频特性);sound(m);clear;x=wavread(C:UserslenovoDesktop录音.wav);yk=abs(fft(x);yk1=20*log10(yk); N=132;M=N-1;n=0:M;a=(N-1)/2;wd=0.275*pi;win=hamming(N);h=fir1(M,wd,low,win);m=fftfilt(h,x);M=20*log10(abs(m);subplot(2,1,1);plot(m);title(FIR低通滤波器之后的时域波形);subplot(2,1,2);plot(M);title(FIR低通滤波器之后对数幅频特性);sound(m);