文档介绍：TYYGROUP system office room 【TYYUA16H-TYY-TYYYUA8Q8-TYYUA162】
Matlab频谱分析程序
Matlab 信号处理工具箱 谱估计专题
频谱分析
Sp % Sinusoid frequencies (column vector)
xn = A*sin(2*pi*f*t) + *randn(size(t));
注意：最后三行表明了一个方便的表示正弦之和的方法，它等价于：
xn = sin(2*pi*150*t) + 2*sin(2*pi*140*t) + *randn(size(t));
对这个PSD的周期图估计可以通过产生一个周期图对象（periodogram object）来计算
Hs = ('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 =
你还可以用单边PSD去计算平均功率
[Pxxo,F] = psd(Hs,xn,fs,'onesided');
Pow = (fs/(2*length(Pxxo))) * sum(Pxxo)
Pow =
周期图性能
下面从四个角度讨论周期图法估计的性能：泄漏，分辨率，偏差和方差。
频谱泄漏
考虑有限长信号，把它表示成无限长序列乘以一个有限长矩形窗的乘积的形式经常很有用：
因为时域的乘积等效于频域的卷积，所以上式的傅立叶变换是
前文中导出的表达式
说明卷积对周期图有影响。
正弦数据的卷积影响最容易理解。假设是M个复正弦的和
其频谱是
对一个有限长序列，就变成了
所以在有限长信号的频谱中，Dirac函数被替换成了形式为的项，该项对应于矩形窗的中心在的频率响应。
一个矩形窗的频率响应形状是一个sinc信号，如下所示
该图显示了一个主瓣和若干旁瓣，最大旁瓣大约在主瓣下方处。这些旁瓣说明了频谱泄漏效应。无限长信号的功率严格的集中在离散频率点处，而有限长信号在离散频率点附近有连续的功率。
因为矩形窗越短，它的频率响应对Dirac冲击的近似性越差，所以数据越短它的频谱泄漏越明显。考虑下面的100个采样的序列
randn('state',0)
fs = 1000; % Sampling frequency
t = (0:fs/10)/fs; % One-tenth of a second worth of samples
A = [1 2]; % Sinusoid amplitudes
f = [150;140]; % Sinusoid frequencies
xn = A*sin(2*pi*f*t) + *randn(size(t));
Hs = ;
psd(Hs,xn,'Fs',fs,'NFFT',1024)
注意到频谱泄露只视数据长度而定。周期图确实只对有限数据样本进行计算，但是这和频谱泄露无关。
分辨率
分辨率指的是区分频谱特征的能力，是分析谱估计性能的关键概念。
要区分两个在频率上离得很近的正弦，要求两个频率差大于任何一个信号泄漏频谱的主瓣宽度。主瓣宽度定义为主瓣上峰值功率一半的点间的距离（3dB带宽）。该宽度近似等于
两个频率为的正弦信号，可分辨条件是
上例中频率间隔10Hz，数据长度要大于100抽才能使得周期图中两个频率可分辨。下图是只有67个数据长度的情况
randn('state',0)
fs = 1000; % Sampling frequency
t = (0:fs/15)./fs; % 67 samples
A = [1 2]; % Sinusoid amplitudes
f = [150;140]; % Sinusoid frequencies
xn = A*sin(2*pi*f*t) + *randn(size(t));
Hs=;
psd(Hs,xn,'Fs',fs,'NFFT',1024)
上述对分辨率的讨论都是在高信噪比的情况进行的，因此没有考虑噪声。当信噪比低的时候，谱特征的分辨更难，而且周期图上会出现一些噪声的伪像，如下所示
randn('state',0)
fs = 1000;