频谱分析(完整版)

Matlab信号处理工具箱帮助文档谱估计专题翻译：无名网友&Lyra频谱分析Spectralestimation（谱估计）的目标是基于一个有限的数据集合描述一个信号的功率（在频率上的）分布。功率谱估计在很多场合下都是有用的，包括对宽带噪声湮没下的信号的检测。从数学上看，一个平稳随机过程nx的powerspectrum（功率谱）和correlationsequence（相关序列）通过discrete-timeFouriertransform（离散时间傅立叶变换）构成联系。从normalizedfrequency（归一化角频率）角度看，有下式jmxxxxmSRme注：2xxSX，其中/2/21limNjnnNnNXxeN。其matlab近似为X=fft(x,N)/sqrt(N)，在下文中LXf就是指matlabfft函数的计算结果了使用关系2/sff可以写成物理频率f的函数，其中sf是采样频率2/sjfmfxxxxmSfRme相关序列可以从功率谱用IDFT变换求得：/22//22sssfjfmfjmxxxxxxsfSeSfeRmddff序列nx在整个Nyquist间隔上的平均功率可以表示为/2/202ssfxxxxxxsfSSfRddff上式中的2xxxxSP以及xxxxsSfPff被定义为平稳随机信号nx的powerspectraldensity(PSD)（功率谱密度）一个信号在频带1212,,0上的平均功率可以通过对PSD在频带上积分求出211212,xxxxPPdPd从上式中可以看出xxP是一个信号在一个无穷小频带上的功率浓度，这也是为什么它叫做功率谱密度。PSD的单位是功率（e.g瓦特）每单位频率。在xxP的情况下，这是瓦特/弧度/抽或只是瓦特/弧度。在xxPf的情况下单位是瓦特/赫兹。PSD对频率的积分得到的单位是瓦特，正如平均功率12,P所期望的那样。对实信号，PSD是关于直流信号对称的，所以0的xxP就足够完整的描述PSD了。然而要获得整个Nyquist间隔上的平均功率，有必要引入单边PSD的概念：0020onesidedxxPP信号在频带1212,,0上的平均功率可以用单边PSD求出2121,onesidedPPd频谱估计方法Matlab信号处理工具箱提供了三种方法Nonparametricmethods（非参量类方法）PSD直接从信号本身估计出来。最简单的就是periodogram（周期图法），一种改进的周期图法是Welch'smethod。更现代的一种方法是multitapermethod（多椎体法）。Parametricmethods（参量类方法）这类方法是假设信号是一个由白噪声驱动的线性系统的输出。这类方法的例子是Yule-Walkerautoregressive(AR)method和Burgmethod。这些方法先估计假设的产生信号的线性系统的参数。这些方法想要对可用数据相对较少的情况产生优于传统非参数方法的结果。Subspacemethods（子空间类）又称为high-resolutionmethods（高分辨率法）或者super-resolutionmethods（超分辨率方法）基于对自相关矩阵的特征分析或者特征值分解产生信号的频率分量。代表方法有multiplesignalclassification(MUSIC)method或eigenvector(EV)method。这类方法对线谱（正弦信号的谱）最合适，对检测噪声下的正弦信号很有效，特别是低信噪比的情况。方法描述函数周期图PSD估计spectrum.periodogram,periodogramWelch重叠，加窗的信号段的平均周期图spectrum.welch,pwelch,cpsd,tfestimate,mscohere多椎体多个正交窗（称为锥）的组合做谱估计spectrum.mtm,pmtmYule-WalkerAR时间序列的估计的自相关函数计算自回归（AR）谱估计spectrum.yulear,pyulearBurg通过最小化线性预测误差计算自回归（AR）谱估计spectrum.burg,pburgCovariance（协方差）通过最小化前向预测误差做时间序列的自回归（AR）谱估计spectrum.cov,pcov修正协方差通过最小化前向及后向预测误差做时间序列的自回归（AR）谱估计spectrum.mcov,pmcovMUSIC多重信号分类spectrum.music,pmusic特征向量法虚谱估计spectrum.eigenvector,peigNonparametricMethods非参数法下面讨论periodogram,modifiedperiodogram,Welch,和multitaper法。同时也讨论CPSD函数，传输函数估计和相关函数。Periodogram周期图法一个估计功率谱的简单方法是直接求随机过程抽样的DFT，然后取结果的幅度的平方。这样的方法叫做周期图法。一个长L的信号Lxn的PSD的周期图估计是2ˆLxxsXfPffL注：这里LXf运用的是matlab里面的fft的定义不带归一化系数，所以要除以L其中12/0sLjfnfLLnXfxne实际对LXf的计算可以只在有限的频率点上执行并且使用FFT。实践上大多数周期图法的应用都计算N点PSD估计2ˆLkxxksXfPffL，,0,1,,1skkffkNN其中12/0LjknNLkLnXfxne选择N是大于L的下一个2的幂次是明智的，要计算LkXf我们直接对Lxn补零到长度为N。假如LN，在计算LkXf前，我们必须绕回Lxn模N。作为一个例子，考虑下面1001元素信号nx，它包含了2个正弦信号和噪声randn('state',0);fs=1000;%Samplingfrequencyt=(0:fs)/fs;%OnesecondworthofsamplesA=[12];%Sinusoidamplitudes(rowvector)f=[150;140];%Sinusoidfrequencies(columnvector)xn=A*sin(2*pi*f*t)+0.1*randn(size(t));注意：最后三行表明了一个方便的表示正弦之和的方法，它等价于：xn=sin(2*pi*150*t)+2*sin(2*pi*140*t)+0.1*randn(size(t));对这个PSD的周期图估计可以通过产生一个周期图对象（periodogramobject）来计算Hs=spectrum.periodogram('Hamming');估计的图形可以用psd函数显示。psd(Hs,xn,'Fs',fs,'NFFT',1024,'SpectrumType','twosided')00.10.20.30.40.50.60.70.80.9-80-70-60-50-40-30-20-100Frequency(kHz)Power/frequency(dB/Hz)PowerSpectralDensityEstimateviaPeriodogram平均功率通过用下述求和去近似积分求得[Pxx,F]=psd(Hs,xn,fs,'twosided');Pow=(fs/length(Pxx))*sum(Pxx)Pow=2.5059你还可以用单边PSD去计算平均功率[Pxxo,F]=psd(Hs,xn,fs,'onesided');Pow=(fs/(2*length(Pxxo)))*sum(Pxxo)Pow=2.5011周期图性能下面从四个角度讨论周期图法估计的性能：泄漏，分辨率，偏差和方差。频谱泄漏考虑有限长信号Lxn，把它表示成无限长序列xn乘以一个有限长矩形窗Rwn的乘积的形式经常很有用：LRxnxnwn因为时域的乘积等效于频域的卷积，所以上式的傅立叶变换是/2/21ssfLRsfXfXWfdf前文中导出的表达式2ˆLxxsXfPffL说明卷积对周期图有影响。正弦数据的卷积影响最容易理解。假设xn是M个复正弦的和1kMjnkkxnAe其频谱是1MskkkXffAff对一个有限长序列，就变成了/211/21ssfMMLskkRkRkkksfXffAfWfdAWfff所以在有限长信号的频谱中，Dirac函数被替换成了形式为RkWff的项，该项对应于矩形窗的中心在kf的频率响应。一个矩形窗的频率响应形状是一个sinc信号，如下所示-500-400-300-200-1000100200300400500-80-70-60-50-40-30-20-100矩形窗在物理频率上的功率谱密度frequency/HzPSDdBwatt/Hz该图显示了一个主瓣和若干旁瓣，最大旁瓣大约在主瓣下方13.5dB处。这些旁瓣说明了频谱泄漏效应。无限长信号的功率严格的集中在离散频率点kf处，而有限长信号在离散频率点kf附近有连续的功率。因为矩形窗越短，它的频率响应对Dirac冲击的近似性越差，所以数据越短它的频谱泄漏越明显。考虑下面的100个采样的序列randn('state',0)fs=1000;%Samplingfrequencyt=(0:fs/10)/fs;%One-tenthofasecondworthofsamplesA=[12];%Sinusoidamplitudesf=[150;140];%Sinusoidfrequenciesxn=A*sin(2*pi*f*t)+0.1*randn(size(t));Hs=spectrum.periodogram;psd(Hs,xn,'Fs',fs,'NFFT',1024)00.050.10.150.20.250.30.350.40.450.5-80-70-60-50-40-30-20-100Frequency(kHz)Power/frequency(dB/Hz)PowerSpectralDensityEstimateviaPeriodogram注意到频谱泄露只视数据长度而定。周期图确实只对有限数据样本进行计算，但是这和频谱泄露无关。分辨率分辨率指的是区分频谱特征的能力，是分析谱估计性能的关键概念。要区分两个在频率上离得很近的正弦，要求两个频率差大于任何一个信号泄漏频谱的主瓣宽度。主瓣宽度定义为主瓣上峰值功率一半的点间的距离（3dB带宽）。该宽度近似等于/sfL两个频率为1f2f的正弦信号，可分辨条件是12sffffL上例中频率间隔10Hz，数据长度要大于100抽才能使得周期图中两个频率可分辨。下图是只有67个数据长度的情况randn('state',0)fs=1000;%Samplingfrequencyt=(0:fs/15)./fs;%67samplesA=[12];%Sinusoidamplitudesf=[150;140];%Sinusoidfrequenciesxn=A*sin(2*pi*f*t)+0.1*randn(size(t));Hs=spectrum.periodogram;psd(Hs,xn,'Fs',fs,'NFFT',1024)00.050.10.150.20.250.30.350.40.450.5-70-60-50-40-30-20-100Frequency(kHz)Power/frequency(dB