二：周期图法

经典谱估计中的周期图法是用得较多且最具代表性的方法。我们先复习基本的周期图法，接着针对它谱分辨率比较低等缺点，利用 Matlab 实现了几种改进的周期图法。

1、基本的周期图法

基本的周期图法可以提高计算效率，不需要计算自相关函数，但谱分辨率较低。基本周期图的基本原理是对观测到的数据直接进行傅立叶变换，然后取模的平方就是功率谱。取平稳随机信号 x(n)的有限个观察点 x(0)、x(1)、…x(n-1)，则傅立叶变换

matlab实现代码如下

clc;

clear;

t = 0:0.001:0.6;

y = sin(2*pi*50*t)+sin(2*pi*120*t);

x = y + randn(size(t));

N=length(x)；%观察数据点数

nfft=1024;%傅里叶变换点数

Xk=fft(x,nfft);

Pxx=abs(Xk).^2/N;

index=0:round(nfft/2-1);

k=index*N/nfft；% 频率（Hz）

plot_Pxx=10*log10（Pxx（index+1））；%功率谱 dB

plot(k,plot_Pxx);

xlabel('Hz');

ylabel（'功率谱 dB'）；

title（'用基本周期图法（直接计算）计算功率谱'）；

grid on;

2、Welch 法

P.O.韦尔奇提出一种把加窗处理与平均处理结合起来的方法。先把分段的数据乘以窗函数（进行加窗处理），分别计算其周期图，然后进行平均，得到的功率谱值为：

为了得到较好的功率谱估值，加窗和平均处理均应兼顾减小随机起伏和保证有足够的谱分辨率两个方面。

matlab程序

clc;

clear;

t = 0:0.001:0.6；

y = sin（2*pi*50*t）+sin（2*pi*120*t）；%信号

x = y+randn（size（t））；

N=length（x）；%观察数据点数

n=1:1:N；

%参数设置

Fs=1000；%每百年采样 100 次

nfft=512；%傅里叶变换点数

%window=hamming（25）；

window=boxcar（100）；

noverlap=20

range='half'；

%计算序列的 psd

[Pxx，f]=pwelch（x，window，noverlap，nfft，Fs，range）；

plot_Pxx=10*log10（Pxx）；%功率谱 dB

plot（f，plot_Pxx）

xlabel（'Hz'）；

ylabel（'功率谱 dB'）；

title（'用 welch 法计算的功率谱'）；

grid on；

3、多窗口的周期图法

由于普通的周期图功率谱估计只利用单一窗口，因此在序列始端和末端均会丢失相关信息。多窗口法利用多个正交窗口获得各自独立的近似功率谱估计，然后综合这些估计，最终得到平稳信号序列的功率谱估计。

Matlab 程序如下：

%信号

t = 0:0.001:0.6；

y = sin（2*pi*50*t）+sin（2*pi*120*t）；

x = y+randn（size（t））；

N=length（x）；%观察数据点数

n=1:1:N；

%参数设置

Fs=1000；%每百年采样 100 次

nfft=512；%傅里叶变换点数

%window=hamming（25）；

window=boxcar（100）；

noverlap=20

range='half'；

%计算序列的 psd

nw=3.5；p=0.99；

[Pxx，Pxxc，f]=pmtm（x，nw，nfft，Fs，p）；

plot_Pxx=10*log10（Pxx）；%功率谱 dB

plot（f，plot_Pxx）

xlabel（'Hz'）；

ylabel（'功率谱 dB'）；

title（'用多窗口法计算的功率谱'）；

grid on；

总结

通过实验可以看出，基本的周期图法谱分辨率较低，会有很大的失真。 Welch 法可以减小随机起伏，收敛性较好，曲线平滑，估计的结果方差较小，但是功率谱主瓣较宽，分辨率低。但是，如果信号序列不是足够长，由于每段序列长度变短，功率谱估值对不同频率成分的分辨能力也随之下降。多窗口的周期图法功率谱比 Welch 法的功率谱谱峰窄一些，谱峰高度也有所增大，因此多窗口的周期图法的功率谱表现得更加准确。普通周期图法的功率谱估计利用单一窗口，所以在序列的始端和末端均会丢失相关信息，而多窗口的周期图法以增加窗口来充分利用这些相关信息。