实验5抽样定理.docx
《实验5抽样定理.docx》由会员分享，可在线阅读，更多相关《实验5抽样定理.docx（26页珍藏版）》请在冰豆网上搜索。
实验5抽样定理
实验5抽样定理
一、实验目的：
1、了解用MATLAB语言进行时域、频域抽样及信号重建的方法。
2、进一步加深对时域、频域抽样定理的基本原理的理解。
3、观察信号抽样与恢复的图形，掌握采样频率的确定方法和内插公式的编程方法。
二、实验原理：
1、时域抽样与信号的重建
（1）对连续信号进行采样
例5-1已知一个连续时间信号
，取最高有限带宽频率fm=5f0，分别显示原连续时间信号波形和Fs>2fm、Fs=2fm、Fs<2fm三情况下抽样信号的波形。
程序清单如下：
%分别取Fs=fm，Fs=2fm，Fs=3fm来研究问题
dt=0.1;f0=1;T0=1/f0;
fm=5*f0;Tm=1/fm;
t=-2:
dt:
2;
f=sin（2*pi*f0*t）+1/3*sin（6*pi*f0*t）;
subplot（4,1,1）;plot（t,f）;
axis（[min（t）,max（t）,1.1*min（f）,1.1*max（f）]）;
title（'原连续信号和抽样信号'）;
fori=1:
3;
fs=i*fm;Ts=1/fs;
n=-2:
Ts:
2;
f=sin（2*pi*f0*n）+1/3*sin（6*pi*f0*n）;
subplot（4,1,i+1）;stem（n,f,'filled'）;
axis（[min（n）,max（n）,1.1*min（f）,1.1*max（f）]）;
end
程序运行结果如图5-1所示。
图5-1
（2）连续信号和抽样信号的频谱
由理论分析可知，信号的频谱图可以很直观地反映出抽样信号能否恢复原模拟信号。
因此，我们对上述三种情况下的时域信号求幅度谱，来进一步分析和验证时域抽样定理。
例5-2编程求解例5-1中连续信号及其三种抽样频率（Fs>2fm、Fs=2fm、Fs<2fm）下的抽样信号的幅度谱。
程序清单如下：
dt=0.1;f0=1;T0=1/f0;
fm=5*f0;Tm=1/fm;
t=-2:
dt:
2;
N=length（t）;
f=sin（2*pi*f0*t）+1/3*sin（6*pi*f0*t）;
wm=2*pi*fm;
k=0:
N-1;
w1=k*wm/N;
F1=f*exp（-j*t'*w1）*dt;
subplot（4,1,1）;plot（w1/（2*pi）,abs（F1））;
axis（[0,max（4*fm）,1.1*min（abs（F1））,1.1*max（abs（F1））]）;
fori=1:
3;
ifi<=2c=0;elsec=1;end
fs=（i+c）*fm;Ts=1/fs;
n=-2:
Ts:
2;
N=length（n）;
f=sin（2*pi*f0*n）+1/3*sin（6*pi*f0*n）;
wm=2*pi*fs;
k=0:
N-1;
w=k*wm/N;
F=f*exp（-j*n'*w）*Ts;
subplot（4,1,i+1）;plot（w/（2*pi）,abs（F））;
axis（[0,max（4*fm）,1.1*min（abs（F））,1.1*max（abs（F））]）;
end
程序运行结果如图5-2所示。
由图可见，当满足Fs≥2fm条件时，抽样信号的频谱没有混叠现象；当不满足Fs≥2fm条件时，抽样信号的频谱发生了混叠，即图5-2的第二行Fs＜2fm的频谱图，，在fm=5f0的范围内，频谱出现了镜像对称的部分。
图5-2
（3）由内插公式重建信号
信号重建一般采用两种方法：
一是用时域信号与理想滤波器系统的单位冲激响应进行卷积积分；二是用低通滤波器对信号进行滤波。
本实验只讨论第一种方法。
由理论分析可知，理想低通滤波器的单位冲激响应为
抽样信号
通过滤波器输出，其结果应为
与h（t）的卷积积分：
该式称为内插公式。
由式可见，xa（t）信号可以由其抽样值xa（nT）及内插函数重构。
MATLAB中提供了sinc函数，可以很方便地使用内插公式。
例5-3用上面推导出的内插公式重建例5-1给定的信号。
程序清单如下：
dt=0.01;f0=1;T0=1/f0;
fm=5*f0;Tm=1/fm;
t=0:
dt:
3*T0;
x=sin（2*pi*f0*t）+1/3*sin（6*pi*f0*t）;
subplot（4,1,1）;plot（t,x）;
axis（[min（t）,max（t）,1.1*min（x）,1.1*max（x）]）;
title（'用时域卷积重建抽样信号'）;
fori=1:
3;
fs=i*fm;Ts=1/fs;
n=0:
（3*T0）/Ts;
t1=0:
Ts:
3*T0;
x1=sin（2*pi*f0*n/fs）+1/3*sin（6*pi*f0*n/fs）;
T_N=ones（length（n）,1）*t1-n'*Ts*ones（1,length（t1））;
xa=x1*sinc（fs*pi*T_N）;
subplot（4,1,i+1）;plot（t1,xa）;
axis（[min（t1）,max（t1）,1.1*min（xa）,1.1*max（xa）]）;
end
程序运行结果如图5-3所示：
图5-3
2、频域抽样与信号恢复
（1）频域抽样定理
从理论学习可知，在单位圆上对任意序列的Z变换等间隔采样N点得到：
k=0,1,…,N-1
该式实现了序列在频域的抽样。
那么由频域的抽样得到的频谱的序列能否不失真地恢复原时域信号呢？
由理论学习又知，频域抽样定理由下列公式表述：
表明对一个频谱采样后经IDFT生成的周期序列
是原非周期序列x（n）的周期延拓序列，其时域周期等于频域抽样点数N。
假定有限长序列x（n）的长度为M，频域抽样点数为N，原时域信号不失真地由频域抽样恢复的条件如下：
①如果x（n）不是有限长序列，则必然造成混叠现象，产生误差；
②如果x（n）是有限长序列，且频域抽样点数N小于序列长度M（即N中不能无失真地恢复出原信号x（n）。
③如果x（n）是有限长序列，且频域抽样点数N大于或等于序列长度M（即N≥M），则从
中能无失真地恢复出原信号x（n），即
（2）从频谱抽样恢复离散时间序列
例5-4已知一个时间序列的频谱为
用IFFT计算并求出其时间序列x（n），并绘图显示时间序列。
分析：
该题使用了数字频率，没有给出采样周期，则默认Ts=1S,另外，从
的解析式可以直接看出时域序列xn=[3,2,1,2,3]。
但为说明问题，仍编写程序求解如下：
程序清单如下：
Ts=1;N0=[3,5,10];
forr=1:
3;
N=N0（r）;
D=2*pi/（Ts*N）;
kn=floor（-（N-1）/2:
-1/2）;
kp=floor（0:
（N-1）/2）;
w=[kp,kn]*D;
X=3+2*exp（-j*w）+1*exp（-j*2*w）+2*exp（-j*3*w）+3*exp（-j*4*w）;
n=0:
N-1;
x=ifft（X,N）
subplot（1,3,r）;stem（n*Ts,abs（x））;
box
end
程序运行结果如图5-4所示：
图5-4
注意：
程序中数字频率的排序进行了处理，这是因为
的排列顺序是从0开始，而不是从-（N-1）/2开始。
程序运行后将显示数据：
x=5.00005.00001.0000
x=3.00002.00001.00002.00003.0000
x=3.0000-0.0000i2.0000+0.0000i1.0000-0.0000i2.0000+0.0000i3.0000-0.0000-0.0000+0.0000i0-0.0000i-0.0000+0.0000i0.0000-0.0000i-0.0000+0.0000i
由
的频谱表达式可知，有限长时间序列x（n）的长度M=5，现分别取频域抽样点数为N=3,5,10，由图5-4显示的结果可以验证：
①当N=5和N=10时，N≥M，能够不失真地恢复出原信号x（n）；
②当N=3时，N＜M，时间序列有泄漏，形成了混叠，不能无失真地恢复出原信号x（n）。
混叠的原因是上一周期的后2点与本周期的前两点发生重叠，如下所示：
32123
32123
例5-5已知一个频率范围在[-62.8,62.8]rad/s间的频谱
，用IFFT计算并求出时间序列x（n），用图形显示时间序列。
分析：
本题给出了模拟频率Ω，其中Ωm=62.8，需将其归一化为数字频率。
根据奈奎斯特定理可知，（1/Ts）=Fs≥（2Ωm/2П）,可以推导出Ts≤（П/Ωm）,取Ts=0.05s，即采样频率Fs为20Hz或40П。
程序清单如下：
wm=62.8;Ts=pi/wm;
N0=[8,20];
forr=1:
2
N=N0（r）;
D=2*pi/（Ts*N）;
k=[0:
N-1]+eps;
omg=k*D;
X=sin（0.275*omg）./sin（0.025*omg）;
n=0:
N-1;
x=abs（ifft（X,N））;
subplot（1,2,r）;stem（n*Ts,abs（x））;
box
end
程序运行结果如图5-5所示：
图5-5
由N=20的结果可知，时间序列x（n）是一个矩形窗。
根据DFT的循环移位性质可知，非零数据存在于n=-5:
5的区域，有限长序列的长度为11。
而N=8小于有限长序列的长度，其结果发生了混叠，不能无失真地恢复出原信号x（n）。
3、从频谱恢复连续时间信号
实际应用中，离散信号往往来源于对连续信号的采样，因此，这里要讨论从频谱如何计算连续时间信号。
从本质上讲，计算机处理的都是离散信号。
当使用计算机处理连续信号时，实际上是用采样周期极小的序列信号来近似为连续信号。
因此在处理时，原来对于离散序列处理的理论依然有效。
（1）选择一个符合奈奎斯特定理的很小的采样周期T，将主要的模拟频谱限制在奈奎斯特频率范围内，
；
（2）在[
]的频率区间取N个频率点Ωk，求出对应的数字频谱：
（3）对
坐IDFT，求xa（t）。
假定没有发生时间混叠，则
（4）作图。
用plot自动进行插值，获得连续信号。
例5-6已知图5-6所示的理想低通滤波器的模拟频率Ωc=3，在|Ω|≤Ωc范围内幅度为1，|Ω|＞Ωc时幅度为0。
要求计算连续脉冲响应xa（t）。
图5-6
分析：
由奈奎斯特定理可知采样频率Ωs≥2Ωc,即采样周期Ts≤П/Ωc恢复原信号时不会发生混叠。
选的再小一些可以增加样点数，因此可以选Ts=0.1П/Ωc=0.1047s。
同时，为使时间信号尽量接近连续信号，需提高N点的个数。
可以由模拟频率的分辨率公式推导：
D=2П/NTs≤0.1Ωc，使频率分辨率小于有效带宽的1/10，得到：
N≥20П/ΩcTs
程序清单如下：
wc=3;Tmax=0.1*pi/wc;
Ts=input（'（TsNmin=20*pi/wc/Ts;
N=input（'（N>Nmin）N='）;
D=2*pi/（Ts*N）;
M=floor（wc/D）;
Xa=[ones（1,M+1）,zeros（1,N-2*M-1）,ones（1,M）];
n=-（N-1）/2:
（N-1）/2;
xa=abs（fftshift（ifft（Xa/Ts）））;
plot（n*Ts,xa）;
程序执行过程中，在MATLAB命令窗口将给出提示：
输入Ts和N的值，再给出绘图结果。
图5-7是分别输入Ts=0.1s，N=300和Ts=0.1s，N=1000两组数据的运行结果。
图5-7
三、实验内容：
1、阅读并输入实验原理中介绍的例题程序，观察输出的数据和图形，结合基本原理理解每一条语句的含义。
2、已知一个连续时间信号f（t）=sinc（t），取最高有限带宽频率fm=1Hz。
（1）分别显示原连续信号波形和Fs=fm、Fs=2fm、Fs=3fm三种情况下抽样信号的波形；
（2）求解原连续信号和抽样信号的幅度谱；
（3）用时域卷积的方法（内插公式）重建信号。
3、已知一个时间序列的频谱为：
分别取频域抽样点数N为3、5和10，用IFFT计算并求出其时间序列x（n），绘图显示个时间序列。
由此讨论由频域抽样不失真地恢复原时域信号的条件。
4、已知一个频率范围在[-6.28,6.28]rad/s的频谱，在模拟频率|Ωc|=3.14处幅度为1，其它范围幅度为0.计算其连续信号xa（t），并绘图显示信号曲线。
四、实验预习：
1、认真阅读实验原理部分，明确实验目的，读懂例题程序，了解实验方法。
2、根据实验任务预先编写实验程序。
3、预习思考题：
①什么是内插公式？
在MATLAB中内插公式可用什么函数来编写？
②从频域抽样序列不失真地恢复离散时域信号的条件是什么？
五、实验报告：
1、列写调试通过的实验程序，打印实验程序产生的曲线图形。
2、给出预习思考题答案。
3、思考题：
①试归纳用IFFT数值计算方法从频谱恢复离散时间序列的方法和步骤。
②从频谱恢复连续时间信号与恢复离散时间序列有何不同？
实验6数字滤波器的网络结构
一、实验目的
1、加深对数字滤波器分类与结构的了解。
2、明确数字滤波器的基本结构及其相互间的转换方法。
3、掌握MATLAB语言进行数字滤波器各种结构相互间转换的子函数及程序编写方法。
二、实验原理、
1、数字滤波器的分类
离散LSI系统对信号的响应过程实际上就是对信号进行滤波的过程。
因此，离散LSI系统又称为数字滤波器。
数字滤波器从滤波功能上可以分为低通、高通、带通、带阻以及全通滤波器；根据单位脉冲响应的特性，又可以分为有限长单位脉冲响应滤波器（FIR）和无限长单位脉冲响应滤波器（IIR）。
一个离散LSI系统可以用系统函数来表示：
也可以用差分方程来表示：
以上两个公式中，当ak至少有一个不为0时，则在有限Z平面上存在极点，表达的是以一个IIR数字滤波器；当ak全都为0时，系统不存在极点，表达的是一个FIR数字滤波器。
FIR数字滤波器可以看成是IIR数字滤波器的ak全都为0时的一个特例。
IIR数字滤波器的基本结构分为直接Ⅰ型、直接Ⅱ型、直接Ⅲ型、级联型和并联型。
FIR数字滤波器的基本结构分为横截型（又称直接型或卷积型）、级联型、线性相位型及频率采样型等。
本实验对线性相位型及频率采样型不做讨论，见实验8。
另外，滤波器的一种新型结构——格型结构也逐步投入应用，有全零点FIR系统格型结构、全极点IIR系统格型结构以及全零极点IIR系统格型结构。
2、IIR数字滤波器的基本结构与实现
（1）直接型与级联型、并联型的转换
例6-1已知一个系统的传递函数为
将其从直接型（其信号流图如图6-1所示）转换为级联型和并联型。
图6-1
分析：
从直接型转换为级联型，就是将系统的传递函数（tf）模型转换为二次分式（sos）模型；从直接型转换为并联型，就是将系统函数的传递函数（tf）模型转换为极点留数（rpk）模型。
程序清单如下：
b=[8,-4,11,-2];
a=[1,-1.25,0.75,-0.125];
[sos,g]=tf2sos（b,a）
[r,p,k]=residuez（b,a）
运行结果如下：
sos=
1.0000-0.190001.0000-0.25000
1.0000-0.31001.31611.0000-1.00000.5000
g=8
r=
-8.0000-12.0000i
-8.0000+12.0000i
8.0000
p=
0.5000+0.5000i
0.5000-0.5000i
0.2500
k=16
由sos和g的数据，可以列写出级联型的表达式：
信号流图如图6-2所示：
图6-2
由r、p、k的数据，可以列写出并联型的表达式：
上式中出现了复系数，可采用二阶分割将共轭极点组成分母上的实习顺二阶环节。
这里使用自定义函数dir2par可以实现滤波器结构从直接型向并联型的转换，且用实系数二阶环节表示。
在使用该函数时，调用了另一个自定义函数cplxcomp以进行复共轭对的正确排序，保证系统二阶环节的分子、分母一定是实数。
dir2par函数和cplxcomp函数定义如下：
functionI=cplxcomp（p1,p2）%按共轭条件排列极点留数对
%比较两个包含同样标量元素但（可能）具有不同下标的复数对
%本语句必须用在p2=cplxpair（p1）语句之后，以重新排序对应的留数向量
I=[];
forj=1:
length（p2）
fori=1:
length（p1）
if（abs（p1（i）-p2（j））<0.0001）
I=[I,i];
end
end
end
I=I';
function[C,B,A]=dir2par（num,den）%直接型到并联型的转换
M=length（num）;N=length（den）;
[r1,p1,C]=residuez（num,den）;%先求系统的单根p1对应的留数r1及直接项C
p=cplxpair（p1,10000000*eps）;
I=cplxcomp（p1,p）;
r=r1（I）;
K=floor（N/2）;B=zeros（K,2）;A=zeros（K,3）;
ifK*2==N;
fori=1:
2:
N-2;
Brow=r（i:
1:
i+1,:
）;
Arow=p（i:
1:
i+1,:
）;
[Brow,Arow]=residuez（Brow,Arow,[]）;
B（fix（（i+1）/2）,:
）=real（Brow）;
A（fix（（i+1）/2）,:
）=real（Arow）;
end
[Brow,Arow]=residuez（r（N-1）,p（N-1）,[]）;
B（K,:
）=[real（Brow）,0];A（K,:
）=[real（Arow）,0];
else
fori=1:
2:
N-1;
Brow=r（i:
1:
i+1,:
）;
Arow=p（i:
1:
i+1,:
）;
[Brow,Arow]=residuez（Brow,Arow,[]）;
B（fix（（i+1）/2）,:
）=real（Brow）;
A（fix（（i+1）/2）,:
）=real（Arow）;
end
end
将例6-1从直接型转换为并联型的程序改写如下：
b=[8,-4,11,-2];
a=[1,-1.25,0.75,-0.125];
[C,B,A]=dir2par（b,a）
运行结果如下：
C=16
B=
-16.000020.0000
8.00000
A=
1.0000-1.00000.5000
1.0000-0.25000
由A,B,C的数据可以直接写出并联型的表达式：
信号流图如图6-3所示：
图6-3
例6-2已知一个系统的级联型系数公式为
将其从级联型（信号流图如图6-4所示）转换为直接型和并联型结构。
图6-4
分析：
从级联型转换为直接型，就是将二次分式（sos）模型转换为系统传递函数（tf）模型；再使用dirpar.m和cplxcomp.m函数将直接型转换为并联型。
程序清单如下：
sos=[10.901-0.250
1-32110.5];
g=0.5;
[b,a]=sos2tf（sos,g）
[C,B,A]=dir2par（b,a）
程序运行结果如下：
b=0.5000-1.0500-0.35000.9000
a=1.00000.75000.2500-0.1250
C=-7.2000
B=3.98461.6308
3.71540
A=1.00001.00000.5000
1.0000-0.25000
由b,a的数据可以直接写出直接型的表达式：
信号流图如图6-5所示：
图6-5
由A,B,C的数据可以写出并联型的表达式：
信号流图如图6-6所示：
图6-
（2）直接型转换为全零极点IIR系统的格型结构
例6-3将例6-1给定的系统传递函数
从直接型转换为格型。
程序清单如下：
b=[8,-4,11,-2];
a=[1,-1.25,0.75,-0.125];
[K,C]=tf2latc（b,a）
[b,a]=latc2tf（K,C）
程序运行结果如下：
K=
-0.7327
0.6032
-0.1250
C=
8.1064
7.4841
8.5000
-2.0000
b=8.0000-4.000011.0000-2.0000
a=1.0000-1.25000.7500-0.1250
由K、C参数可以画出格型结构图，如同图6-7所示：
图6-7
（3）直接型转换为全极点IIR系统的格型结构
例6-4将一个全极点IIR系统的传递函数
从直接型转换为格型结构。
程序清单如下：
b=[1];
a=[1,-1.25,0.75,-0.125];
K=tf2latc（b,a）
[b,a]=latc2tf（K）
程序运行结果如下：
K=
-0.7327
0.6032
-0.1250
b=1.0000-1.25000.7500-0.1250
a=1
格型结构信号流图如图6-8所示：
图6-8
3、FIR数字滤波器的基本结构与实现
（1）横截型与级联型之间的转换
例6-5已知一个FIR系统的传递函数为
将其从横截型（信号流图如图6-9所示）转换为级联型。
图6-9
分析：
从横截型转换为级联型就是将系统传递函数（tf）模型转换为二次分式（sos）模型。
程序清单如下：
b=[2,0.9,1.55,2.375];
a=[1];
[sos,g]=tf2sos（b,a）
[b,a]=sos2tf（sos,g）
程序运行结果如下：
sos=
1.00000.950001.000000
1.0000-0.50001.25001.000000
g=2
b=2.00000.90001.55002.3750
a=1000
由sos和g的数据可以写出级联型的表达式：
信号流图如图6-10所示：
图6-10
（2）横截型转换为全零点FIR系统的格型结构
例6-6已知一个FIR系统的传递函数为
将其从横截型转换为全零点FIR系统的格型结构
程序清单如下：
b=[1,2.7917,2,1.375,0.3333];
a=[1];
K=tf2latc（b,a）
[b,a]=latc2tf（K）
程序运行结果如下：
K=
2.0004
0.2498
0.5001
升级会员 