数字信号处理实验指导书.docx
《数字信号处理实验指导书.docx》由会员分享，可在线阅读，更多相关《数字信号处理实验指导书.docx（76页珍藏版）》请在冰豆网上搜索。
数字信号处理实验指导书
数字信号处理
实验指导书
王玉富编著
实验一抽样定理
一、实验题目
抽样定理
二、实验目的
⒈掌握MATLAB相关函数的运用；
⒉掌握用MATLAB软件来演示抽样定理；
⒊加深对抽样定理的理解；
三、具体实验内容及相关结论
⒈模拟信号与抽样
一般的MATLAB不能用来处理和分析模拟信号（具有符号工具箱的MATLAB除外），但是可以用近似理想抽样（
远小于T）的方法处理。
⑴已知
，求出并绘制傅立叶变换图。
利用
，可以用一个在
范围之间的有限长信号近似。
因为，当
时，有
，则可选
这样就可以使用MATLAB解题了。
程序Example1_1如下：
%模拟信号
Dt=0.00005;
t=-0.005:
Dt:
0.005;
xa=exp（-1000*abs（t））;
%连续时间傅立叶变换
Wmax=2*pi*2000;
K=500;
k=0:
1:
K;
W=k*Wmax/K;
Xa=xa*exp（-j*t'*W）*Dt;
Xa=real（Xa）;
W=[-fliplr（W）,W（2:
501）];%频率从-Wmax到Wmax
Xa=[fliplr（Xa）,Xa（2:
501）];%Xa介于-Wmax到Wmax之间
subplot（211）,plot（t*1000,xa）;
xlabel（'t毫秒'）;
ylabel（'xa（t）'）
title（'模拟信号'）
subplot（212）,plot（W/（2*pi*1000）,Xa*1000）;
xlabel（'频率（kHz）'）;
ylabel（'Xa（jW）*1000'）;
title（'连续时间傅立叶变换'）
实验结果如图1-1所示。
图1-1example1-1运行结果
⑵为了研究抽样结果，
用两种不同的抽样频率进行抽样。
①抽样频率为
，对
抽样得到
，求出并画出
；
②抽样频率为
，对
抽样得到
，求出并画出
；
①因为
的带宽为
，两倍频率小于
，所以不存在混叠现象。
程序Example1_2_1如下：
%模拟信号
Dt=0.00005;
t=-0.005:
Dt:
0.005;
xa=exp（-1000*abs（t））;
%连续时间傅立叶变换
Ts=0.0002;
n=-25:
1:
25;
x=exp（-1000*abs（n*Ts））;
%离散傅立叶变换
K=500;
k=0:
1:
K;
w=pi*k/K;
X=x*exp（-j*n'*w）;
X=real（X）;
w=[-fliplr（w）,w（2:
K+1）];%频率从-Wmax到Wmax
X=[fliplr（X）,X（2:
K+1）];%Xa介于-Wmax到Wmax之间
subplot（211）,plot（t*1000,xa）;
xlabel（'n毫秒'）;
ylabel（'x1（n）'）;
title（'离散信号'）;
holdon
stem（n*Ts*1000,x）;
gtext（'Ts=0.2毫秒'）;
holdoff
subplot（212）,plot（w/pi,X）;
xlabel（'频率（pi）'）;
ylabel（'X1（w）'）;
title（'离散时间傅立叶变换'）
实验结果如图1-2所示。
离散信号叠加在模拟信号图上，可以清楚地看出抽样过程，很显然不存在混叠现象。
图1-2example1_2-1运行结果
②因为
的带宽为
，两倍频率大于
，所以存在混叠现象。
因为抽样频率不同，实现时只需将Example1_2_1中的条件改为Ts=0.001；n=-5:
1:
5；gtext（'Ts=1毫秒'）就可以了。
运行结果如图1-3所示，其形状不同于图1-2就是因为存在混叠现象。
图1-3混叠现象
⒉信号恢复
理想低通滤波器的时域输出为：
给出的是一无限阶内插系统，在实际中是无法实现的，因此只能用有限阶内插逼近。
MATLAB提供了在相邻点内插的几种方法。
⑴
函数：
根据这个公式就可以用矩阵与向量的相乘运算实现内插了。
利用
函数恢复出
，程序Example1_3如下：
%由x1恢复xa
Ts=0.0002;
n=-25:
1:
25;
nTs=n*Ts;
x1=exp（-1000*abs（nTs））;
subplot（221）,stem（nTs,x1）
title（'x1'）;
%恢复模拟信号
Dt=0.00005;
t=-0.005:
Dt:
0.005;
Fs=1/Ts;
xa=x1*sinc（Fs*（ones（length（n）,1）*t-nTs'*ones（1,length（t））））;
subplot（222）,plot（t,xa）
title（'由x1恢复xa'）;
%校验
error1=max（abs（xa-exp（-1000*abs（t））））
%由x2恢复xa
Ts=0.001;
n=-5:
1:
5;
nTs=n*Ts;
x=exp（-1000*abs（nTs））;
subplot（223）,stem（nTs,x）
title（'x2'）;
%模拟信号恢复
Dt=0.00005;
t=-0.005:
Dt:
0.005;
Fs=1/Ts;
xa=x*sinc（Fs*（ones（length（n）,1）*t-nTs'*ones（1,length（t））））;
subplot（224）,plot（t,xa）
title（'由x2恢复xa'）;
%校验
error2=max（abs（xa-exp（-1000*abs（t））））
校验结果：
error1=0.0363；error2=0.1852，运行结果如图1-4所示。
图1-4sinc（x）函数恢复原始信号
由程序结果可见，恢复信号与原始信号存在误差，这是由于采用了有限个样本逼近无限阶内插所造成的，对于
恢复的
，误差只有0.0363，应该说是效果不错的。
而由样本
恢复出
，结果为error2=0.1852，显然信号恢复效果比较差，这是因为抽样频率不满足奈奎斯特抽样定理，内插后存在混叠的缘故。
⑵图解法：
图解法就是给定样本值，而后利用stairs函数画出阶梯形的模拟信号逼近，或者利用plot函数画出内插曲线图。
程序Example1_4如下：
Ts=0.0002;
n=-25:
1:
25;
nTs=n*Ts;
x=exp（-1000*abs（nTs））;
%利用stairs函数恢复模拟信号
subplot（211）,stairs（nTs*1000,x）;
xlabel（'t毫秒'）;
ylabel（'xa（t）'）
title（'利用stairs函数恢复模拟信号'）;
holdon
stem（n*Ts*1000,x）;
holdoff%画出样本值进行比较
%利用plot函数恢复模拟信号
subplot（212）,plot（nTs*1000,x）;
xlabel（'t毫秒'）;
ylabel（'xa（t）'）
title（'利用plot函数恢复模拟信号'）;
holdon
stem（n*Ts*1000,x）;
holdoff%画出样本值进行比较
图1-5图解法恢复原始信号
程序运行结果如图1-5所示，由图可见，利用plot函数绘出的图形比用stairs函数绘出的图形好，只要抽样符合奈奎斯特抽样定理，则用图解法可以得到近似解。
⑶spline函数：
这种方法是采用一组称为3次样条函数的分段多项式实现，3次样条多项式为
，
其中
是多项式系数，是由样本值的最小二乘分析获得（严格说来，这不是一个因果性运算），这在MATLAB中是很方便得出的。
利用spline函数由
恢复出
。
程序Example1_5如下：
s=0.0002;
n=-25:
1:
25;
nTs=n*Ts;
x=exp（-1000*abs（nTs））;
%模拟信号恢复
Dt=0.00005;
t=-0.005:
Dt:
0.005;
xa=spline（nTs,x,t）;
plot（t,xa）
title（'由x1恢复xa'）
%校验
error=max（abs（xa-exp（-1000*abs（t））））
图1-6spline函数恢复原始信号
程序运行结果：
error=0.0317，由图1-6可以看出，使用spline函数恢复模拟信号仍然存在误差，但其效果要好于其它方法。
实验二序列的傅立叶变换
一、实验题目
序列的傅立叶变换
二、实验目的
⒈掌握相关的MATLAB函数；
⒉用MATLAB来实现序列的傅立叶变换；
⒊利用MATLAB进行序列的对称性分析。
三、具体实验内容及相关结论
⒈已知
，用图解法讨论其傅立叶变换
的对称性。
程序Example2_1如下：
图2-1程序Example2_1运行结果
subplot（1,1,1）
n=-5:
5;
x=（-0.9）.^n;
k=-200:
200;
w=（pi/100）*k;
X=x*（exp（-j*pi/100））.^（n'*k）;
magX=abs（X）;angX=angle（X）;
subplot（211）,plot（w/pi,magX）;
axis（[-2,2,0,15]）
xlabel（'频率（以pi为单位）'）;
ylabel（'X的绝对值'）;
title（'幅度'）;
subplot（2,1,2）;plot（w/pi,angX）;
axis（[-2,2,-4,4]）
xlabel（'频率（以pi为单位）'）;
ylabel（'弧度/pi'）
title（'相角'）
程序运行结果如图2-1所示从中可以看出，
对
对称，因此，分析时只看到0到π之间的模和相位就可以了。
⒉验证实序列的对称特性，令
程序Example2_1如下：
n=-5:
10;x=sin（pi*n/2）;
k=-100:
100;
w=（pi/100）*k;%-pi到pi之间的频率
X=x*（exp（-j*pi/100））.^（n'*k）;%x的DTFT
%信号分解
[xe,xo,m]=evenodd（x,n）;%利用evenodd分解奇偶部
XE=xe*（exp（-j*pi/100））.^（m'*k）;%xe的DTFT
XO=xo*（exp（-j*pi/100））.^（m'*k）;%xo的DTFT
%校验
XR=real（X）;%X的实部
error1=max（abs（XE-XR））%XE与XR是否相等
XI=imag（X）;%X的虚部
error2=max（abs（XO-j*XI））%XO与XI是否相等
%用图形验证
subplot（221）,plot（w/pi,XR）;
gridon
axis（[-1,1,-2,2]）
xlabel（'频率（以pi为单位）'）;
ylabel（'Re（X）'）;
title（'X的实部'）;
subplot（222）,plot（w/pi,XI）;
gridon
axis（[-1,1,-10,10]）
xlabel（'频率（以pi为单位）'）;
ylabel（'Im（X）'）;
title（'X的虚部'）
subplot（223）,plot（w/pi,real（XE））;
gridon
axis（[-1,1,-2,2]）
xlabel（'频率（以pi为单位）'）;
ylabel（'XE'）;
title（'x偶部的DTFT'）
subplot（224）,plot（w/pi,imag（XO））;
gridon
axis（[-1,1,-10,10]）
xlabel（'频率（以pi为单位）'）;
ylabel（'XO'）;
title（'x奇部的DTFT'）
程序运行结果为error1=6.4666e-014，error2=6.4661e-014，如图2-2所示，可见
的实部等于
的DTFT，而
的虚部等于
的DTFT
图2-2程序Example2_2运行结果
⒊evenodd函数
function[xe,xo,m]=evenodd（x,n）
%将实信号分解为奇偶两部分
%检查序列是否为实序列
ifany（imag（x）~=0）
error（'x不是实序列'）
end
m=-fliplr（n）;%fliplr:
Flipmatrixinleft/rightdirection.
m1=min（[m,n]）;m2=max（[m,n]）;m=m1:
m2;
nm=n
（1）-m
（1）;n1=1:
length（n）;
x1=zeros（1,length（m））;
x1（n1+nm）=x;x=x1;
xe=0.5*（x+fliplr（x））;
xo=0.5*（x-fliplr（x））;
实验三IIR数字滤波器设计
一、实验题目
IIR数字滤波器设计
二、实验目的
⒈利用MATLAB软件设计模拟低通滤波器；
⒉利用MATLAB软件实现冲激响应不变法设计IIR数字滤波器；
⒊利用MATLAB软件实现双线性变换法设计IIR数字滤波器；
⒋利用MATLAB软件实现IIR数字滤波器的频带转换。
三、具体实验内容及相关结论
⒈模拟低通滤波器设计
MATLAB提供了N阶归一化（即
）巴特沃斯模拟滤波器设计函数[z,p,k]=buttap（N），其返回值为数组z和p（零点和极点）以及增益k。
MATLAB扩展函数提供的[b,a]=u_buttap函数提供了一个直接型滤波器结构，而sdir2cas给除了适合于模拟滤波器的直接型到级联型的转换。
巴特沃斯低通模拟滤波器的设计过程可以使用MATLAB扩展函数中的afd_butt实现，用freq_m函数显示滤波器的频域图。
程序Example3_1如下：
Wp=0.2*pi;
Ws=0.3*pi;
Rp=7;
Ar=16;%设计指标
[b,a]=afd_butt（Wp,Ws,Rp,Ar）;%模拟滤波器设计
wmax=0.5*pi;
[db,mag,pha,w]=freqs_m（b,a,wmax）;%计算频率响应和相位响应
[ha,x,t]=impulse（b,a）;
[C,B,A]=sdir2cas（b,a）%直接型到级联型转换
%画图
figure
（1）
subplot（331）
plot（mag）
axis（[050001.05]）;
title（'频率响应'）
subplot（332）
plot（pha）
axis（[0500-3.23.2]）;
title（'相位响应'）
subplot（333）
plot（ha）
axis（[0100-0.0250.21]）;
title（'脉冲响应'）
滤波器响应如图3-1所示。
图3-1程序Example3_1运行结果
常用的模拟滤波器设计方法除了巴特沃斯滤波器外，还有切比雪夫Ⅰ型、切比雪夫Ⅱ型、椭圆滤波器等。
MATLAB提供了设计N阶归一化切比雪夫Ⅰ型滤波器的cheb1ap函数，设计N阶归一化切比雪夫Ⅱ型滤波器的cheb2ap函数，设计N阶归一化椭圆滤波器的ellipap函数。
MATLAB扩展函数中提供了设计N阶未归一化切比雪夫Ⅰ型滤波器的u_cheb1ap函数和根据给定指标设计切比雪夫Ⅰ型滤波器的afd_cheb1函数；设计N阶未归一化切比雪夫Ⅱ型滤波器的u_cheb2ap函数和根据给定指标设计切比雪夫滤Ⅱ型波器的afd_cheb2ap函数；设计N阶椭圆模拟滤波器的u_elipap函数和根据给定指标设计椭圆模拟滤波器的afd_elip函数。
⑴afd_butt函数
function[b,a]=afd_butt（Wp,Ws,Rp,Ar）
%按给定的技术指标设计模拟巴特沃思低通滤波器
%b直接型滤波器分子系数；a直接型滤波器分母系数
%Wp通带截止频率，单位弧度/秒；Ws阻带截止频率，单位弧度/秒；
%Rp通带中的振幅波动系数（dB）；Ar阻带衰减，（dB）；
ifWp<=0
error（'通带边缘必须大于0'）
end
ifWs<=Wp
error（'阻带边缘必须大于通带边缘'）
end
if（Rp<=0）|（Ar<0）
error（'通带波动或阻带衰减必须大于0'）
end
N=ceil（log10（（10^（Rp/10）-1）/（10^（Ar/10）-1））/（2*log10（Wp/Ws）））;
fprintf（'\nButterworth滤波器阶数=%2.0f\n',N）
OmegaC=Wp/（（10^（Rp/10）-1）^（1/（2*N）））;
[b,a]=u_buttap（N,OmegaC）;
⑵u_buttap函数
function[b,a]=u_buttap（N,Omegac）
%未归一化巴特沃思模拟滤波器设计函数
%b分子多项式系数；a分母多项式系数；N阶数；Omegac截止频率，单位为弧度秒
[z,p,k]=buttap（N）;
p=p*Omegac;
k=k*Omegac^N;
B=real（poly（z））;
b0=k;
b=k*B;
a=real（poly（p））;
⑶freqs_m函数
function[db,mag,pha,w]=freqs_m（b,a,wmax）
%s域频率响应计算freqs的改进版本
%wmax希望获得的最大频率响应；b直接型滤波器分子系数；a直接型滤波器分母系数
%mag为[0wmax]区间的绝对值；pha为[0wmax]区间的相位响应；w为[0wmax]区间内的500个频率样本数组
w=[0:
1:
500]*wmax/500;
H=freqs（b,a,w）;
mag=abs（H）;
db=20*log10（（mag+eps）/max（mag））;
pha=angle（H）;
⑷sdir2cas函数
function[C,B,A]=sdir2cas（b,a）
%s平面上直接型到级联型的转换
%C增益系数；B级联子系统分子矩阵；A级联子系统分母矩阵；
%b直接型分子多项式系数；a直接型分母多项式系数；
Na=length（a）-1;Nb=length（b）-1;
%计算增益系数C
b0=b
（1）;b=b/b0;
a0=a
（1）;a=a/a0;
C=b0/a0;
%分母的二阶因子部分
p=cplxpair（roots（a））;K=floor（Na/2）;
ifK*2==Na%当Na为偶时计算
A=zeros（K,3）;
forn=1:
2:
Na
Arow=p（n:
1:
n+1,:
）;
Arow=poly（Arow）;
A（fix（（n+1）/2）,:
）=real（Arow）;
end
elseifNa==1%当Na=1时计算
A=[0real（poly（p））];
else%当Na为奇及>1时计算
A=zeros（K+1,3）;
forn=1:
2:
2*K
Arow=p（n:
1:
n+1,:
）;
Arow=poly（Arow）;
A（fix（（n+1）/2）,:
）=real（Arow）;
end
A（K+1,:
）=[0real（poly（p（Na）））];
end
%分子的二阶因子部分
z=cplxpair（roots（b））;K=floor（Nb/2）;
ifNb==0%当Nb=0时计算
B=[00poly（z）];
elseifK*2==Nb%当Nb为偶时计算
B=zeros（K,3）;
forn=1:
2:
Nb
Brow=z（n:
1:
n+1,:
）;
Brow=poly（Brow）;
B（fix（（n+1）/2）,:
）=real（Brow）;
end
elseifNb==1%当Nb=1时计算
B=[0real（poly（z））];
else%当Nb为奇且>1时计算
B=zeros（K+1,3）;
forn=1:
2:
2*K
Brow=z（n:
1:
n+1,:
）;
Brow=poly（Brow）;
B（fix（（n+1）/2）,:
）=real（Brow）;
end
B（K+1,:
）=[0real（poly（z（Nb）））];
end
⒉脉冲响应不变法
利用MATLAB实现脉冲响应不变法，即已知系统函数
，利用residue函数得到系统函数零极点，然后根据
将每一个模拟极点映射为数字极点。
最后利用residuez函数将
转换成滤波器直接形式，这一过程用imp_invr函数实现。
⑴滤波器设计指标为
，用3种方法设计数字滤波器。
程序Example3_2如下：
Wp=0.2*pi;
Ws=0.3*pi;
Rp=1;
Ar=15;
T=1;%数字滤波器设计指标
OmegaP=Wp*T;%模拟原型滤波器通带截止频率
OmegaS=Ws*T;%模拟原型滤波器阻带截止频率
%
[cs1,ds1]=afd_butt（OmegaP,OmegaS,Rp,Ar）;%巴特沃思模拟滤波器设计
[b1,a1]=imp_invr（cs1,ds1,T）;%用脉冲响应不变法实现模拟滤波器到数字滤波器的变换
[db1,mag1,pha1,grd1,w1]=freqz_m（b1,a1）;%计算频率响应和相位响应
[C1,B1,A1]=dir2par（b1,a1）%直接型到并联型转换
%画图
figure
（1）
subplot（341）
plot（w1/pi,mag1）;
axis（[0101.1]）;
ylabel（'（a）'）
title（'幅度响应'）;
subplot（342）;
plot（w1/pi,db1）;
axis（[01-405]）;
title（'幅度响应（dB）'）;
subplot（343）;
plot（w1/pi,pha1/pi）;
axis（[01-11]）;
title（'相位响应'）;
subplot（344）;
plot（w1/pi,grd1）;
axis（[01010]）;
title（'群延迟'）;
%
[cs2,ds2]=afd_chb1（OmegaP,OmegaS,Rp,Ar）;%未归一化切比雪夫Ι型模拟低通滤波器设计
[b2,a2]=imp_invr（cs2,ds2,T）;%用脉冲响应不变法实现模拟滤波器到数字滤波器的变换
[db2,mag2,pha2,grd2,w2]=freqz_m（b2,a2）;%计算频率响应和相位响应
[C2,B2,A2]=dir2par（b2,a2）%直接型到并联型转换
%画图
subplot（345）
plot（w2/pi,mag2）;
axis（[0101.1]）;
ylabel（'（b）'）
subplot（346）;
plot（w2/pi,db2）;
axis（[01-405]）;
subplot（347）;
plot（w2/pi,pha2/pi）;
axis（[01-11
升级会员 