1
Matlab 信号处理工具箱 谱估计专题
频谱分析
Spectral estimation（谱估计）的目标是基于一个有限的数据集合描述一个信号的功率（在频率上的）分布。功率谱估计在很多场合下都是有用的，包括对宽带噪声湮没下的信号的检测。
从数学上看，一个平稳随机过程xn的power spectrum（功率谱）和correlation sequence（相关序列）通过discrete-time Fourier transform（离散时间傅立叶变换）构成联系。从normalized frequency（归一化角频率）角度看，有下式
?2
3 4 5
6 7 8 9
10
Sxx????m????R?m?exx?j?m
11
注：Sxx????X???21，其中X????limN??Nn??N/2?N/2xnej?n??????。其matlab
12 13
近似为X=fft(x,N)/sqrt(N)，在下文中XL?f?就是指matlab fft函数的计算结果了
使用关系??2?f/fs可以写成物理频率f的函数，其中fs是采样频率
?14
15
Sxx?f??m????R?m?exx?2?jfm/fs
16
相关序列可以从功率谱用IDFT变换求得：
?sSxx???ej?mSxx?f?e2?jfm/fsd???df
2?fs?fs/2f/217
Rxx?m?????1
18
序列xn在整个Nyquist间隔上的平均功率可以表示为
?sSxx???S?f?d???xxdf 2?fs?fs/2f/219
Rxx?0?????20
上式中的
Sxx???S?f?以及Pxx?f??xx 2?fs21
Pxx????22
被定义为平稳随机信号xn的power spectral density (PSD)（功率谱密度） 一个信号在频带??1,?2?,0??1??2??上的平均功率可以通过对PSD在频带上积分求出
?2??123 24
25
P??1,?2??P???d???P???d? ???xxxx1?226 27
从上式中可以看出Pxx???是一个信号在一个无穷小频带上的功率浓度，这也是为什么它叫做功率谱密度。
PSD的单位是功率（e.g 瓦特）每单位频率。在Pxx???的情况下，这是瓦特/弧度/抽或只是瓦特/弧度。在Pxx?f?的情况下单位是瓦特/赫兹。PSD对频率的积分得到的单位是瓦特，正如平均功率P??1,?2?所期望的那样。
对实信号，PSD是关于直流信号对称的，所以0????的Pxx???就足够完整的描述PSD了。然而要获得整个Nyquist间隔上的平均功率，有必要引入单边PSD的概念：
28 29 30
31 32 33
34
?????0?0Ponesided?????
?2Pxx???0????2
35
信号在频带??1,?2?,0??1??2??上的平均功率可以用单边PSD求出
?236
P??1,?2??P??1onesided???d?
37
频谱估计方法
Matlab 信号处理工具箱提供了三种方法
38 39 40 41 42
PSD直接从信号本身估计出来。最简单的就是periodogram（周期图法），一种改进的周期图法是Welch's method。更现代的一种方法是multitaper method（多椎体法）。
Parametric methods （参量类方法）
这类方法是假设信号是一个由白噪声驱动的线性系统的输出。这类方法的例子是Yule-Walker autoregressive (AR) method和Burg method。这些方法先估计假设的产生信号的线性系统的参数。这些方法想要对可用数据相对较少的情况产生优于传统非参数方法的结果。 Subspace methods （子空间类）
又称为high-resolution methods（高分辨率法）或者super-resolution methods（超分辨率方法）基于对自相关矩阵的特征分析或者特征值分解产生信号的频率分量。代表方法有multiple signal classification (MUSIC) method或eigenvector (EV) method。这类方法对线谱（正弦信号的谱）最合适，对检测噪声下的正弦信号很有效，特别是低信噪比的情况。 Nonparametric Methods非参数法
下面讨论periodogram, modified periodogram, Welch, 和 multitaper法。
3
43
44 45 46 47
48
49 50 51 52 53
54
55
56
同时也讨论CPSD函数，传输函数估计和相关函数。 Periodogram周期图法
一个估计功率谱的简单方法是直接求随机过程抽样的DFT，然后取结果的幅度的平方。这样的方法叫做周期图法。
一个长L的信号xL?n?的PSD的周期图估计是
XL?f?fsL257
58 59
60
61
??f??Pxx 62 63
注：这里XL?f?运用的是matlab里面的fft的定义不带归一化系数，所以要除以L 其中
L?164
65
XL?f???xL?n?e?2?jfn/fs
n?066 67
实际对XL?f?的计算可以只在有限的频率点上执行并且使用FFT。实践上大多数周期图法的应用都计算N点PSD估计
XL?fk?fsL268
??f??Pxxk，fk?kfs,k?0,1,N,N?1 69
其中
L?170
XL?fk???xL?n?e?2?jkn/N
n?071 72
选择N是大于L的下一个2的幂次是明智的，要计算XL?fk?我们直接对xL?n?补零到长度为N。假如L>N，在计算XL?fk?前，我们必须绕回xL?n?模N。
4
73
作为一个例子，考虑下面1001元素信号xn，它包含了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') 74
5
最新Matlab频谱分析程序
