DSP综合实验.docx
《DSP综合实验.docx》由会员分享，可在线阅读，更多相关《DSP综合实验.docx（15页珍藏版）》请在冰豆网上搜索。
DSP综合实验
一、实验目的
1．学会MATLAB的使用，掌握MATLAB的程序设计方法；
2．掌握在Windows环境下语音信号采集的方法；
3．掌握数字信号处理的基本概念、基本理论和基本方法；
4．掌握MATLAB设计FIR和IIR数字滤波器的方法；
5．学会用MATLAB对信号进行分析和处理。
二．实验原理
1．用窗函数法设计FIR滤波器的基本思想是在时域逼近理想滤波器的单位脉冲响应。
首先根据待逼近理想频率滤波器的频率响应为
。
由IDFT求出理想滤波器的单位脉冲响应h[k]。
为了获得线性相位FIR滤波器在窗函数法设计FIR滤波器的过程中，需要将线性相位因子
加入理想滤波器的频率响应
。
2．利用双线性变换设计IIR滤波器，首先要设计出满足指标要求的模拟滤波器的传递函数
，然后由
通过双线性变换可得所要设计的IIR滤波器的系统函数
。
如果给定的指标为数字滤波器的指标，则首先要转换成模拟域指标。
三．实验内容
1．语音信号的采集
要求利用windows下的录音机或其他软件，录制一段自己的话音，时间控制在1秒左右。
然后在MATLAB软件平台下，利用函数wavread对语音信号进行采样，记住采样频率和采样点数。
通过wavread函数的使用，要求理解采样频率、采样位数等概念。
wavread函数调用格式：
y=wavread（file），读取file所规定的wav文件，返回采样值放在向量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）低通滤波器性能指标，fp=1000Hz，fc=1200Hz，As=100dB，Ap=1dB；2）高通滤波器性能指标，fc=4800Hz，fp=5000HzAs=100dB，Ap=1dB；3）带通滤波器性能指标，fp1=1200Hz，fp2=3000Hz，fc1=1000Hz，c2=3200Hz，As=100dB，Ap=1dB。
要求学生首先用窗函数法设计上面要求的三种滤波器，在MATLAB中，可以利用函数fir1设计FIR滤波器；然后在用双线性变换法设计上面要求的三种滤波器，在MATLAB中，可以利用函数butte、cheby1和ellip设计IIR滤波器；最后，利用MATLAB中的函数freqz画出各滤波器的频率响应。
4．用滤波器对信号进行滤波
比较两种滤波器的性能，然后用性能好的各滤波器分别对采集的信号进行滤波，在MATLAB中，FIR滤波器利用函数fftfilt对信号进行滤波，IIR滤波器利用函数filter对信号进行滤波。
5．比较滤波前后语音信号的波形及频谱要求在一个窗口同时画出滤波前后的波形及频谱。
6．回放语音信号
在MATLAB中，函数sound可以对声音进行回放。
其调用格式：
sound（x,fs,bits）；
可以感觉滤波前后的声音有变化。
四、实验结果及实验分析
1.语音信号的采集及频谱分析
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=1200Hz，As=100dB，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]）;gridon
实验分析：
由图象可以看出，窗函数低通滤波器具有很大的阻带波纹，在通带基本没有波纹。
（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））;
axis（[0,3000,0,1.5]）
grid;
xlabel（'频率/HZ'）;
ylabel（'幅度'）;
title（'Batterworth低通IIR数字滤波器'）;
实验分析：
由实验可看出巴特沃斯滤波器的过渡带很宽。
2）高通滤波器性能指标，fc=4800Hz，fp=5000HzAs=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.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]）;gridon
实验分析：
由图像可知，Kaiser窗函数高通滤波器在通带基本没波纹，过渡带较窄，基本上能符合要求。
（2）用Batterworth设计IIR高通滤波器：
fp=5000;
fc=4800;
Fs=15000;
wc=fc*pi;
Rp=1;
Rs=100;
e=fzero（inline（'（10^（0.1*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=1200Hz，fp2=3000Hz，fc1=1000Hz，c2=3200Hz，As=100dB，Ap=1dB。
（1）用Kaiser窗函数设计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,deltaw,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）用ChebyshevI设计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,wn]=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滤波器的更加接近理想矩形，过渡带窄，通带基本无波纹，衰减基本上能符合要求，但是在阻带的波纹值较大。
巴特沃兹滤波器的过渡带很宽，很难满足要求。
切比雪夫I型滤波器过渡带较窄，但是在通带上有时有较大波纹值。
综合考虑，选择基于Kaiser窗函数的FIR低通滤波器、基于Kaiser窗函数的FIR高通滤波器、基于Kaiser窗函数的FIR带通滤波器为最佳。
（1）使用Kaiser窗函数FIR低通滤波器对声音信号进行滤波
（2） 使用Kaiser窗函数FIR高通滤波器对声音信号进行滤波
（3）使用Kaiser窗函数FIR带通滤波器对声音信号进行滤波
由滤波前后的对比我们可以发现，基于Kaiser窗函数的滤波器基本上能让通带内的信号通过，滤掉阻带内的信号。
5.回放语音信号
回放语音信号之后可以发现，声音发生了变化。
证明了滤波的效果。
经过低通滤波器之后的声音变得低沉了很多，经过高通滤波器之后声音消失不见。
经过带通滤波器之后的声音也发生变化，跟原声也比较接近。
这是由声音本身的频谱特性和滤波器的特性决定的。
五、实验思考
1．双线性变换法中Ω和ω之间的关系是非线性的，在实验中你注意到这种非线性关系了吗?
从那几种数字滤波器的幅频特性曲线中可以观察到这种非线性关系？
答：
在双线性变换法中，模拟频率与数字频率不再是线性关系，所以一个线性相位模拟器经过双线性变换后得到的数字滤波器不再保持原有的线性相位了。
如以上实验过程中，采用双线性变化法设计的butter和cheby1数字滤波器，从图中可以看到这种非线性关系。
2．能否利用公式完成脉冲响应不变法的数字滤波器设计？
为什么？
答：
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/fs;
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]）;gridon
x=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,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=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,-200,0]）;gridon
x=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'）;
sound（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.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]）;
gridon
x=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）;
升级会员 