DSP综合实验.docx

1、DSP综合实验一、实验目的1学会 MATLAB 的使用，掌握 MATLAB 的程序设计方法；2掌握在 Windows 环境下语音信号采集的方法；3掌握数字信号处理的基本概念、基本理论和基本方法；4掌握 MATLAB 设计 FIR 和 IIR 数字滤波器的方法；5学会用 MATLAB 对信号进行分析和处理。二实验原理 1用窗函数法设计FIR滤波器的基本思想是在时域逼近理想滤波器的单位脉冲响应。首先根据待逼近理想频率滤波器的频率响应为。由IDFT求出理想滤波器的单位脉冲响应hk。为了获得线性相位FIR滤波器在窗函数法设计FIR滤波器的过程中，需要将线性相位因子加入理想滤波器的频率响应。 2利用双线

2、性变换设计IIR滤波器，首先要设计出满足指标要求的模拟滤波器的传递函数，然后由通过双线性变换可得所要设计的IIR滤波器的系统函数。如果给定的指标为数字滤波器的指标，则首先要转换成模拟域指标。三实验内容 1语音信号的采集要求利用 windows 下的录音机或其他软件，录制一段自己的话音，时间控制在 1 秒左右。然后在 MATLAB 软件平台下，利用函数 wavread 对语音信号进行采样，记住采样频率和采样点数。通过 wavread 函数的使用，要求理解采样频率、采样位数等概念。wavread 函数调用格式：y=wavread(file)，读取 file 所规定的 wav 文件，返回采样值放在向

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

4、器性能指标， fp=1000Hz，fc=1200 Hz， As=100dB， Ap=1dB； 2）高通滤波器性能指标， fc=4800 Hz， fp=5000 HzAs=100dB， Ap=1dB； 3）带通滤波器性能指标， fp1=1200 Hz， fp2=3000 Hz， fc1=1000 Hz，c2=3200 Hz， As=100dB， Ap=1dB。要求学生首先用窗函数法设计上面要求的三种滤波器，在 MATLAB 中，可以利用函数 fir1 设计 FIR 滤波器；然后在用双线性变换法设计上面要求的三种滤波器，在 MATLAB 中，可以利用函数 butte、 cheby1 和 ellip

5、 设计 IIR 滤波器；最后，利用 MATLAB 中的函数 freqz画出各滤波器的频率响应。 4用滤波器对信号进行滤波比较两种滤波器的性能，然后用性能好的各滤波器分别对采集的信号进行滤波，在MATLAB 中， FIR 滤波器利用函数 fftfilt 对信号进行滤波， IIR 滤波器利用函数 filter 对信号进行滤波。 5比较滤波前后语音信号的波形及频谱要求在一个窗口同时画出滤波前后的波形及频谱。 6回放语音信号在 MATLAB 中，函数 sound 可以对声音进行回放。其调用格式：sound(x,fs,bits)；可以感觉滤波前后的声音有变化。四、实验结果及实验分析1.语音信号的采集及频

6、谱分析Matlab程序代码如下：y,fs,bits=wavread(dsp.wav)y=y(:,1);sound(y,fs,bits);Y=fft(y,fs);subplot(211);plot(y);title(音频信号的时域波形 );subplot(212);plot(abs(Y);title(音频信号的频域波形);axis(0,5000,0,1.2*max(abs(Y);1.得到的结果如下：说明录音文件的采样频率为44100HZ，采样位数为16位。2.信号的时域和频域特性：3.设计数字滤波器和画出频率响应1）低通滤波器性能指标， fp=1000Hz，fc=1200 Hz， As=100d

7、B， Ap=1dB。（1）用Kaiser窗函数设计FIR低通滤波器：fp=1000;fc=1200;As=100;Ap=1;wc=2*pi*fc/fs;wp=2*pi*fp/fs;ws=wc-wp;deltaw=(wc+wp)/(2*pi);n0=ceil(As-8)/2.285/ws);beta=0.1102*(As-8.7);wdkai=(kaiser(n0+1,beta);b0=fir1(n0,deltaw,wdkai);figure(1);h,w=freqz(b0,1);plot(w*fs*0.5/pi,20*log10(abs(h);axis(0,4000,-200,0);grid

8、on 实验分析：由图象可以看出，窗函数低通滤波器具有很大的阻带波纹，在通带基本没有波纹。（2）用Batterworth设计IIR低通滤波器： fp=1000;fc=1200;Fs=22050;wc=fc*pi;Rp=1;Rs=100;e=fzero(inline(10(0.1*1)-1)x-(1000/1200),1);N=1/(2*e);N=ceil(N);Z,P,K=buttap(N);Bap,Aap=zp2tf(Z,P,K);b,a=lp2lp(Bap,Aap,wc);bz,az=bilinear(b,a,Fs);H,W=freqz(bz,az);plot(W*Fs/pi,abs(H);

9、axis(0,3000,0,1.5)grid;xlabel(频率/HZ);ylabel(幅度);title(Batterworth低通IIR数字滤波器 ); 实验分析：由实验可看出巴特沃斯滤波器的过渡带很宽。2）高通滤波器性能指标， fc=4800 Hz， fp=5000 HzAs=100dB， Ap=1dB。 （1）用Kaiser窗函数设计FIR高通滤波器： fc=4800;fp=5000;As=100;Ap=1;wc=2*pi*fc/fs;wp=2*pi*fp/fs;deltaw=(wc+wp)/(2*pi); ws=wp-wc;n0=ceil(As-8)/2.285/ws);beta=0

10、.1102*(As-8.7);wdkai=(kaiser(n0,beta);b0=fir1(n0-1,deltaw,high,wdkai);figure(1);h,w=freqz(b0,1);plot(w*fs*0.5/pi,20*log10(abs(h);axis(0,8000,-200,0);grid on 实验分析：由图像可知，Kaiser窗函数高通滤波器在通带基本没波纹，过渡带较窄，基本上能符合要求。（2）用Batterworth设计IIR高通滤波器： fp=5000;fc=4800;Fs=15000;wc=fc*pi;Rp=1;Rs=100;e=fzero(inline(10(0.1

11、*1)-1)x-(4800/5000),1);N=1/(2*e);N=ceil(N)b,a=butter(N,wc,high,s);bz,az=bilinear(b,a,Fs);H,W=freqz(bz,az);plot(W*Fs/pi,abs(H);grid;xlabel(频率/HZ);ylabel(幅度);title(Batterworth高通IIR数字滤波器);实验分析：由图中可知，巴特沃兹滤波器具有很宽的过渡带。3）带通滤波器性能指标， fp1=1200 Hz， fp2=3000 Hz， fc1=1000 Hz，c2=3200 Hz， As=100dB， Ap=1dB。(1) 用Kai

12、ser窗函数设计FIR高通滤波器：fp1=1200;fp2=3000;fc1=1000;fc2=3200;As=100;Ap=1;wc1=2*pi*fc1/fs;wp1=2*pi*fp1/fs;wc2=2*pi*fc2/fs;wp2=2*pi*fp2/fs;deltaw1=(wc1+wp1)/(2*pi);deltaw2=(wc2+wp2)/(2*pi);ws=wp1-wc1;deltaw=deltaw1,deltaw2;n0=ceil(As-7.95)/2.285/ws);beta=0.1102*(As-8.7);wdkai=(kaiser(n0+1,beta);b0=fir1(n0,del

13、taw,wdkai);figure(1);h,w=freqz(b0,1);plot(w*fs*0.5/pi,20*log10(abs(h);axis(0,4000,-200,0);实验分析：由图像可知，Kaiser窗函数带通滤波器具有较窄的过渡带，通带上基本没波纹，阻带波纹较大。（2）用Chebyshev I设计IIR带通滤波器： fp1=1200;fp2=3000;fc1=1000;fc2=3200;Fs=22050;wc1=fc1*2/Fs;wc2=fc2*2/Fs;wp1=fp1*2/Fs;wp2=fp2*2/Fs;Rp=1;Rs=10;wp=wp1,wp2;wc=wc1,wc2;N,w

14、n=cheb1ord(wp,wc,Rp,Rs,s);bz,az=cheby1(N,Rp,wc,bandpass);H,W=freqz(bz,az);plot(W*Fs/(2*pi),abs(H);grid;xlabel(频率/HZ);ylabel(幅度);title(chebyshev1带通IIR数字滤波器); 实验分析：由图像可知，切比雪夫I型带通滤波器过渡带较窄，通带上有波纹但是较小。 4.用滤波器对信号进行滤波 对以上的图像分析我们发现，FIR滤波器的更加接近理想矩形，过渡带窄，通带基本无波纹，衰减基本上能符合要求，但是在阻带的波纹值较大。巴特沃兹滤波器的过渡带很宽，很难满足要求。切比雪

15、夫I型滤波器过渡带较窄，但是在通带上有时有较大波纹值。综合考虑，选择基于Kaiser窗函数的FIR低通滤波器、基于Kaiser窗函数的FIR高通滤波器、基于Kaiser窗函数的FIR带通滤波器为最佳。（1）使用Kaiser窗函数FIR低通滤波器对声音信号进行滤波（2）使用Kaiser窗函数FIR高通滤波器对声音信号进行滤波（3）使用Kaiser窗函数FIR带通滤波器对声音信号进行滤波 由滤波前后的对比我们可以发现，基于Kaiser窗函数的滤波器基本上能让通带内的信号通过，滤掉阻带内的信号。 5.回放语音信号回放语音信号之后可以发现，声音发生了变化。证明了滤波的效果。经过低通滤波器之后的声音变得

16、低沉了很多，经过高通滤波器之后声音消失不见。经过带通滤波器之后的声音也发生变化，跟原声也比较接近。这是由声音本身的频谱特性和滤波器的特性决定的。五、实验思考1双线性变换法中 和 之间的关系是非线性的，在实验中你注意到这种非线性关系了吗?从那几种数字滤波器的幅频特性曲线中可以观察到这种非线性关系？答：在双线性变换法中，模拟频率与数字频率不再是线性关系，所以一个线性相位模拟器经过双线性变换后得到的数字滤波器不再保持原有的线性相位了。如以上实验过程中，采用双线性变化法设计的butter和cheby1数字滤波器，从图中可以看到这种非线性关系。 2能否利用公式完成脉冲响应不变法的数字滤波器设计？为什么？

17、答：IIR数字滤波器的设计实际上是求解滤波器的系数和，它是数学上的一种逼近问题，即在规定意义上（通常采用最小均方误差准则）去逼近系统的特性。如果在S平面上去逼近，就得到模拟滤波器；如果在z平面上去逼近，就得到数字滤波器。但是它的缺点是，存在频率混迭效应，故只适用于阻带的模拟滤波器。附：Kaiser窗滤波器对语音信号滤波代码：低通滤波：y,fs,bits=wavread(dsp.wav);sound(y,fs,bits);n=length(y);Y=fft(y,n);fp=1000;fc=1200;As=100;Ap=1;%给定指标wc=2*pi*fc/fs;%归一化频率wp=2*pi*fp/f

18、s;ws=wc-wp;deltaw=(wc+wp)/(2*pi);%归一化过渡带宽计算n0=ceil(As-8)/2.285/ws);%计算所需滤波器长度beta=0.1102*(As-8.7);%计算kaiser窗的值wdkai=(kaiser(n0+1,beta);%kaiser窗函数b0=fir1(n0,deltaw,wdkai);figure(1);h,w=freqz(b0,1);plot(w*fs*0.5/pi,20*log10(abs(h);axis(0,4000,-200,0);grid onx=fftfilt(b0,y);X=fft(x,fs);figure(2);subplo

19、t(2,2,1);plot(abs(Y);axis(0,5000,0,200);title(滤波前信号频谱);subplot(2,2,2);plot(abs(X);axis(0,5000,0,200);title(滤波后信号频谱);subplot(2,2,3);plot(y);title(滤波前信号波形);subplot(2,2,4);plot(x);title(滤波前信号波形);sound(x,fs,bits);高通滤波：y,fs,bits=wavread(dsp.wav);sound(y,fs,bits);n=length(y);Y=fft(y,n);fc=4800;fp=5000;As=

20、100;Ap=1;%给定指标wc=2*pi*fc/fs;%归一化频率wp=2*pi*fp/fs;deltaw=(wc+wp)/(2*pi);%截止频率计算ws=wp-wc;%归一化带宽计算n0=ceil(As-8)/2.285/ws);%计算所需滤波器长度beta=0.1102*(As-8.7);%计算kaiser窗的值wdkai=(kaiser(n0,beta);%kaiser窗函数b0=fir1(n0-1,deltaw,high,wdkai);figure(1);h,w=freqz(b0,1);plot(w*fs*0.5/pi,20*log10(abs(h);axis(0,8000,-20

21、0,0);grid onx=fftfilt(b0,y);X=fft(x,fs);figure(2);subplot(2,2,1);plot(abs(Y);axis(0,5000,0,200);title(滤波前信号频谱);subplot(2,2,2);plot(abs(X);axis(0,8000,0,50);title(滤波后信号频谱);subplot(2,2,3);plot(y);title(滤波前信号波形);subplot(2,2,4);plot(x);title(滤波前信号波形);sound(x,fs,bits);带通滤波：y,fs,bits=wavread(dsp.wav);soun

22、d(y,fs,bits);n=length(y);Y=fft(y,n);fp1=1200;fp2=3000;fc1=1000;fc2=3200;As=100;Ap=1;%给定指标wc1=2*pi*fc1/fs;%归一化频率wp1=2*pi*fp1/fs;wc2=2*pi*fc2/fs;wp2=2*pi*fp2/fs;deltaw1=(wc1+wp1)/(2*pi);%截止频率计算deltaw2=(wc2+wp2)/(2*pi);ws=wp1-wc1;%归一化带宽deltaw=deltaw1,deltaw2;n0=ceil(As-7.95)/2.285/ws);%计算所需滤波器长度beta=0.

23、1102*(As-8.7);%计算kaiser窗的值wdkai=(kaiser(n0+1,beta);%kaiser窗函数b0=fir1(n0,deltaw,wdkai);figure(1);h,w=freqz(b0,1);plot(w*fs*0.5/pi,20*log10(abs(h);axis(0,4000,-200,0);grid onx=fftfilt(b0,y);X=fft(x,fs);figure(2);subplot(2,2,1);plot(abs(Y);axis(0,5000,0,150);title(滤波前信号频谱);subplot(2,2,2);plot(abs(X);axis(0,5000,0,150);title(滤波后信号频谱);subplot(2,2,3);plot(y);title(滤波前信号波形);subplot(2,2,4);plot(x);title(滤波前信号波形);sound(x,fs,bits);