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

-1-实验1.1：用MATLAB进行FFT频谱分析假设一信号：292.7/2cos1.0996.2/2sin1.06.0ttR画出其频谱图。分析：首先，连续周期信号截断对频谱的影响。DFT变换频谱泄漏的根本原因是信号的截断。即时域加窗，对应为频域卷积，因此，窗函数的主瓣宽度等就会影响到频谱。实验表明，连续周期信号截断时持续时间与信号周期呈整数倍关系时，利用DFT变换可以得到精确的模拟信号频谱。举一个简单的例子：2.0100costY其周期为0.02。截断时不同的持续时间影响如图一.1：（对应程序shiyan1ex1.m）00.0020.0040.0060.0080.010.0120.0140.0160.0180.02-1-0.500.51函数y=cos(100t+0.2)01002003004005006007008009001000020406080100截断时，时间间期为周期整数倍，频谱图00.0050.010.0150.020.0250.03-1-0.500.51函数y=cos(100t+0.2)01002003004005006007008009001000020406080100截断时，时间间期不为周期整数倍，频谱图图错误!文档中没有指定样式的文字。.1-2-其次，采样频率的确定。根据Shannon采样定理，采样带限信号采样频率为截止频率的两倍以上，给定信号的采样频率应1/7.92，取16。再次，DFT算法包括时域采样和频域采样两步，频域采样长度M和时域采样长度N的关系要符合M≧N时，从频谱X(k)才可完全重建原信号。实验中信号R经采样后的离散信号不是周期信号，但是它又是一个无限长的信号，因此处理时时域窗函数尽量取得宽一些已接近实际信号。实验结果如图一.2：其中，0点位置的冲激项为直流分量0.6造成（对应程序为shiyan1.m）0204060801001201401601802000.40.50.60.70.800.050.10.150.20.250.30.350.40.450.5050100150图错误!文档中没有指定样式的文字。.2♣ARMA（AutoRecursiveMovingAverage）模型：将平稳随机信号x(n)看作是零均值，方差为σu2的白噪声u(n)经过线性非移变系统H(z)后的输出，模型的传递函数为-3-PkkkQrrrzazbzAzBzH111用差分方程表示为QrrPkkrnubknxanx01AR（AutoRecursive）自回归模型，即ARMA模型中系数b只有在r=0的情况下为1，其余都是零，获得一个全极点模型：PkkkzazAzBzH111差分方程表示为：nuknxanxPkk1AR模型的功率谱估计为：jjujxeAeAeS12（摘自《数字信号处理技术的算法分析与应用》祁才君编著）实验3.1：用MATLAB进行AR模型功率谱分析随机信号序列x(n)是均值为0方差为1的高斯型白噪声经过AR模型43219606.01697.29403.22137.211zzzzzH后的输出，采样长度为512，AR模型阶次取3，4，5，用L-D算法估计功率谱密度。分析：MATLAB函数pyulear()的用法-4-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）：00.20.40.60.81050100150200PSDestimatedforAR(3)NormalizedFrequency/RelativePowerMagnititude00.20.40.60.810100200300400PSDestimatedforAR(4)NormalizedFrequency/RelativePowerMagnititude00.10.20.30.40.50.60.70.80.9102004006008001000PSDestimatedforAR(5)NormalizedFrequency/RelativePowerMagnititude图错误!文档中没有指定样式的文字。.3短时傅里叶变换（ShortTimeFourierTransform,STFT）法，在MATLAB中做短时傅里叶变换的函数为spectrogram：-5-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对给定信号做短时傅里叶分析。tttx1502cos2.0502cosstsst2110采样频率取1024Hz，分析窗之间的滑动因子为20，分析窗用kaiser窗，长度分别取60和260时的STFT变换。绘制的频谱图如下图五.1：（对应程序为shiyan5.m）图错误!文档中没有指定样式的文字。.4-6-功率谱估计（PowerSpectrumDensity,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画出相应的频谱图。假设nvnnnx4.0cos35.0cos其中，v(n)为均值0方差1的随机噪声信号，用welch法估计其功率谱密度。在使用pwelch时一定要注意w或f与pxx的对应关系，当采样频率fs改变时，w及f要做出相应的调整。结果如下图六.1（对应程序为shiyan6.m）。实验中可以看出，随分段长度的增加，即时域窗函数长度的增加，功率谱估计的方差特性变差，为改变周期图法功率谱估计方差特性差的缺点，Thomson提出了功率谱的多窗估计法（MultitaperMethod），图六.2为用多窗估计法对上例中的x(n)做功率谱估计（对应程序仍为shiyan6.m）。MATLAB中使用多窗估计法估计功率谱的函数为pmtm()：[pxx,pxxc,w]=pmtm(x,nw,nfft,fs)pxx为功率谱估计，nw为离散扁长球体序列（DiscreteProlateSpheroidal-7-Sequences,DPSS,orSlepianSequences），默认值为4，其余参数的含义与以上相同。实验中可看出，随窗函数的增加，方差特性得到进一步改善，但是频率分辨率下降。00.10.20.30.40.50.60.70.80.91-20-1001020NormalizedFrequency(rad/sample)Power/Frequency(dB/rad/sample)00.10.20.30.40.50.60.70.80.91-20-1001020NormalizedFrequency(rad/sample)Power/Frequency(dB/rad/sample)PowerSpectralDensityEstimateviaWelch(modifiedperiodogrammethod)图错误!文档中没有指定样式的文字。.500.10.20.30.40.50.60.70.80.91-20-1001020NormalizedFrequency(rad/sample)Power/Frequency(dB/rad/sample)PowerSpectralDensityEstimateviaThomsonMultitaper00.10.20.30.40.50.60.70.80.91-20-1001020NormalizedFrequency(rad/sample)Power/Frequency(dB/rad/sample)图错误!文档中没有指定样式的文字。.6