实验四离散信号频谱分析

59实验四离散信号的频域分析一、实验目的1.掌握序列的傅里叶变换、离散傅里叶级数、离散傅里叶变换、快速傅里叶变换的Matlab实现；2.学习用FFT对连续信号和离散信号进行谱分析的方法，了解可能出现的分析误差及其原因，以便在实际中正确应用FFT。二、实验原理及方法1.离散信号的谱分析（1）序列的傅里叶变换对于满足绝对可和的序列，即|)(|nx，其傅里叶变换和反变换的定义为nnjjenxeX)()(（4.1）deeXnxnjj)(21)(（4.2）序列)(nx是离散的，但)(jeX是以2为周期的的连续函数，为了能够在计算机上处理，需要对)(nx进行截断，对频域进行离散化，近似处理后21()()kknjjnnnXexne（4.3）其中2kkM，M是对在一个周期内的采样，k的取值由读者确定，若想观察一个周期内的频谱，0~1kM，若观察两个周期，0~21kM，以此类推。序列傅里叶变换的Matlab实现：clc;clear;n=0:3;x=[1,1,1,1];M=200;k=0:M-1;W=2*pi/M*k;X=x*(exp(-j*2*pi/M)).^(n'*k);%序列的傅里叶变换magX=abs(X);phaX=angle(X);subplot(131);stem(x);subplot(132);plot(W,magX);subplot(133);plot(W,phaX);60对4()Rn进行序列的傅里叶变换得到图4-1。02400.10.20.30.40.50.60.70.80.91时间序列051000.511.522.533.54DTFT幅度谱0510-4-3-2-10123DTFT相位谱图4-1信号及信号的幅度谱和相位谱（2）离散傅里叶变换（DFT）如果序列)(nx是有限长的，序列的谱分析可以采用离散傅里叶变换，其定义为：10,)()]([)(10NkWnxnxDFTkXNnknN（4.4）10,)(1)]([)(10NnWkXNkXIDFTnxNkknN（4.5）因为)(nx与)(kX都是离散的，所以可以利用计算机进行数值计算。从数学观点看，DFT表示的是对序列)(nx或)(kX的线性运算。此处应用DFT变换近似分析采样序列的频谱。设时域序列用x表示，长度为M；x的DFT变换为X，变换区间长度为N（MN）。10)()(MnknNWnxkX1,,1,0Nk将)(kX展开，得：61)1)(1(2)1(1)1(0)1()1(2222120)1(1121110)1(0020100)1()2()1()0()1()1()2()1()0()2()1()2()1()0()1()1()2()1()0()0(MNNNNNNNNMNNNNMNNNNMNNNNWMxWxWxWxNXWMxWxWxWxXWMxWxWxWxXWMxWxWxWxX将上式表示成矩阵的形式：)1)(1()1(1)1(02)1(2021)1(11010)1(1000)]1()1()0([)]1()1()0([MNNMNMNNNNNNNNNNNNN　　　　　　　　　　　　　　　　　　　　　　　　　１DFT变换的Matlab实现：clc;clear;M=4;N=8;x=[1,1,1,1];n=0:M-1;k=0:N-1;WN=exp(-j*2*pi/N);kn=n'*k;WNkn=WN.^kn;X=x*WNkn;magX=abs(X);phaX=angle(X);k=0:7;subplot(131);stem(x);subplot(132);stem(magX);subplot(133);stem(phaX);对4()Rn进行离散傅里叶变换得到图4-2。6202400.10.20.30.40.50.60.70.80.91时间序列051000.511.522.533.54幅度谱0510-2.5-2-1.5-1-0.500.511.5相位谱图4-2信号及信号的离散傅里叶变换（3）快速傅里叶变换（FFT）快速傅里叶变换并不是一种新的变换，只是离散傅里叶变换的快速算法，常用的有按时间抽取的基-2FFT算法和按频率抽取的基-2FFT算法。在Matlab中对离散信号进行FFT，可以直接调用函数。快速傅里叶变换的原理及子程序见附录。fft(x)：利用快速算法计算x的M点DFT，其中M是x的长度。fft(x,N)：利用快速算法计算x的N点DFT，其中N是用户指定的长度。分两种情况：①若x的长度MN，则将x截短为N点序列，再作N点DFT；②若x的长度MN，则将x补零至N点，再作N点DFT。ifft(X)：利用快速算法计算X的M点IDFT，其中M是X的长度。ifft(X,N)：利用快速算法计算X的N点IDFT，其中N是用户指定的长度。同样分两种情况，同fft(x,N)。对4()Rn进行16点的快速傅里叶变换得到图4-3。clc;clear;x=[1,1,1,1];X=fft(x,16);magX=abs(X);phaX=angle(X);k=0:15;W=pi/8*k;subplot(131);stem(x);subplot(132);stem(magX);63subplot(133);stem(phaX);02400.10.20.30.40.50.60.70.80.91时间序列0102000.511.522.533.54幅度谱01020-2-1.5-1-0.500.511.52相位谱图4-3信号及信号的快速傅里叶变换观察图4-1、4-2、4-3可以得到，DFT和FFT都是对序列的傅里叶变换)(jeX在[0,2]上的离散化，离散化的采样点数即为做DFT和FFT的点数。而FFT仅为DFT的快速运算，当做DFT和FFT的点数相同时，两者的结果是一样的。2.利用FFT进行谱分析的误差分析下面分析利用FFT对信号进行傅里叶分析时可能造成的误差。（1）频谱混叠失真利用FFT对连续信号进行傅里叶分析时首先需要进行时域采样。按照时域采样定理，为了不产生混叠，必须满足2scff（4.16）也就是采样间隔T满足112scTff（4.17）一般应取(3~5)scff（4.18）如果不满足2shff，就会产生频谱混叠失真。64（2）栅栏效应利用FFT计算频谱，只能给出离散点2kkN或2kkNT上的频谱采样值，而不可能得到连续频谱函数，这就像通过一个“栅栏”观看信号频谱，只能在离散点上看到信号频谱，称之为“栅栏效应”。这时，如果在两个离散的谱线之间有一个特别大的频谱分量，就无法检测出来了。减小栅栏效应的一个方法就是要使频域采样更密，即增加频域采样点数N，在不改变时域数据的情况下，必然是在数据末端添加一些零值点，使一个周期内的点数增加，但并不改变原有的记录数据。频谱采样间距为2N，N增加，必然使样点间距更近（单位圆上样点更多），谱线更密，谱线变密后原来看不到的谱分量就有可能看到了。必须指出，补零以改变计算FFT的周期时，不能提高频率分辨率，这是因为数据的实际长度仍为补零前的数据长度。（3）频谱泄漏与谱间干扰对信号进行FFT计算，首先必须使其变成有限时宽的信号，这就相当于信号在时域乘一个窗函数如矩形窗，窗内数据并不改变。加窗即为时域相乘，对频域的影响可用卷积公式表示()1()()()21()()2jjjjjVeXeWeXeWed（4.19）式（4.19）中，()jXe是信号()xn的频谱，()jWe是窗函数的频谱。卷积的结果使加窗后得到的频谱()jVe与原来的频谱()jXe不相同，产生失真。这种失真造成频谱的“扩散”，这就是所谓的“频谱泄漏”。三、实验内容及步骤1.计算序列的DTFT和DFT，观察栅栏效应设)()(8nRnx，要求用MATLAB实现：（1）计算)(nx的傅里叶变换)(jeX，并绘出其幅度谱；（2）分别计算)(nx的8点DFT和16点DFT，绘出其幅度谱。并说明它们和)(jeX的关系。（提示：DFT变换可用MATLAB提供的函数fft实现，也可以自己用C语言或matlab编写）2.计算序列的FFT，观察频谱泄漏65已知周期为16的信号)1612cos()1610cos()(nnnx,（1）截取一个周期长度M=16点，计算其16点FFT，并绘出其幅度谱；（2）截取序列长度M=10点，计算其16点FFT，绘出其幅度谱，并与（1）的结果进行比较，观察频谱泄漏现象，说明产生频谱泄漏的原因。四、实验报告要求1.结合实验中所得给定典型序列幅频特性曲线，与理论结果比较，并分析说明误差产生的原因以及用FFT作谱分析时有关参数的选择方法。2.总结实验所得主要结论。