FFT+STFT+经典功率谱+现代功率谱 信号分析实验报告+matlab程序.doc

实验1.1：用MATLAB进行FFT频谱分析假设一信号：画出其频谱图。分析：首先，连续周期信号截断对频谱的影响。DFT变换频谱泄漏的根本原因是信号的截断。即时域加窗，对应为频域卷积，因此，窗函数的主瓣宽度等就会影响到频谱。实验表明，连续周期信号截断时持续时间与信号周期呈整数倍关系时，利用DFT变换可以得到精确的模拟信号频谱。举一个简单的例子：其周期为0.02。截断时不同的持续时间影响如图一.1：（对应程序shiyan1ex1.m）图 Error! No text of specified style in document.1其次，采样频率的确定。根据Shannon采样定理，采样带限信号采样频率为截止频率的两倍以上，给定信号的采样频率应1/7.92，取16。再次，DFT算法包括时域采样和频域采样两步，频域采样长度M和时域采样长度N的关系要符合MN时，从频谱X(k)才可完全重建原信号。实验中信号R经采样后的离散信号不是周期信号，但是它又是一个无限长的信号，因此处理时时域窗函数尽量取得宽一些已接近实际信号。实验结果如图一.2：其中，0点位置的冲激项为直流分量0.6造成（对应程序为shiyan1.m）图 Error! No text of specified style in document.2ARMA（Auto Recursive Moving Average）模型：将平稳随机信号x(n)看作是零均值，方差为u2的白噪声u(n)经过线性非移变系统H(z)后的输出，模型的传递函数为用差分方程表示为AR（Auto Recursive）自回归模型，即ARMA模型中系数b只有在r=0的情况下为1，其余都是零，获得一个全极点模型：差分方程表示为：AR模型的功率谱估计为：（摘自数字信号处理技术的算法分析与应用祁才君编著）实验3.1：用MATLAB进行AR模型功率谱分析随机信号序列x(n)是均值为0方差为1的高斯型白噪声经过AR模型后的输出，采样长度为512，AR模型阶次取3，4，5，用L-D算法估计功率谱密度。分析：MATLAB函数pyulear()的用法pyulear()是基于自相关法、利用Levesion-Durbin算法估计功率谱密度。px,w=pyulear(x,p,nfft,range)x为随即信号序列，是由白噪声经AR模型产生的，在MATLAB中可以由白噪声序列u经过表示AR模型的数字滤波器后得到，使用的是filter函数；p为AR模型阶次；nfft为由模型参数计算频谱时的频域采样点数，默认为256；range用于选择输出是为单边0,，还是双边0,2；w的范围0,，还是 0,2由range确定或由nfft的奇偶性确定；该函数返回实际频率w下的功率谱密度向量，w的单位即为rad/sample，默认sample为1Hz，若要转化为归一化频率，只需用w/即可。实验结果如图三.1（对应程序为shiyan3.m）：图 Error! No text of specified style in document.3短时傅里叶变换（Short Time Fourier Transform, STFT）法，在MATLAB中做短时傅里叶变换的函数为spectrogram：spectrogram(x,window,overlap,f,fs)s,f,t,p=spectrogram(x,window,overlap,f,fs)x为被分析序列，window为窗函数及长度，默认为hamming窗，overlap为相邻两个短时序列之间重叠的数据点数，f为一向量，确定在某一个频率范围内做短时傅里叶变换，fs为采样频率。实验5.1：用MATLAB对给定信号做短时傅里叶分析。 采样频率取1024Hz，分析窗之间的滑动因子为20，分析窗用kaiser窗，长度分别取60和260时的STFT变换。绘制的频谱图如下图五.1：（对应程序为shiyan5.m）图 Error! No text of specified style in document.4功率谱估计（Power Spectrum Density, PSD）是用有限长的数据来估计信号的功率谱，方法有经典功率谱估计法和现代功率谱估计法，前者常用的为周期图法，后者常用的是AR模型法。实验6.1：周期图法估计平稳随机信号功率谱的MATLAB实现分析：使用周期图法估计功率谱的MATLAB函数为psd()（已被pwelch()代替了，实验中使用后者）：welch方法为改进的功率谱法，是用相干平均对周期图法功率谱估计进行二次处理，改善功率谱的估计方差，welch方法对取样序列采取重叠分段方法，MATLAB实现函数为pwelch()：pxx,w=pwelch(x,window,noverlap,nfft)pxx,f=pwelch(x,window,noverlap,nfft,fs)各参数的含义和其余谱估计方法相似，window默认为hamming窗，使用前者为产生在w下的谱密度，w的单位为rad/sample，实信号w范围为0,复信号为0,2；后者fs为采样频率，获得单位Hz下的谱密度，f的范围为0,f/2或0,fs，判断与w相同。无输出参数时，MATLAB画出相应的频谱图。假设其中，v(n)为均值0方差1的随机噪声信号，用welch法估计其功率谱密度。在使用pwelch时一定要注意w或f与pxx的对应关系，当采样频率fs改变时，w及f要做出相应的调整。结果如下图六.1（对应程序为shiyan6.m）。实验中可以看出，随分段长度的增加，即时域窗函数长度的增加，功率谱估计的方差特性变差，为改变周期图法功率谱估计方差特性差的缺点，Thomson提出了功率谱的多窗估计法（Multitaper Method），图六.2为用多窗估计法对上例中的x(n)做功率谱估计（对应程序仍为shiyan6.m）。MATLAB中使用多窗估计法估计功率谱的函数为pmtm()：pxx,pxxc,w=pmtm(x,nw,nfft,fs)pxx为功率谱估计，nw为离散扁长球体序列（Discrete Prolate Spheroidal Sequences, DPSS, or Slepian Sequences），默认值为4，其余参数的含义与以上相同。实验中