基于Matlab的相关频谱分析程序教程

Matlab 信号处理工具箱 谱估计专题 频谱分析 Spectral estimation(谱估计)的目标是基于一个有限的数据集合描述一个信号的功率(在频率上的)分布。功率谱估计在很多场合下都是有用的，包括对宽带噪声湮没下的信号的检测。 从数学上看，一个平稳随机过程的power spectrum(功率谱)和correlation sequence(相关序列)通过discrete-time Fourier trans(离散时间傅立叶变换)构成联系。从normalized frequency(归一化角频率)角度看，有下式 注：，其中。其matlab近似为X=fft(x,N)/sqrt(N)，在下文中就是指matlab fft函数的计算结果了 使用关系可以写成物理频率的函数，其中是采样频率 相关序列可以从功率谱用IDFT变换求得： 序列在整个Nyquist间隔上的平均功率可以表示为 上式中的 以及 被定义为平稳随机信号的power spectral density (PSD)(功率谱密度) 一个信号在频带上的平均功率可以通过对PSD在频带上积分求出 从上式中可以看出是一个信号在一个无穷小频带上的功率浓度，这也是为什么它叫做功率谱密度。 PSD的单位是功率(e.g 瓦特)每单位频率。在的情况下，这是瓦特/弧度/抽或只是瓦特/弧度。在的情况下单位是瓦特/赫兹。PSD对频率的积分得到的单位是瓦特，正如平均功率所期望的那样。 对实信号，PSD是关于直流信号对称的，所以的就足够完整的描述PSD了。然而要获得整个Nyquist间隔上的平均功率，有必要引入单边PSD的概念： 信号在频带上的平均功率可以用单边PSD求出 频谱估计方法 Matlab 信号处理工具箱提供了三种方法 Nonparametric s(非参量类方法) PSD直接从信号本身估计出来。最简单的就是periodogram(周期图法)，一种改进的周期图法是Welch s 。更现代的一种方法是multitaper (多椎体法)。 Parametric s (参量类方法) 这类方法是假设信号是一个由白噪声驱动的线性系统的输出。这类方法的例子是Yule-Walker autoregressive (AR) 和Burg 。这些方法先估计假设的产生信号的线性系统的参数。这些方法想要对可用数据相对较少的情况产生优于传统非参数方法的结果。 Subspace s (子空间类) 又称为high-resolution s(高分辨率法)或者super-resolution s(超分辨率方法)基于对自相关矩阵的特征分析或者特征值分解产生信号的频率分量。代表方法有multiple signal classification (MUSIC) 或eigenvector (EV) 。这类方法对线谱(正弦信号的谱)最合适，对检测噪声下的正弦信号很有效，特别是低信噪比的情况。 方法 描述 函数 周期图 PSD 估计 spectrum.periodogram, periodogram Welch 重叠，加窗的信号段的平均周期图 spectrum.welch, pwelch, cpsd, tfestimate, mscohere 多椎体 多个正交窗(称为锥)的组合做谱估计 spectrum.mtm, pmtm Yule-Walker AR 时间序列的估计的自相关函数计算自回归(AR)谱估计 spectrum.yulear, pyulear Burg 通过最小化线性预测误差计算自回归(AR)谱估计 spectrum.burg, pburg Covariance(协方差) 通过最小化前向预测误差做时间序列的自回归(AR)谱估计 spectrum.cov, pcov 修正协方差 通过最小化前向及后向预测误差做时间序列的自回归(AR)谱估计 spectrum.mcov, pmcov MUSIC 多重信号分类 spectrum.music, pmusic 特征向量法 虚谱估计 spectrum.eigenvector, peig Nonparametric s非参数法 下面讨论periodogram, modified periodogram, Welch, 和 multitaper法。同时也讨论CPSD函数，传输函数估计和相关函数。 Periodogram周期图法 一个估计功率谱的简单方法是直接求随机过程抽样的DFT，然后取结果的幅度的平方。这样的方法叫做周期图法。 一个长L的信号的PSD的周期图估计是 注：这里运用的是matlab里面的fft的定义不带归一化系数Matlab FFT函数未做归一化 ，所以要除以L 其中 实际对的计算可以只在有限的频率点上执行并且使用FFT。实践上大多数周期图法的应用都计算N点PSD估计 ， 其中 选择N是大于L的下一个2的幂次是明智的，要计算我们直接对补零到长度为N。假如L>N，在计算前，我们必须绕回模N。 作为一个例子，考虑下面1001元素信号，它包含了2个正弦信号和噪声 randn( state ,0); fs = 1000; % Sampling frequency t = (0:fs)/fs; % One second worth of samples A = [1 2]; % Sinusoid amplitudes (row vector) f = [150;140]; % Sinusoid frequencies (column vector) 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的周期图估计可以通过产生一个周期图对象(periodogram object)来计算 Hs = spectrum.periodogram( Hamming ); 估计的图形可以用psd函数显示。 psd(Hs,xn, Fs ,fs, NFFT ,1024, SpectrumType , twosided ) 平均功率通过用下述求和去近似积分 求得 [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(