文档介绍：Matlab信号处理工具箱谱估计专题
频谱分析
Spectralestimation（谱估计）的目标就是基于一个有限的数据集合描述一个信号的功率（在频率上的）分布。功率谱估计在很多场合下都就是有用的，包括对宽带噪声湮没下的
信号的检测 N点PSD估计
FFT。实践上大多数周期
2
Xl fk
fsL
,fk
kfs
s,k 0,1,K ,N 1
N
其中
Xl fk
Xl n 0
2 jkn/N
n e
选才iN就是大于L的下一个2的哥次就是明智的，要计算Xlfk我们直接对xLn补
零到长度为No假如L>N,在计算Xlfk前,我们必须绕回Xln模No
作为一个例子，考虑下面1001元素信号Xn,它包含了2个正弦信号与噪声randn('state',0);
fs = 1000;
t = (0:fs)/fs;
A = [1 2];
f = [150;140];
%Samplingfrequency
%Onesecondworthofsamples
%Sinusoidamplitudes(rowvector)
%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')
0
-10
Power Spectral Density Estimate via Periodogram
o O 2 3 - -
o O 4 5 - -
O
6 -
-80

Frequency(kHz)
平均功率通过用下述求与去近似积分求得
[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
周期图性能
卜面从四个角度讨论周期图法估计的性能：泄漏，分辨率，偏差与方差。
频谱泄漏
考虑有限长信号xLn，把它表示成无限长序列xn乘以一个有限长矩形窗wRn的
乘积的形式经常很有用
xLnxnwRn
因为时域的乘积等效于频域的卷积，所以上式的傅立叶变换就是
fs/2
・1s•
XLfXWRf
f
sfs/2
前文中导出的表达式
Px f
Xl f fsL
说明卷积对周期图有影响。
正弦数据的卷积影响最容易理解。假设
x n就是M个复正弦的与
M
xnAkejkn
k1
其频谱就是
M
XffsAkffk
k1
对一个有限长序列
，就变成了
Xl f
f fs/2 M
；fs A
fs fs/2 k 1
fk Wr f
M
d AW f fk k 1
所以在有限长信号的频谱中，Dirac函数被替换成了形式为WRffk的项，该项对应
于矩形窗的中心在fk的频率响应。
一个矩形窗的频率响应形状就是一个
sinc信号，如下所示
zuprtawBa D§1
，所以数据越短它的频谱泄
randn('state',0) fs = 1000;
t = (0:fs/10)/fs;
A = [1 2];
f = [150;140];
该图显示了一个主瓣与若干旁瓣，最大旁瓣大约在主瓣下方13、5dB处。这些旁瓣说
明了频谱泄漏效应。无限长信号的功率严格的集中在离散频率点fk处,而有限长信号在离散
频率点fk附近有连续的功率。
因为矩形窗越短，它的频率响应对Dirac冲击的近似性越差
漏越明显。考虑下面的100个采样的序列
%Samplingfrequency
%One-tenthofasecondworthofsamples
%Sinusoidamplitudes
%Sinusoidfrequencies