中南大学数字信号处理实验报告Word下载.doc
《中南大学数字信号处理实验报告Word下载.doc》由会员分享，可在线阅读，更多相关《中南大学数字信号处理实验报告Word下载.doc（21页珍藏版）》请在冰点文库上搜索。
如果对此连续周期信号进行抽样，其抽样时间间隔为T，抽样后信号以表示，则有，如果令为数字频率，满足，其中是抽样重复频率，简称抽样频率。
为了在数字计算机上观察分析各种序列的频域特性，通常对在上进行M点采样来观察分析。
对长度为N的有限长序列x（n）,有
其中
通常应取得大一些，以便观察谱的细节变化。
取模可绘出幅频特性曲线。
（3）用DFT进行普分析的三种误差
三种误差：
混叠现象、泄露现象、栅栏效应
a）混叠现象
当采样频率小于两倍信号（这里指是信号）最大频率时，经过采样就会发生频谱混叠，这使得采样后的信号序列频谱不能真实地反映原信号的频谱。
所以在利用DFT分析连续信号的频谱时，必须注意这一问题。
避免混叠现象的唯一方法是保证采样速率足够高，使频谱交叠现象不致出现。
也就是说，在确定采样频率之前，必须对信号的性质有所了解，一般在采样前，信号通过一个防混叠低通滤波器。
b）泄漏现象
实际中的信号序列往往很长,为了方便我们往往用截短的序列来近似它们,这样可以使用较短的DFT来对信号进行频谱分析,这种截短等价于给原信号序列乘以一个矩形窗函数。
泄漏是不能与混叠完全分离开的,因为泄漏导致频谱的扩散，从而造成混叠。
为了减小泄漏的影响，可以选择适当的窗函数，使频谱的扩散减到最小。
c）栅栏效应
因为DFT是对单位圆上Z变换的均匀采样，所以他不可能将频谱视为一个连续函数。
这样就产生了栅栏效应，就一定意义上看，DFT来观看频谱就好像通过一个尖桩的栅栏来观看一个图景一样，只能在离散点上看到真实频谱，这样就可能发生一些频谱的峰点或谷点被“尖桩的栅栏”所挡住，不能被我们观察到。
减小栅栏效应的一个方法就是借助在原序列的末端添补一些零值，从而变动DFT的点数。
这一方法实际上是人为地改变了对真实谱采样的点数和位置，相当于搬动了每一根“尖桩栅栏”的位置，从而使得频谱的峰点或者谷点暴露出来。
当然，这是每根谱线所对应的频率和原来的不同了。
综上所述，DFT可以用于信号的频谱分析，但必须注意可能产生的误差，在应用过程中要尽可能减少和消除这些误差的影响。
三、实验内容及要求
（1）复习常用离散时间信号的有关内容；
（2）用MATLAB编程产生上述任意3种序列（长度可输入确定，对（d）（e）（f）中的参数可自行选择），并绘出其图形；
（3）混叠现象
对连续信号其中，进行采样，分别取采样频率，观察的变化，并做记录（打印曲线）,观察随着采样频率降低频谱混叠是否明显存在，说明原因。
（4）截断效应
给定，截取一定长度的信号，为窗函数，长度为N，。
做2N点DFT变换，分析当N逐渐增大时，分析是否有频谱泄露现象、主瓣的宽度变化？
如何减小泄露？
（5）栅栏效应
给定，分别计算在频率区间上的16点、32点、64点等间隔采样，绘制采样的幅频特性图，分析栅栏效应，如何减小栅栏效应？
（6）提高题：
给定信号，，。
分别计算在频率区间上的频谱，观察其幅频特性图，分析是否存在谱间干扰，如何减小谱间干扰？
四、实验用MATLAB函数介绍
（1）数字信号处理中常用到的绘图指令（只给出函数名，具体调用格式参看help）
figure（）;
plot（）;
stem（）;
axis（）;
gridon;
title（）;
xlabel（）;
ylabel（）;
text（）;
holdon;
subplot（）
（2）离散时间信号产生可能涉及的函数
zeros（）;
ones（）;
exp（）;
sin（）;
cos（）;
abs（）;
angle（）;
real（）;
imag（）;
五、实验结果与分析
在时间上依次出现的数值序列,例如，{…，0.5，1，2，-1，0，5，…}。
相邻两个数之间的时间间隔可以是相等的，也可以是不等的。
在前一情况下，设时间间隔为T秒,则离散信号可用符号x（nT）来表示。
在间隔T归一化为1的条件下,T可以省略,即将x（nT）表示为x（n）。
x（n）既可表示整个序列,也可表示离散信号在nT瞬间的值。
单位阶跃序列：
n=-20:
20;
xn=heaviside（n）;
xn（n==0）=1;
plot（n,xn）;
stem（n,xn）;
axis（[-202001.2]）;
title（'
单位阶跃序列'
）;
xlabel（'
n'
ylabel（'
u（n）'
boxon
矩阵序列：
N=5;
xn=heaviside（n）-heaviside（n-N）;
xn（n==N）=0;
矩阵序列'
R_{N}（n）'
正弦序列：
n=-40:
40;
A=2;
w=pi/8;
f=pi/4;
xn=A*sin（w.*n+f）;
axis（[-4040-4.24.2]）
正弦序列'
x（n）'
subplot（2,2,2）;
f01=500;
fs=1000;
n=1/fs;
N=length（t）;
t=0:
n:
0.1;
x=sin（2*pi*f01*t）;
plot（t,x）;
X=fft（x,N）;
plot（fs/N*（0:
（N/2-1））,abs（X（1:
（N/2））））;
取样频率为1000'
subplot（2,2,1）;
fs=2000;
N/2-1））,abs（X（1:
N/2）））;
取样频率为2000'
subplot（2,2,3）;
fs=1200;
取样频率为1200'
subplot（2,2,4）;
fs=800;
取样频率为800'
结论：
连续信号的最高频率为，根据取样定理，采样频率必须满足Fs>
=2fc，否则会在折叠频率Fs/2处出现频谱混叠。
通过实验，当Fs为1200HZ，2000HZ时，峰值在500HZ处，没有发生混叠。
但当取样频率为800HZ时，峰值在300HZ处，产生较大的误差。
n=50;
Rn=[ones（1,n）];
wn=Rn;
xn=cos（（pi./4）*（0:
n-1））;
yn=xn.*wn;
N=20;
Y=fft（yn,2*N）;
plot（2*pi/N*（0:
N/2-1）,abs（Y（（1:
N/2））））;
截断效应N=20'
N=40;
截断效应N=40'
N=100;
截断效应N=100'
N=1000;
截断效应N=1000'
由上图可知，随着矩形窗的N增大，主瓣宽度变窄，分辨率提高，泄露也相继减少，峰起值占总比越来越接近9%。
随N减少到20时，即取样长度比较小时，波形泄露比较严重，无法反映实际波形。
另外，为了减小谱间干扰，应用其他形状的窗函数w（n）代替矩形窗会降低泄露程度。
%
n=0:
1:
10;
xn=[ones（1,4）,zeros（1,7）];
%输入时域序列向量xn=R4（n）
Xk16=fft（xn,16）;
%计算xn的16点DFT
Xk32=fft（xn,32）;
%计算xn的32点DFT
Xk64=fft（xn,64）;
%计算xn的64点DFT
stem（n,xn,'
.'
（a）x_1（n）'
x_1（n）'
k=0:
15;
wk=2*k/16;
%产生16点DFT对应的采样点频率（关于π归一化值）
stem（wk,abs（Xk16）,'
%绘制16点DFT的幅频特性图
（c）16点DFT的幅频特性图'
\omega/\pi'
幅度'
31;
wk=2*k/32;
%产生32点DFT对应的采样点频率（关于π归一化值）
stem（wk,abs（Xk32）,'
%绘制32点DFT的幅频特性图
（d）32点DFT的幅频特性图'
63;
wk=2*k/64;
%产生64点DFT对应的采样点频率（关于π归一化值）
stem（wk,abs（Xk64）,'
%绘制64点DFT的幅频特性图
（d）64点DFT的幅频特性图'
由栅栏效应可知，采样点的间隔是看不到的。
所以对于有限长序列，我们在原序列尾部补零，从而增大频域采样点数和采样点位置。
只要采样频率Fs足够高，且采样点数足够多，分辨率就会提高到认为DFT的离散谱近似代表原信号频谱。
首先取N=100时
n=-50:
49;
N=length（n）
x1=5*sin（pi/5*n）;
x2=cos（pi/4*n）;
x3=x1+x2;
subplot（2,2,1）
plot（（0:
N-1）*2*pi/N/pi,abs（fftshift（fft（x1,N））））;
X1频谱'
\omega/pi'
|X（e^{j\omega}）|'
axis（[020400]）
subplot（2,2,2）
N-1）*2*pi/N/pi,abs（fftshift（fft（x2,N））））;
X2频谱'
axis（[020100]）
subplot（2,2,3）
N-1）*2*pi/N/pi,abs（fftshift（fft（x3,N））））;
X3频谱'
由N=100时，主瓣两边形成了一定的旁瓣，引起不同频率分量间的干扰。
下面提高N=2000
N=2000时，
n=-1000:
999;
axis（[0205000]）
axis（[0201000]）
当N=2000，谱间干扰程度继续加强，甚至可以误认为强信号谱的旁瓣是另一个频率的信号谱线，从而造成假信号，引起较大偏差。
减小谱间干扰，是用其他的窗函数代替矩形窗，例如三角窗、汉宁窗、哈明窗等。
实验二数字滤波器的设计
四、实验目的
（1）熟悉用双线性变换法设计IIR数字滤波器的原理与方法；
（2）学会调用MATLAB信号处理工具中滤波器设计函数，设计各种IIR滤波器，学会根据滤波需求确定滤波器指标参数；
（3）掌握用窗函数法设计FIR数字滤波器的原理和方法。
五、实验原理
（1）利用双线性变换设计IIR滤波器，首先要设计出满足指标要求的模拟滤波器的传递函数，然后由通过双线性变换可得所要设计的IIR滤波器的系统函数。
如果给定的指标为数字滤波器的指标，则首先要转换成模拟域指标。
（2）模拟滤波器设计
巴特沃兹滤波器的振幅平方函数为
（1）
其传输函数为
（2）
（3）
首先确定技术指标：
i.通带中允许的最大衰减和通带截止频率；
ii.阻带允许的最小衰减和阻带起始频率。
由式（8-11）可得：
（4）
（5）
得到
（6）
（7）
再利用上面两式得到
令
，
则
（8）
已知，可由式（8）求出滤波器的阶数N。
求出的N可能有小数部分一般取大于等于N的最小整数。
关于3dB截止频率，有时在技术指标中给出，如果没有给出可以按照式（6）或式（7）求出。
根据以上所述，巴特沃兹滤波器的设计步骤为：
i.根据要求，由式（8）求出阶数N；
ii.由式（6）或式（7）求出3dB截止频率；
iii.由式（3）求出N个极点；
iv.由式
（2）写出传递函数。
实际设计中，第ⅲ、ⅳ两步由以下两步代替：
v.由N可查下表，得归一化低通巴特沃兹滤波器，
vi.去归一化，即将用代替，得到实际。
（3）设所希望得到的滤波器的理想频率响应为。
那么FIR滤波器的设计就在于寻找一个传递函数去逼近。
在这种逼近中最直接的一种方法是从单位取样响应序列着手，使逼近理想的单位取样响应。
我们知道可以从理想频率响应通过傅里叶反变换来得到，即：
（9）
但是一般来说，这样得到的单位取样响应往往都是无限长序列；
而且是非因果的。
以一个截止频率为的线性相应位理想低通为例来说明。
设低通滤波器的时延为，即：
（10）
这是一个以为中心的偶对称的无限长非因果序列。
这样一个无限长的序列怎样用一个有限长序列去近似呢？
最简单的办法就是直接截取它的一段来代替它。
例如把到的一段截取来作为，但是为要保证所得到的是线性相位滤波器。
必须满足的对称性，所以时延应该取长度的一半，即
这种直接截取的办法可以形象地想象为，好比是通过一个“窗口”所看到的一段。
中表达为和一个“窗口函数”的乘积。
在这里，窗口函数就是矩形脉冲函数，即
但是一般来说，窗口函数并不一定是矩形函数，可以在矩形以内还对作一定的加权处理，因此，一般可以表示为
这里就是窗口函数。
这种对理想单位取样响应加窗的处理对频率响应会产生以下三点影响：
a）使理想特性不连续的边沿加宽，形成一过渡带，过渡带的宽度取决于窗口频谱的主瓣宽度。
b）在过渡带两旁产生肩峰和余振，它们取决于窗口频谱的旁瓣；
旁瓣越多，余振也越多；
旁瓣相对值越大，肩峰则越强。
c）增加截取长度N，只能缩小窗口频谱的主瓣宽度而不能改变旁瓣的相对值；
旁瓣与主瓣的相对关系只决定于窗口函灵敏的形状。
因此增加N，只能相对应减小过渡带宽。
而不能改变肩峰值。
肩峰值的大小直接决定通带内的平稳和阻带的衰减，对滤波器性能有很大关系。
例如矩形窗的情况下，肩峰达8.95%，致使阻带最小衰减只有21分贝，这在工程上往往是不够的。
怎样才能改善阻带的衰减特性呢？
只能从改善窗口函数的形状上找出路，所以希望的窗口频谱中应该减少旁瓣，使能量集中在主瓣，这样可以减少肩峰和余振，提高阻带的衰减。
而且要求主瓣宽度尽量窄，以获得较陡的过渡带，然而这两个要求总不能兼得，往往需要用增加主瓣宽度带换取决瓣的抑制，于是提出了海明窗、凯塞-贝塞尔窗、切比雪夫窗等窗口函数。
六、实验内容及要求
（1）用双线性变换法设计Butterworth数字滤波器，绘制幅频响应特性曲线。
a）Butterworth低通数字滤波器：
要求其通带截至频率100Hz，阻带截至频率200Hz，通带衰减小于2dB，阻带衰减大于15dB，采样频率Fs=500HZ。
实验代码：
Rp=2;
Rs=15;
Fs=500;
Ts=1/Fs;
%双线性变换时频率的预畸变
wp=100*2*pi*Ts;
%利用公式，求数字域通带截止频率wp
ws=200*2*pi*Ts;
%利用公式，求数字域通带截止频率ws
wp1=2*tan（wp/2）/Ts;
%利用公式，进行预畸变
ws1=2*tan（ws/2）/Ts;
[N,Wn]=buttord（wp1,ws1,Rp,Rs,'
s'
%选择滤波器的最小阶数N
[Z,P,K]=buttap（N）;
%计算Butterworth低通原型模拟滤波器
[Bap,Aap]=zp2tf（Z,P,K）;
%零极点增益到传递函数模型的转变
[b,a]=lp2lp（Bap,Aap,Wn）;
%把模拟滤波器原型转换成截止频率为Wn的低通滤波器
[bz,az]=bilinear（b,a,Fs）;
%用双线性变换法实现模拟滤波器到数字滤波器的转换
[H,W]=freqz（bz,az）;
plot（W/pi,abs（H））;
grid;
频率w/pi'
频率相应幅度/dB'
title（
升级会员 