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

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

-1-实验1.1：用MATLAB进行FFT频谱分析假设一信号：29.7/2cos1.096.2/sin1.06.ttR画出其频谱图。分析：首先，连续周期信号截断对频谱的影响。DFT变换频谱泄漏的根本原因是信号的截断。即时域加窗，对应为频域卷积，因此，窗函数的主瓣宽度等就会影响到频谱。实验表明，连续周期信号截断时持续时间与信号周期呈整数倍关系时，利用DFT变换可以得到精确的模拟信号频谱。举一个简单的例子：2.01cos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不不不不不不不不不不不不不不不不不不不图Error!Notextofspecifiedstyleindocument1-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图Error!Notextofspecifiedstyleindocument2♣ARMA(AutoRecursiveMovingAverage)模型：将平稳随机信号x(n)看作是零均值，方差为σu2的白噪声u(n)经过线性非移变系统H(z)后的输出，模型的传递函数为-3-PkkQrrzabzABH11用差分方程表示为QrrPkrnubknxanx01AR(AutoRecursive)自回归模型，即ARMA模型中系数b只有在r=0的情况下为1，其余都是零，获得一个全极点模型：PkkzazABH1差分方程表示为：nuknxanxPk1AR模型的功率谱估计为：jjujxeAeS12(摘自《数字信号处理技术的算法分析与应用》祁才君编著)实验3.1：用MATLAB进行AR模型功率谱分析随机信号序列x(n)是均值为0方差为1的高斯型白噪声经过AR模型4321960.697.943.237.1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图Error!Notextofspecifiedstyleindocument3短时傅里叶变换(ShortTimeFourierTrans,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对给定信号做短时傅里叶分析。ttx1502cos.0sts2110采样频率取1024Hz，分析窗之间的滑动因子为20，分析窗用kaiser窗，长度分别取60和260时的STFT变换。绘制的频谱图如下图五.1：(对应程序为shiyan5.m)图Error!Notextofspecifiedstyleindocument4-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，实