河海大学数字信号处理实验 综合实验.docx

河海大学数字信号处理实验 综合实验.docx

《河海大学数字信号处理实验 综合实验.docx》由会员分享，可在线阅读，更多相关《河海大学数字信号处理实验 综合实验.docx（29页珍藏版）》请在冰豆网上搜索。

河海大学数字信号处理实验综合实验

数字信号处理综合实验

班级：

姓名：

学号：

一、实验目的

1．掌握MATLAB的程序设计方法；

2．掌握数字信号处理的基本理论和基本方法；

3．掌握语音信号的采集与处理方法；

4．掌握用MATLAB设计FIR和IIR数字滤波器的方法；

5．掌握用MATLAB对信号进行分析和处理的方法；

二、实验原理

y=exp（x）:

以e为底的指数。

conj（x）:

取x的共轭，即改变x的虚部符号。

real（x）:

取复数x的实部。

rand（1,N）:

生成在0和1之间均匀分布的随机序列，长度为N。

randn（1,N）:

生成正态分布（高斯分布）的随机序列，长度为N。

sound（f,fs）:

输入参量是音频数据向量、采样频率和转换位数。

X=fft（x,N）:

计算序列x的N点FFT。

如果x的长度小于N，则在x后面补零；如果x的长度大于N，则对x进行截取；如果不指定参数N，则以x的实际长度作为FFT的点数。

x=ifft（X,N）:

计算序列X的N点IFFT。

Y=fftshift（X）:

将序列X分成左右两部分并交换位置。

[N,Wn]=buttord（Wp,Ws,Rp,Rs,'s'）:

求巴特沃思滤波器的阶数N和3dB边界角频率Wn。

Rp和Rs分别为通带波动δ和阻带衰减At，单位均为dB；'s'表示模拟滤波器设计，如无此参数，则表示数字滤波器设计；Wp和Ws分别表示通带边界角频率Ωc和阻带边界角频率Ωr，其值为标量（低通和高通）或双元素向量（带通和带阻）。

[b,a]=butter（N,Wn,'ftype','s'）:

根据N和Wn求巴特沃思滤波器的系统函数H（s），b和a分别为H（s）分子和分母多项式的系数。

‘ftype’指定滤波器类型，值为‘low’、‘high’或‘stop’，‘low’为缺省值，表示设计低通或带通滤波器（取决于Wn是标量还是双元素向量），‘high’或‘stop’分别表示设计高通和带阻滤波器。

低通和高通滤波器的阶数为N，带通和带阻滤波器的阶数为2N。

[z,p,k]=butter（N,Wn,'ftype','s'）:

根据N和Wn求巴特沃思滤波器的零、极点和增益。

z、p和k分别为H（s）的零点向量、极点向量和增益。

系统函数H（s）与b和a或z、p和k的关系为：

[N,Wp]=cheb1ord（Wp,Ws,Rp,Rs,'s'）:

求切比雪夫I型滤波器的阶数N。

输入参数的含义与巴特沃思滤波器相同，输出参数Wp等于输入参数Wp。

[b,a]=cheby1（N,Rp,Wp,'ftype','s'）:

求切比雪夫I型滤波器的系统函数H（s），b和a分别为H（s）分子和分母多项式的系数。

[z,p,k]=cheby1（N,Rp,Wp,'ftype','s'）:

求切比雪夫I型滤波器的零、极点和增益。

[N,Wp]=ellipord（Wp,Ws,Rp,Rs,'s'）:

求椭圆滤波器的阶数N。

[b,a]=ellip（N,Rp,Rs,Wp,'ftype','s'）:

求椭圆滤波器的系统函数H（s），b和a分别为H（s）分子和分母多项式的系数。

[z,p,k]=ellip（N,Rp,Rs,Wp,'ftype','s'）:

求椭圆滤波器的零、极点和增益。

[bz,az]=impinvar（b,a,fs）:

通过脉冲响应不变法求数字滤波器的系统函数H（z）。

b和a分别为模拟滤波器系统函数H（s）的分子和分母多项式系数；fs为采样频率；bz和az分别为H（z）的分子和分母多项式系数。

[bz,az]=bilinear（b,a,fs）:

通过双线性变换法求数字滤波器的系统函数H（z）。

输入、输出参数的含义与impinvar函数相同。

系统函数H（z）与bz和az的关系为：

window=ones（1,N）:

产生N点矩形窗，行向量。

window=hann（N）:

产生N点汉宁窗，列向量。

window=hanning（N）:

产生N点非零汉宁窗，列向量。

等价于去除hann（N+2）的第一个零元素和最后一个零元素，得到的N点非零窗函数。

window=hamming（N）:

产生N点海明窗，列向量。

window=blackman（N）:

产生N点布莱克曼窗，列向量。

window=kaiser（N,beta）:

产生参数为beta的N点凯塞窗，列向量。

[M,Wd,beta,ftype]=kaiserord（f,a,dev,fs）:

凯塞窗参数估计。

f为一组边界频率，最高频率为fs/2。

a为f中各个频带的幅度值，通带取1，阻带取0。

如果f中有2个元素，则形成3个频带，其中第1个和第3个是通带或阻带，第2个是过渡带，a中也有2个元素，指明第1个和第3个频带是通带还是阻带；如果f中有4个元素，则形成5个频带，其中1，3和5是通带或阻带，2和4是过渡带，a中有3个元素，指明1，3和5是通带还是阻带。

dev的维数与a相同，指明每个频带上的波动值。

fs为采样频率。

M为FIR滤波器的阶数，M=N-1。

Wd为归一化边界频率，等于数字边界角频率除以π，或者边界频率除以fs/2。

beta就是凯塞窗的参数β。

ftype为滤波器的类型。

b=fir1（M,Wd,'ftype',window）:

用窗函数法求FIR滤波器的系数b（单位脉冲响应）。

M为滤波器的阶数，M=N-1。

Wd为一组归一化边界频率，通带和阻带间隔分布，无过渡带；只有一个元素，表示低通或高通滤波器；有两个元素表示带通和带阻滤波器；有三个及以上元素，表示多带滤波器。

'ftype'表示滤波器类型，'high'表示高通滤波器，'stop'表示带阻滤波器，'DC-0'表示多带滤波器的第一个频带为阻带，'DC-1'表示多带滤波器的第一个频带为通带。

window为窗口类型，缺省为海明窗。

b=fir2（M,f,m,window）:

用频率采样法求FIR滤波器的系数b。

M为滤波器的阶数，M=N-1。

f为一组归一化频率，第一个元素必须为0，最后一个元素必须为1（对应奈奎斯特频率，即采样频率的一半），中间的元素按升序排列。

m的维数与f相同，指明f中每个频率上的理想幅度。

window为窗口类型，缺省为海明窗。

Fir2可以实现任意幅度特性的滤波器。

三、实验内容

1、连续信号的频率分析

一个连续信号

含两个频率分量：

，用

Hz的采样器对该信号进行采样，对采样后的序列分别用长度为16和128的FFT观察其频谱，你观察到什么现象？

请解释原因。

2、数字滤波器的设计

分别用双线性变换法和窗口设计法设计IIR和FIR数字滤波器，包括低通、高通和带通三种类型。

三种滤波器的性能指标如下：

（1）低通滤波器：

通带和阻带边界频率分别为1000Hz和1200Hz，通带波动为1dB，阻带最小衰减为50dB，采样频率为8000Hz；

（2）高通滤波器：

通带和阻带边界频率分别为3000Hz和2800Hz，通带波动为1dB，阻带最小衰减为50dB，采样频率为8000Hz；

（3）带通滤波器：

通带边界频率分别为1200Hz和3000Hz，阻带边界频率分别为1000Hz和3200Hz，通带波动为1dB，阻带最小衰减为50dB，采样频率为8000Hz。

请给出每种滤波器的阶数、系统函数和对数幅频响应。

3、语音信号的采集与处理

（1）用Windows自带的录音机或其他录音软件录制一段语音。

录制格式如下：

采样频率8000Hz，8位，单声道。

录音内容为你的姓名和学号，时长控制在5秒左右，说话不要有过多的停顿。

（2）对语音信号进行频谱分析，请画出时域波形和对数幅频特性。

（3）分别用“2、数字滤波器的设计”中的三种FIR滤波器对语音信号进行滤波，比较滤波前后语音信号的时域波形与幅频特性。

并用MATLAB的sound函数回放语音，试听滤波前后语音的变化。

请在实验结果与分析中解释实验现象。

（4）分别用“2、数字滤波器的设计”中的IIR低通滤波器和线性相位FIR低通滤波器对语音信号进行滤波，比较两种滤波器滤波后的时域波形与幅频特性。

并用MATLAB的sound函数回放语音，试听两种滤波器滤波后语音的差异。

请在实验结果与分析中解释实验现象。

四、实验结果与分析

1.

程序：

n1=0:

1:

15;

x=sin（0.25*pi*n1）+sin（0.28125*pi*n1）;

subplot（2,2,1）;

plot（n1,x,'-*'）;

xlabel（'n'）;

ylabel（'x1（n）'）;

title（'N=16,时域特性'）;

subplot（2,2,2）;

xk1=abs（fft（x））;

stem（n1,xk1）

xlabel（'k'）;

ylabel（'X1（k）'）;

title（'N=16,幅频特性'）;

n2=0:

1:

127;

x=sin（0.25*pi*n2）+sin（0.28125*pi*n2）;

subplot（2,2,3）;

plot（n2,x,'-*'）;

xlabel（'n'）;

ylabel（'x2（n）'）;

title（'N=128,时域特性'）;

subplot（2,2,4）;

xk2=abs（fft（x））;

stem（n2,xk2）

xlabel（'k'）;

ylabel（'X2（k）'）;

title（'N=128,幅频特性'）;

结果：

分析：

当采样点N2=128时，幅频特性的拟合性更好。

因为所取的N值越大，即采样的点越多，就会越贴近实际，误差也就越小。

2.

（1）低通滤波器：

IIR数字滤波器：

程序：

clear;

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）;

[N,wn]=ellipord（wp,ws,rp,rs,'s'）;

[ba]=ellip（N,rp,rs,wn,'s'）;

[bz,az]=bilinear（b,a,fs）;

[h,w]=freqz（bz,az）;

f=w/（2*pi）*fs;

plot（f,20*log10（abs（h）））

axis（[0,3000,-120,10]）;

grid;xlabel（'频率/Hz'）;ylabel（'幅度'）;

legend（'双线性变换法'）;

title（'椭圆低通滤波器'）;

结果：

系统函数为：

FIR数字滤波器：

通带和阻带的数字边界频率:

求理想低通滤波器的边界频率用通带和阻带边界频率的中心近似:

滤波器的过渡带宽为

，

阻带衰减不小于50dB，因此选择海明窗，因此窗口长度为:

线性相位延迟常数为:

根据理想边界频率

和线性相位延迟常数

，求理想单位脉冲响应

:

窗函数与理想单位脉冲响应相乘，得到线性相位FIR低通滤波器的单位脉冲响应:

程序：

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;

figure;stem（n,h）;

xlabel（'n'）;ylabel（'h（n）'）;grid;

title（'线性相位FIR低通滤波器的单位脉冲响应h（n）'）;

[H,w]=freqz（h,1）;H=20*log10（abs（H））;

figure;plot（w/pi,H）;

axis（[01-10010]）;xlabel（'频率/Hz'）;ylabel（'幅度/dB'）;grid;

title（'线性相位FIR低通滤波器，海明窗，N=132'）;

结果：

（2）高通滤波器：

IIR数字滤波器：

程序：

clear

fs=8000;fr=2800;fc=3000;rs=50;rp=1;

wp=2*fs*tan（2*pi*fc/fs/2）;

ws=2*fs*tan（2*pi*fr/fs/2）;

[N,wp]=ellipord（wp,ws,rp,rs,'s'）;

[b,a]=ellip（N,rp,rs,wp,'high','s'）

[bz,az]=bilinear（b,a,fs）;

[h,w]=freqz（bz,az）;

f=w/（2*pi）*fs;

plot（f,20*log10（abs（h）））;

axis（[0,4000,-100,10]）;grid;xlabel（'频率/Hz'）;ylabel（'幅度'）;

legend（'双线性变换法'）;

title（'椭圆高通滤波器'）;

结果：

系统函数：

N=6

FIR滤波器：

通带和阻带的数字边界频率:

理想高通滤波器的边界频率用通带和阻带边界频率的中心近似:

滤波器的过渡带宽为

，

阻带衰减不小于50dB，因此选择海明窗因此窗口长度为:

线性相位延迟常数为:

根据理想边界频率

和线性相位延迟常数

，求理想单位脉冲响应

:

窗函数与理想单位脉冲响应相乘，得到线性相位FIR低通滤波器的单位脉冲响应:

程序：

N=132;n=0:

N-1;

a=（N-1）/2;

wn=0.725*pi;

hd=（sin（pi*（n-a））-sin（0.725*pi*（n-a）））./（pi*（n-a））;

win=0.54-0.46*cos（n*2*pi/（N-1））;%win=hamming（N）';

h=win.*hd;

figure;stem（n,h）;

xlabel（'n'）;ylabel（'h（n）'）;grid;

title（'线性相位FIR高通滤波器的单位脉冲响应h（n）'）;

[H,w]=freqz（h,1）;H=20*log10（abs（H））;

figure;plot（w/pi,H）;

axis（[01-10010]）;xlabel（'频率/Hz'）;ylabel（'幅度/dB'）;grid;

title（'线性相位FIR高通滤波器，海明窗，N=132'）;

结果：

（3）带通滤波器

IIR：

程序：

clear;

fs=8000;fc=[12003000];fr=[10003200];rp=1;rs=50;

wp=2*fs*tan（2*pi*fc/fs/2）;

ws=2*fs*tan（2*pi*fr/fs/2）;

[N,wn]=cheb1ord（wp,ws,rp,rs,'s'）;

[b2a2]=cheby1（N,rp,wp,'s'）;

[bz2,az2]=bilinear（b2,a2,fs）;

[h2,w]=freqz（bz2,az2）;

f=w/（2*pi）*fs;

figure;plot（f,20*log10（abs（h2））,'-b'）;

grid;xlabel（'频率/Hz'）;ylabel（'幅度'）;

legend（'双线性变换法'）;

title（'切比雪夫带通滤波器'）;

系统函数：

阶数N=9

FIR：

求理想通带和阻带的数字边界频率和过渡带宽:

滤波器的过渡带宽为

阻带衰减不小于50dB，因此选择汉明窗，因此窗口长度为:

线性相位延迟常数为:

根据理想边界频率

和线性相位延迟常数

，求理想单位脉冲响应

:

窗函数与理想单位脉冲响应相乘，得到线性相位FIR低通滤波器的单位脉冲响应:

程序：

N=132;n=0:

N-1;

a=（N-1）/2;

fs=20000;

w1=0.275*pi;

w2=0.775*pi;

hd=（sin（w2*（n-a））-sin（w1*（n-a）））./（pi*（n-a））;

win=0.54-0.46*cos（n*2*pi/（N-1））;%window=hamming（N）;

h=win.*hd;

figure;stem（n,h）;

xlabel（'n'）;ylabel（'h（n）'）;grid;

title（'线性相位FIR带通滤波器的单位脉冲响应h（n）'）;

[H,w]=freqz（h,1）;H=20*log10（abs（H））;

figure;plot（w/（2*pi）*fs,H）;

axis（[00.5*fs-10010]）;xlabel（'频率/Hz'）;ylabel（'幅度/dB'）;grid;

title（'线性相位FIR带通滤波器，海明窗，N=132'）;

结果：

3.

（2）

程序：

[xfsbits]=wavread（'F:

/姓名学号.wav'）

figure

（1）

plot（x）

xlabel（'t'）;ylabel（'x（t）'）

title（'语音信号时域图像'）

X=fft（x）

figure

（2）

stem（20*log10（abs（X）））

xlabel（'频率/Hz'）;ylabel（'幅度/dB'）;

title（'语音信号对数频谱'）

结果：

（3）

FIR低通滤波：

程序：

[xfsbits]=wavread（'F:

/姓名学号/.wav'）

X=fft（x）

N=132;n=0:

N-1;

a=（N-1）/2;

wn=0.275*pi;

hd=sin（wn*（n-a））./（pi*（n-a））;

win=0.54-0.46*cos（n*2*pi/（N-1））;

h=win.*hd;

y=fftfilt（h,x）

Y=fft（y）

figure

（1）

subplot（211）

plot（x）;title（'原始信号时域'）

subplot（212）

plot（y）;title（'低通滤波后时域'）

figure

（2）

subplot（211）

stem（20*log10（abs（X）））

title（'原始信号对数幅频'）

subplot（212）

stem（20*log10（abs（Y）））

title（'低通滤波后对数幅频'）

sound（y）

结果：

FIR高通滤波：

程序：

[xfsbits]=wavread（'F:

//姓名学号.wav'）

X=fft（x）

N=132;n=0:

N-1;

a=（N-1）/2;

wn=0.725*pi;

hd=sin（pi*（n-a））./（pi*（n-a））-sin（wn*（n-a））./（pi*（n-a））;

win=0.54-0.46*cos（n*2*pi/（N-1））;

h=win.*hd;

y=fftfilt（h,x）

Y=fft（y）

figure

（1）

subplot（211）

plot（x）;title（'原始信号时域'）

subplot（212）

plot（y）;title（'高通滤波后时域'）

figure

（2）

subplot（211）

stem（20*log10（abs（X）））

title（'原始信号对数幅频'）

subplot（212）

stem（20*log10（abs（Y）））

title（'高通滤波后对数幅频'）

sound（y）

结果：

FIR带通滤波：

程序:

[xfsbits]=wavread（'F:

/姓名学号/.wav'）

X=fft（x）

N=132;n=0:

N-1;

hd=（sin（0.775*pi*（n-65.5））-sin（0.275*pi*（n-65.5）））./（pi*（n-65.5））;

window=hamming（N）

h=window'.*hd;

y=fftfilt（h,x）

Y=fft（y）

figure

（1）

subplot（211）

plot（x）;title（'原始信号时域'）

subplot（212）

plot（y）;title（'带通滤波后时域'）

figure

（2）

subplot（211）

stem（20*log10（abs（X）））

title（'原始信号对数幅频'）

subplot（212）

stem（20*log10（abs（Y）））

title（'带通滤波后对视幅频'）

sound（y）

结果：

分析：

低通滤波后声音较低沉；高通滤波后声音较尖锐；带通滤波后声音较贴近正常声音。

但滤波后的声音都没有原声音清亮因为都只是部分声音不完整。

（4）

程序：

X=fft（x）

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）;

[N,wn]=ellipord（wp,ws,rp,rs,'s'）;

[ba]=ellip（N,rp,rs,wn,'s'）;

[bz,az]=bilinear（b,a,fs）;

y1=filter（bz,az,x）

Y1=fft（y1）

N=132;n=0:

N-1;

a=（N-1）/2;

wn=0.275*pi;

hd=sin（wn*（n-a））./（pi*（n-a））;

win=0.54-0.46*cos（n*2*pi/（N-1））;

h=win.*hd;

y2=fftfilt（h,x）

Y2=fft（y2）

figure

（1）

subplot（211）

plot（y1）;title（'IIR低通滤波后时域'）

subplot（212）

plot（y2）;title（'FIR低通滤波后时域'）

figure

（2）

subplot（211）

stem（20*log10（abs（Y1）））

title（'IIR低通滤波后对数幅频'）

subplot（212）

stem（20*log10（abs（Y2）））

title（'FIR低通滤波后对数幅频'）

sound（y1）

sound（y2）

分析：

经过比较可以发现FIR滤波器效果优于IIR滤波器，这是因为其线性更好，可以更好地滤波。

五、实验小结

本次综合实验为既是对前几个实验的总结,而且更贴近于实践，通过对实际语音信号的采集与用MATLAB进行分析与处理，增强了我对用matlab处理信号方面了解。