基于Matlab的频谱分析

基于Matlab的频谱分析一、实验目的1、掌握时域抽样定理。2、通过实验加深对FFT的理解。3、熟悉应用FFT对典型信号进行频谱分析的方法。二、实验原理1、时域抽样定理时域抽样定理给出了连续信号抽样过程中信号不失真的约束条件：对于基带信号，信号抽样频率大于等于2倍的信号最高频率，即。时域抽样是把连续信号变成适于数字系统处理的离散信号。对连续信号以间隔T抽样，则可得到的离散序列为。图1连续信号抽样的离散序列若，则信号与的频谱之间存在：其中：的频谱为，的频谱为。可见，信号时域抽样导致信号频谱的周期化。(rad/s)x(t)x[k]tk0T2T012)e(jXnnXT)(j1sam)e(jXkX)e(jX)j(XTsam/2kXtXkTtkTXXktXkXtXkTtkTXXkmsamff2samfmf为抽样角频率，为抽样频率。数字角频率Ω与模拟角频率ω的关系为：Ω=ωT。2、离散傅立叶变换（DFT）有限长序列)(nx的离散傅立叶变换（DFT）为10,)()]([)(10NnWnxnxDFTkXNnknN逆变换为10,)(1)]([)(10NnWkXNkXIDFTnxNnknN3、快速傅立叶变换（FFT）在各种信号序列中，有限长序列占重要地位。对有限长序列可以利用离散傅立叶变换(DFT)进行分析。DFT不但可以很好的反映序列的频谱特性，而且易于用快速算法(FFT)在计算机上进行分析。有限长序列的DFT是其z变换在单位圆上的等距离采样，或者说是序列傅立叶的等距离采样，因此可以用于序列的谱分析。FFT是DFT的一种快速算法，它是对变换式进行一次次分解，使其成为若干小数据点的组合，从而减少运算量。MATLAB为计算数据的离散快速傅立叶变换，提供了一系列丰富的数学函数，主要有Fft、Ifft、Fft2、Ifft2,Fftn、ifftn和Fftshift、Ifftshift等。当所处理的数据的长度为2的幂次时，采用基-2算法进行计算，计算速度会显著增加。所以，要尽可能使所要处理的数据长度为2的幂次或者用添零的方式来添补数据使之成为2的幂次。Fft函数调用方式：○1、Y=fft(X)；○2、Y＝fft(X,N)；○3、Y＝fft(X,[],dim)或Y＝fft(X,N,dim)。Tfsam1函数Ifft的参数应用与函数Fft完全相同。三.实验内容与步骤1、为了说明高密度频谱和高分辨率频谱之间的区别，考察序列x(n)有限个样本的频谱。)52.0cos()48.0cos()(nnnx当100n时，求x(n)的DFT.当1000n时，求x(n)的DFT.程序:functionDFTX1N=input('请输入离散信号的长度：');n=0:N-1;k=0:N-1;xn=cos(0.48*n*pi)+cos(0.52*n*pi);WN=exp(-j*2*pi/N);nk=n'*k;WNnk=WN.^nk;Xk=xn*WNnk;subplot(2,1,1)stem(n,xn);title('x(n):');gridon;subplot(2,1,2)stem(k,abs(Xk));title('x(n)的DFT:');gridon;DFTX1请输入离散信号的长度：10DFTX1请输入离散信号的长度：100频谱图：0123456789-2-1012当0≤n≤10时，序列x(n)图：012345678902468当0≤n≤10时，序列x(n)的DFT图:0102030405060708090100-2-1012当0≤n≤100时，序列x(n)图:01020304050607080901000204060当0≤n≤100时，序列x(n)的DFT图:2、已知一模拟信号()()taxteut，现在采样率20sfHz进行采样。用DFT计算当序列长度分别为L=100,L=20时，N=200点的幅度频谱样值并通过作图与理论上准确的频谱样值进行比较。程序：L1=100;L2=20;t1=0:L1-1;x1=exp(-t1/20);t2=0:L2-1;x2=exp(-t2/20);N=0:199;WN=exp(-j*2*pi/length(N));nk1=t1'*N;nk2=t2'*N;WNnk1=WN.^nk1;WNnk2=WN.^nk2;Xk1=x1*WNnk1;Xk2=x2*WNnk2;subplot(4,1,1)stem(t1,x1);title('序列长度L=100时x(n)图:');gridon;subplot(4,1,3)stem(t2,x2);title('序列长度L=20时x(n)图:');gridon;subplot(4,1,2)stem(N,abs(Xk1));title('序列长度L=100时x(n)的DFT图:');gridon;subplot(4,1,4)stem(N,abs(Xk2));title('序列长度L=20时x(n)的DFT图:');gridon;010203040506070809010000.51序列长度L=100时x(n)图:010203040506070809010000.51序列长度L=20时x(n)图:0204060801001201401601802000102030序列长度L=100时x(n)的DFT图:020406080100120140160180200051015序列长度L=20时x(n)的DFT图:实验图：3、一个连续信号含三个频谱分量，经采样得以下序列：))215.0(2cos())15.0(2cos()15.02sin()(ndfndfnnx（1）N＝64，df分别为161、1/64，观察其频谱；程序：functiondftx3(df)N=64;k=0:127;n=0:N-1;xn=sin(2*pi*0.15*n)+cos(2*pi*(0.15+df)*n)+cos(2*pi*(0.15+2*df)*n);WN=exp(-j*2*pi/length(k));nk=n'*k;WNnk=WN.^nk;Xk=xn*WNnk;subplot(2,1,1)plot(n,xn);title('df=1/16,长度为N=64的有限序列x(n)图:');gridon;subplot(2,1,2)plot(abs(Xk(1:64)));title('取上图x(n)做128点的DFT图:');gridon;dftx3(1/16)dftx3(1/64)实验图：010203040506070-4-2024df=1/16,长度为N=64的有限序列x(n)图:010203040506070010203040取上图x(n)做128点的DFT图:010203040506070-4-2024df=1/64,长度为N=64的有限序列x(n)图:010203040506070010203040取上图x(n)做128点的DFT图:020406080100120140-505df=1/64时，有限序列x(n):01020304050607002040取上图有限序列x(n)长度N=64的128点FFT图:010203040506070050100取上上图有限序列x(n)长度N=128的128点FFT图:（2）N＝64、128，df为1/64，做128点得FFT，其结果有何不同？程序：n1=0:63;n2=0:127;df=1/64;x1=sin(2*pi*0.15*n1)+cos(2*pi*(0.15+df)*n1)+cos(2*pi*(0.15+2*df)*n1);x2=sin(2*pi*0.15*n2)+cos(2*pi*(0.15+df)*n2)+cos(2*pi*(0.15+2*df)*n2);y1=fft(x1,128);y2=fft(x2,128);subplot(3,1,1)plot(n2,x2);title('有限序列x(n):');gridon;subplot(3,1,2)plot(abs(y1(1:64)));title('取有限序列x(n)长度N=64的128点FFT图:');gridon;subplot(3,1,3)plot(abs(y2(1:64)));title('取有限序列x(n)长度N=128的128点FFT图:');gridon;实验图：4、被噪声污染的信号，较难看出所包含得频率分量，如一个由50Hz和120Hz正弦信号构成的信号，受零均值随机噪声的干扰，数据采样率为1000Hz，试用FFT函数来分析其信号频率成分，要求：（1）画出时域波形；（2）分析信号功率谱密度。注：在MATLAB中，可用函数rand（1，N）产生均值为0，方差为1，长度为N的高斯随机序列。t=0:0.001:1;x=sin(2*pi*50*t)+sin(2*pi*120*t)+rand(size(t));subplot(2,1,1)plot(x(1:50));title('由50Hz和120Hz正弦信号构成且受噪声污染的信号图:');gridon;Y=fft(x,512);subplot(2,1,2)plot(abs(Y).^2);title('受污染信号功率谱密度图:');gridon;05101520253035404550-2-10123由50Hz和120Hz正弦信号构成且受噪声污染的信号图:010020030040050060002468x104受污染信号功率谱密度图: