数字信号处理实验三离散时间信号的频域分析.docx

数字信号处理实验三离散时间信号的频域分析.docx

《数字信号处理实验三离散时间信号的频域分析.docx》由会员分享，可在线阅读，更多相关《数字信号处理实验三离散时间信号的频域分析.docx（10页珍藏版）》请在冰豆网上搜索。

数字信号处理实验三离散时间信号的频域分析

实验三:

离散时间信号的频域分析

一．实验目的

1．在学习了离散时间信号的时域分析的基础上，对这些信号在频域上进行分析，从而进一步研究它们的性质.

2．熟悉离散时间序列的3种表示方法:

离散时间傅立叶变换（DTFT），离散傅立叶变换（DFT）和Z变换.

二．实验相关知识准备

1．用到的MATLAB命令

运算符和特殊字符：

〈>。

*^.^

语言构造与调试：

errorfunctionpause

基本函数:

angleconjrem

数据分析和傅立叶变换函数：

fftifftmaxmin

工具箱：

freqzimpzresiduezzplane

三．实验内容

1．离散傅立叶变换

在MATLAB中，使用fft可以很容易地计算有限长序列x［n］的离散傅立叶变换。

此函数有两种形式：

y=fft（x）

y=fft（x,n）求出时域信号x的离散傅立叶变换

n为规定的点数，n的默认值为所给x的长度。

当n取2的整数幂时变换的速度最快。

通常取大于又最靠近x的幂次。

（即一般在使用fft函数前用n=2^nextpow2（length（x））得到最合适的n）。

当x的长度小于n时，fft函数在x的尾部补0，以构成长为n点数据。

当x的长度大于n时，fft函数将序列x截断，取前n点。

一般情况下，fft求出的函数多为复数，可用abs及angle分别求其幅度和相位。

注意：

栅栏效应，截断效应（频谱泄露和谱间干扰），混叠失真

例3－1：

fft函数最通常的应用是计算信号的频谱。

考虑一个由100hz和200hz正弦信号构成的信号,受零均值随机信号的干扰，数据采样频率为1000hz。

通过fft函数来分析其信号频率成分。

t=0：

0.001:

1;%采样周期为0。

001s，即采样频率为1000hz

x=sin（2＊pi＊100＊t）+sin（2＊pi*200*t）+1。

5*rand（1,length（t））;%产生受噪声污染的正弦波信号

subplot（2，1，1）;

plot（x（1:

50））；%画出时域内的信号

y=fft（x，512）;%对x进行512点的fft

f=1000＊（0：

256）/512；％设置频率轴（横轴）坐标,1000为采样频率

subplot（2,1，2）;

plot（f，y（1:

257））；％画出频域内的信号

实验内容3－2：

频谱泄漏和谱间干扰

假设现有含有三种频率成分的信号x（t）=cos（200πt）+sin（100πt）+cos（50πt）

用DFT分析x（t）的频谱结构。

选择不同的截取长度，观察DFT进行频谱分析十存在的截断效应。

试用加窗的方法减少谱间干扰。

请分析截取长度对频谱泄漏和频率分辨率的影响，分析不同窗函数对谱间干扰的影响。

提示：

截断效应使谱分辨率（能分开的两根谱线间的最小间距）降低,并产生谱间干扰；频谱混叠失真使折叠频率（fs/2）附近的频谱产生较大的失真。

理论和实践都已证明，加大截取长度可提高频率分辨率；选择合适的窗函数可降低谱间干扰；而频谱混叠失真要通过提高采样频率fs和预滤波来改善.

解：

取采样频率fs=400Hz，采样信号序列x（n）=x（t）w（n），n=0,1…….N—1；

N为采样点数，N=fs*T，T为截取时间长度，w（n）为窗函数。

实验取三种长度T1=0.04s，T2=4*0.04s，T3=8＊0。

04s

窗函数分别用矩形窗函数w（n）=RN（n），Hamming窗。

clear;closeall

fs=400;T=1/fs；％采样频率为400

Tp=0.04；N=Tp＊fs;％采样点数

N1=［N，4＊N，8＊N];%设定三种截取长度供调用

st=［’|X1（jf）|’;'|X4（jf）|’;’|X8（jf）|’］;%设定三种标注语句供调用

%矩形窗截断

form=1:

3

n=1：

N1（m）;

xn=cos（200＊pi＊n*T）+sin（100*pi＊n*T）+cos（50*pi＊n＊T）;%产生采样序列

Xk=fft（xn,4096）%4096点DFT，用FFT实现

fk=［0:

4095］/4096/T;

subplot（3，2，2＊m—1）

plot（fk，abs（Xk）/max（abs（Xk）））;

ylabel（st（m，:

））；

ifm==1title（'矩形窗截取'）;end

end；

％加hamming窗改善谱间干扰

form=1:

3

n=1：

N1（m）;

wn=hamming（N1（m））;％调用函数hamming产生N长hamming窗序列wn

xn=（cos（200＊pi*n＊T）+sin（100＊pi＊n＊T）+cos（50＊pi*n＊T））.＊wn’；％产生采样序列

Xk=fft（xn，4096）%4096点DFT，用FFT实现

fk=［0：

4095]/4096/T；

subplot（3，2，2＊m）

plot（fk,abs（Xk）/max（abs（Xk）））；

ylabel（st（m,：

））;

ifm==1title（’hamming窗截取’）；end

end;

2．一维逆快速傅立叶变换

y=ifft（x）

y=ifft（x,n）

实验内容3—3：

频域采样定理的验证

（1）产生一个三角波序列x（n）

（2）对M=40，计算x（n）的64点DFT，并图示x（n）和X（k）=DFT［x（n）］，k=0，1，…，63。

（3）对

（2）中所得X（k）在[0,2pi］上进行32点抽样得

X1（k）=X（2k）,k=0，1,…,31。

并图示。

（4）求X1（k）的32点IDFT,并图示。

即x1（n）=IDFT［X1（k）］，k=0，1，…,31

（5）采用周期延拓的方法绘出

的波形，评述

与x（n）的关系，并根据频域采样理论加以解释。

clear；closeall

M=40;N=64;n=0：

M；

xa=0:

floor（M/2）;xb=ceil（M/2）—1:

-1:

0；

xn=［xa，xb];％产生长度为M的三角波序列xn

Xk=fft（xn，64）；%计算xn的64点FFT

X1k=Xk（1：

2:

N）；％隔点抽取Xk得到X1（k）

x1n=ifft（X1k，N/2）;％计算X1k的32点IFFT得到x1（n）

nc=0：

3*N/2；%取97点进行观察

xc=x1n（mod（nc，N/2）+1）;％x1（n）的周期延拓序列

subplot（3，2,1）

stem（n，xn,’.’）;

title（'40点三角波序列x（n）’）；

xlabel（'n'）；

ylabel（’x（n）'）

subplot（3，2，2）

n1=0：

N/2-1；

stem（n1,x1n,’.’）；

title（'32点IDFT［X1（k）]'）;

xlabel（'n'）;

ylabel（’x1（n）'）；

subplot（3，2,3）

k=0：

N—1；

stem（k，abs（Xk）,'。

’）;

title（'64点DFT[x（n）］'）;

xlabel（'k'）;

ylabel（’X（k）'）；

subplot（3,2,4）

k1=0:

N/2-1；

stem（k1，abs（X1k）,’。

’）；

title（’64点DFT［x（n）］’）;

xlabel（’k’）;

ylabel（'X1（k）'）;

subplot（3，2，5）；

stem（nc,xc,'.’）;

title（'x1（n）的周期延拓序列'）；

xlabel（’n’）;

ylabel（'x（mod（n.32））’）；

由于频域在［0,2

]上的采样点数N小于x（n）的长度M，所以产生时域混叠现象，不能由X1（k）恢复序列x（n）。

只有满足N大于等于M时，可优频域采样X1（k）得到原序列x（n）。

这就是频域采样鼎立,请同学们自己编程验证

。

3．线性调频Z变换

离散傅立叶变换（DFT）可以看作信号在Z域上沿单位圆的均匀采样。

但在实际应用中，并非整个单位圆上的频谱都有意义.一些情况下,如对于窄带信号,只希望分析信号所在的一段频带等，采样点的轨迹是一条弧线或圆周。

这种需求，就导致了线性调频Z变换（Chirpz变换）的出现.

Chirpz变换与DFT计算整个频谱的算法不同，它是一种更为灵活的计算频谱的算法，可以用来计算单位圆上任一段曲线的Z变换，作频谱分析时输入的点数和输出的点数可以不相等,从而达到频域“细化”的目的。

y=czt（x,m，w，a）

y=czt（x）

例3-4：

利用Chirpz变换计算滤波器h（且h=fir1（30,125/500,boxcar（31）））在100Hz~200Hz的频率特性,并用图形文件比较CZT函数和FFT函数.

h=fir1（30,125/500，boxcar（31））;

fs=1000;

f1=100;f2=200；

m=1024;

y=fft（h,1024）；

fy=fs＊（0：

1023）/1024；

subplot（2，1,1）;

plot（fy,abs（y））;

axis（［0,500，0，1.5］）;

w=exp（—j＊2＊pi*（f2—f1）/（m＊fs））;

a=exp（j＊2*pi＊f1/fs）；

z=czt（h，m,w，a）;%产生1024个z值

fz=（f2—f1）＊（0：

1023）/1024+f1;subplot（2，1,2）;

plot（fz,abs（z））

4．实验问题回答

（1）完成例3－1的程序并观察信号序列x（t）在时域和频域上分析的特点.

（2）在实验内容3-2和3—3中按要求完成程序,并回答3—2和3—3中提出的问题。

（2）完成例3－4的程序并观察Chirpz变换与fft变换的特点.

四．实验报告要求

1．按照实验内容要求完成相关实验程序,并得出相关的实验结果（包括图形结果）。

2．回答实验中提出的问题.

3．总结本次实验结果，按照实验报告格式要求，书写实验报告。

五．实验设备

PC机，MATLAB软件

附录AMATLAB系统的常用概念

1、命令窗口

在Windows2000下启动MATLAB系统后，Windows2000的工作平台上会弹出一个窗口，如下图所示,这个窗口称为MATLAB的命令窗口。

MATLAB的命令窗口是用户与MATLAB解释器进行通信的工作环境，提示符‘〉>’表示MATLAB解释器正等待用户输入命令。

所有的MATLAB命令、MATLAB函数，以及MATLAB程序都要在这个窗口下运行。

在命令窗口,用户可以发出MATLAB命令。

例如，为了生成一个3＊3的矩阵,可以在提示符下，键入如下的命令:

A=[123;456;789]

方括号命令表示矩阵,空格或逗号将每行的元素分开，而分号将矩阵的各行数值分开。

再键入Enter后，MATLAB将回显如下的矩阵:

A=

123

456

789

为了求该矩阵的逆矩阵，则只要键入命令

？

B=inv（A）;

MATLAB就将计算出相应的结果。

如果不想在命令窗口中显示计算结果，只要如上所示，在该命令后多键入一个分号即可。

此时，MATLAB系统只完成该命令所要求的计算任务，其计算结果不回显.这项功能在程序设计中是非常必要的.

MATLAB系统也可以说是一种新的语言，该语言十分容易掌握，其结构非常类似于数学式子的书写格式，用户花上很少的时间即可掌握MATLAB的大部分命令。

2、图形窗口

MATLAB系统的强大功能之一是其优秀的图形功能。

对于任何作图命令，MATLAB将打开另一个窗口来绘制输出图形，这样的窗口在MATLAB系统中被成为图形窗口。

在同一个图形窗口中，可以绘制多个图形,也可以生成多个图形窗口，并选择其中的一个图形窗口，在其中绘制图形.生成图形窗口的方法比较多，在没有图形窗口存在时,每个绘图函数都能自动生成一个图形窗口;也可以用figure命令生成一个新的图形窗口；还可以用命令窗口的File菜单的New子菜单的Figure项来打开一个新的图形窗口。

3、搜索路径

MATLAB管理着一条搜索路径，它在搜索路径下寻找与命令相关的函数文件。

例如,如果在MATLAB提示符下输入example，MATLAB解释器将按照下面的步骤来处理这条字符串：

（1）检查example是不是一个变量;

（2）如果不是，检查example是不是一个内部函数；

（3）如果不是,检查在当前文件夹下是否存在名为example.mex,example.dll或example。

m的文件。

MEX文件是MATLAB的执行文件,将优先执行；