数字信号处理实验三离散时间信号的频域分析.docx
《数字信号处理实验三离散时间信号的频域分析.docx》由会员分享，可在线阅读，更多相关《数字信号处理实验三离散时间信号的频域分析.docx（16页珍藏版）》请在冰豆网上搜索。
数字信号处理实验三离散时间信号的频域分析
实验三：
离散时间信号的频域分析
一．实验目的
1．在学习了离散时间信号的时域分析的基础上，对这些信号在频域上进行分析，从而进一步研究它们的性质。
2．熟悉离散时间序列的3种表示方法：
离散时间傅立叶变换（DTFT）,离散傅立叶变换（DFT）和Z变换。
二．实验相关知识准备
1．用到的MATLAB命令
运算符和特殊字符：
< > .* ^ .^
语言构造与调试：
error function pause
基本函数：
angle conj rem
数据分析和傅立叶变换函数：
fft ifft max min
工具箱：
freqz impz residuez zplane
三．实验内容
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软件
附录A MATLAB系统的常用概念
1、命令窗口
在Windows2000下启动MATLAB系统后，Windows2000的工作平台上会弹出一个窗口，如下图所示，这个窗口称为MATLAB的命令窗口。
MATLAB的命令窗口是用户与MATLAB解释器进行通信的工作环境，提示符‘>>’表示MATLAB解释器正等待用户输入命令。
所有的MATLAB命令、MATLAB函数，以及MATLAB程序都要在这个窗口下运行。
在命令窗口，用户可以发出MATLAB命令。
例如，为了生成一个3*3的矩阵，可以在提示符下，键入如下的命令：
A=[123;456;789]
方括号命令表示矩阵，空格或逗号将每行的元素分开，而分号将矩阵的各行数值分开。
再键入Enter后，MATLAB将回显如下的矩阵：
A=
1 2 3
4 5 6
7 8 9
为了求该矩阵的逆矩阵，则只要键入命令
？
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的执行文件，将优先执行；
（4）如果不存在，检查在MATLAB的搜索路径的目录下是否存在名为MEX,example.mex,example.dll,或example.m的文件。
MEX文件优先执行。
用户可以打开路径浏览器（PathBrowse）查看MATLAB系统的当前搜索路径，也可以在其中加入自己的路径。
4、文件类型
在MATLAB系统中，根据功能可将MATLAB系统所使用的外部文件分为几类，并用不同的扩展名作为其标识，我们用的主要是M文件。
M文件以字母m为其扩展名，例如startup.m。
一般说来，M文件是以ASCII码文本文件，可以用任何文本编辑器进行编辑。
在MATLAB系统中，有两类M文件。
一类称为程序M文件，简称M文件；另一类称为函数M文件，或简称为函数，统称为M文件。
M文件的内容是由符合MATLAB语法的语句构成的，函数M文件的第一行必须是以关键字function开始的函数说明语句。
两类M文件的共同特征是：
在MATLAB命令窗口中的命令提示符下键入文件名，来执行M文件中的所有语句规定的计算任务或完成一定的功能。
它们的区别在于以下两方面：
第一，程序M文件中创建的变量都是MATLAB工作空间中的变量，工作空间中的其他程序或函数可以共享，而函数M文件中创建的所有变量除了全程变量外，均局限于函数运行空间内的局部变量；第二，函数M文件的调用式中可以有输入参数和输出参数，而程序M文件则没有这种功能。
5、语言语法要素
MATLAB只管理一种对象——矩阵。
可以使用下列的任何一种方法在MATLAB环境下创建或输入一个矩阵：
1）显示的输入一个元素序列；
2）用MATLAB的内部函数创建一个矩阵；
3）在M文件中用MATLAB语句创建一个矩阵；
4）从一个外部数据文件中装载并创建一个矩阵。
在MATLAB中有两个基本概念：
变量和表达式。
变量由变量名表示，函数名作为特殊的变量名看待，每个变量名由一个字母后面跟随任意个字母或数字（包括下划线）组成，但MATLAB只能分辨前19个字符。
MATLAB能区分组成变量名的大小写字母。
MATLAB的语句则是两种形式之一：
变量名=表达式或者表达式。
在前一种语句形式下，MATLAB将运算的结果赋给“变量名”；而在第二种语句形式下，将运算的结果赋给MATLAB的永久变量ans，每条语句以回车符结束。
一般的，运算的结果在命令窗口中显示出来。
如果语句的最后一个字符是分号“；”，那末，MATLAB仅执行赋值运算，不再显示运算的结果。
与C语言一样，MATLAB将字符串当作数组或矩阵处理。
6、MATLAB的基本运算符
矩阵运算符
矩阵A的转置。
A+B,A-B
矩阵A与B的和与差。
A*B
矩阵A与B的乘法。
A.*B
矩阵A与B的对应元素相乘
关系运算符
<
小于
>
小于或等于
<=
大于
>=
大于或等于
==
等于
~=
不等于
7、特殊运算符
在MATLAB的M文件中，可以加入解释行。
解释行的标识符为“%”，该标识符将被作为注解内容。
程序执行时，注解被忽略。
方括号“[]”用于生成矩阵。
特别的，语句A=[]生成空矩阵A。
行分隔符“；”用于MATLAB语句后时，表示该语句的执行结果不被回显出来，这可避免显示一些不感性趣的结果。
冒号“：
“最主要的作用是生成向量，从下面的例子中可以看出它的使用方法：
j:
k 生成向量[j,j+1,j+2,…,k]
j:
i:
k 生成向量[j,j+i,j+2*i,…,k]，如果j>k，则生成空矩阵
A（:
j） 矩阵A的第j列
A（I,:
） 矩阵A的第I行
A（j:
k） 向量A（j）,A（j+1）,…,A（k）
A（:
j:
k） 从第j列到第k列的矩阵子块
换行连接符“…”，有时一条MATLAB语句会很长，在命令窗口的一行内很可能写不下，此时只要在该语句中加入三连点，再回车即可在下一行接着写该语句。
8、MATLAB的常用数学函数
三角函数
sin
正弦函数
Cos
余弦函数
Tan
正切函数
Asin
反正弦函数
Atan
反正切函数
Sinh
双曲正弦函数
Cosh
双曲余弦函数
Tanh
双曲正切函数
Asinh
反双曲正弦函数
Acosh
反双曲余弦函数
Atanh
反双曲正切函数
Acos
反余弦函数
初等函数
Abs
实数的绝对值、复数的模、字符串的ASIIC值
Angle
复数的幅角
Sqrt
方根函数
Real
复数的实部
Imag
复数的虚部
Conj
复共轭运算
Round
最邻近整数截断（四舍五入）
Ceil
不大于自变量的最大整数
Rem
不小于自变量的最小整数
Exp
自然指数函数（以e为底）
Log
自然对数函数（以e为底）
log10
以10为底的对数函数
9、程序流控制
与其他的程序设计语言一样，MATLAB语言也提供了条件语句。
下面分别予以介绍。
1）for循环语句
MATLAB也有自己的for循环语句。
如果要反复执行的一组语句的循环次数是已知或预定义的，就可以用for循环语句。
例如：
forI=1:
n
x（I）=0;
end
这条语句将向量x的前n个元素赋予零值，这里的变量n必须预先给定。
注意：
每一个for必须与end配对使用。
2）while循环语句
MATLAB提供有while循环语句，它的作用是在一定的逻辑条件的控制下，不断的循环执行一条或一组语句，直到逻辑条件不满足为止。
While语句的一般形式是
While 表达式
语句组
end
3）if条件语句和bread语句
break语句用于退出循环体，if条件语句有两种形式，分别是
if 表达式
语句组1
else
语句组2
end
和
if 表达式1
语句组1
elseif 表达式2
语句组2
else
语句组3
end
10、MATLAB的在线帮助
用户可以随时利用MATLAB的在线帮助查询自己不懂得用法的函数的具体用法，例如：
在命令窗口键入helpabs后的显示如下：
ABSAbsolute value ABS（X） isthe absolutevalueoftheelementsofX.WhenXiscomplex,ABS（X）isthecomplexmodulus（magnitude）of the elementsofX.
SeealsoSIGN,ANGLE,UNWRAP.
Overloadedmethods
Helpsym/abs.m
将abs函数的主要用法和用途都列了出来。
11实验上机的具体过程如下：
1）在windows2000/XP的桌面上找到MATLAB的图标单击进入MATLAB的命令窗口，或在开始菜单里选择程序再找到MATLAB单击也可进入MATLAB的命令窗口。
2）进入命令窗口后，在菜单File中选择open，打开已存在的文件。
如果是新文件，在菜单File中选择New，再找到M-file即可打开MATLAB的编辑窗口，在编辑窗口内输入你的源程序后存盘。
3）
为了使程序能够在MATLAB中运行，需要在搜索路径中加入你的路径，加入路径的过程如下，在MATLAB的命令窗口中选择File菜单的Setpath选项，则打开一个PathBrowser窗口，如下图所示：
在该窗口中单击Browse，选择你的路径加入，然后在File中选择SavePath后退出即可。
这样，你可以在MATLAB的命令窗口中键入你的源程序文件名，或在编辑窗口中选择Tools中的Run即可编译运行。
若编译无错，则可得出结果；若编译有错，则可根据命令窗口的提示进行修改后再编译运行，直至得出正确的结果。
附录B 信号处理工具箱函数
MATLAB包含了进行信号处理的许多工具箱函数，有关这些工具箱函数的使用可通过Help命令得到。
为使用方便，在这里将给出几个常用到的函数的使用说明。
函数形式
函数功能
关于函数参数的说明
X=sawtooth（t,width）
产生锯齿波或三角波。
width用于确定最大值的位置，即从0到2
*width函数从-1上升到+1。
X=square（t,duty）
产生方波
Duty用于指定正半周期的比例
Y=abs（x）
求绝对值
当x为复数时，得到的是复数模（幅值），若x为字符串，得到的是各个字符的ASCII码。
C=conv（a,b）
求卷积
求取矢量a和b的卷积,c的长度为a和b的长度和减去1。
[h,w]=freqs（b,a,n）
模拟滤波器的频率响应
b,a为滤波器的冲击响应s变换的分子和分母多项式的系数，在n个频率点计算频率响应
[h,f]=freqz（b,a,n,Fs）
数字滤波器的频率响应。
Fs为采样频率，b,a为滤波器的冲击响应的Z变换的分子和分母多项式的系数，该函数的作用是在0~Fs/2频率范围内选取n个点（记在f中），并计算相应的频率响应。
[h,t]=impz（b,a,n）
数字滤波器的冲击响应
b,a为滤波器的冲击响应s变换的分子和分母多项式的系数，计算出冲击响应h，取样点数为n.
[n,Wn]=buttord（
Wp,Ws,Rp,Rs,’s’）
Butterworth滤波器阶的选择。
Wp和Ws分别为通带和阻带的截止频率，皆大于0小于1。
Rp和Rs分别是通带和阻带的波纹系数，’s’表示模拟域，也可不加‘s’，则为数字域。
[b,a]=butter（n,Wn,’ftype’,’s’）
Butterworth模拟和数字滤波器设计。
设计阶数为n，3dB截止频率为Wn的滤波器，ftype指滤波器的类型，‘high’是高通，’stop’是带阻，无此参数则是低通，’s’指模拟域，无则表示数字域，b,a是对应变换的分子分母多项式的系数。
[n,Wn]=cheblord（
Wp,Ws,Rp,Rs,’s’）
chebyshevI滤波器阶的选择。
Wp和Ws分别为通带和阻带的截止频率，皆大于0小于1。
Rp和Rs分别是通带和阻带的波纹系数，’s’表示模拟域，也可不加‘s’，则为数字域。
[b,a]=cheby1（n,RpWn,’ftype’,’s’）
Chebyshev（切比雪夫）I型模拟和数字滤波器设计（通带等波纹）。
设计阶数
升级会员 