一、经典功率谱估计方法

对于离散时间平稳随机过程 $u(n)$，它的功率谱 $S(w)$ 描述了随机过程 $u(n)$ 中各频率成分的平均功率的大小。因此，可以通过计算功率谱密度函数 $S(w)$ 来间接的了解随机过程 $u(n)$ 中各频率成分的构成情况。
根据维纳-辛钦定理，平稳随机过程 $u(n)$ 的自相关函数 $r(m)$ 与其功率谱 $S(w)$ 是一对傅里叶变换关系。利用相关函数的傅里叶变换来估计随机过程功率谱的方法，称为经典功率谱估计。

1、周期图法

由于随机过程 $u(n)$ 的 $N$ 个观测值 $u_N(n)$ 是确定信号，对其进行傅里叶变换，得：
$$U_N(w)=\sum _{n=0}^{N-1}{u}_N(n)e^{-jwn}$$
根据帕斯瓦尔关系，上式的模的平方是确定信号 $u_N(n)$ 的能量谱，对能量谱除以持续时间 $N$，其结果则为 $u_N(n)$ 的功率谱估计，将其作为随机过程 $u(n)$ 的功率谱的估计，表示为：
$$S_{per}(w)=\frac{1}{N}|U_N(w){|}^2$$
该方法称为功率谱估计的周期图法(periodogram)。该方法由于是直接通过观测数据的傅里叶变换求得，其又被称为直接法。

2、BT法

1958年，Blackman和Tukey在维纳-辛钦定理的基础上，提出了自相关谱估计法(简称BT法)。其具体实现方法为：先由观测数据估计出自相关函数，然后对其求傅里叶变换，以此结果作为对功率谱的估计。
该方法的计算步骤为以下两个公式：
$$r(m)=\frac{1}{N}\sum _{n=0}^{N-1}{u}_N(n)u_N^{\ast }(n-m)$$
$$S_{BT}(w)=\sum _{m=-M}^{M}r(m)e^{-jwm}，0\le M\le N-1$$
上述方法是通过自相关函数间接得到的，因此又被称为间接法。
当 $M=N-1$ 时，上式的功率谱估计可表示为：
$$S_{BT}(w)=\frac{1}{N}|U_N(w){|}^2$$
比较周期图法和BT法可以看出，周期图法是BT法的一个特例。当 $M=N-1$ 时，周期图法与BT法相同；当 $M<<N-1$ 时，BT法是对周期图法的平滑。

3、Welch法

1967年Welch提出了修正平均周期图法，即Welch法，它是在Bartlett法的基础上改进的。
Welch法估计信号功率谱的计算步骤如下：
（1）将长度为 $N$ 的信号 $u_N(n)$ 进行分段，相邻的两段数据交叠一半。若每段信号的长度为 $M$，信号将被分成 $L$ 段，即：
$$L=\frac{N-M/2}{M/2}$$
（2）将第 $i(1\le i\le L)$ 段信号 $u_N^i(n)$ 与长度为 $M$ 的窗函数 $w(n)$ 相乘。
（3）对加窗后的每段信号，利用周期图法求得其功率谱，即：
$$S_{per}^i(w)=\frac{1}{M}|\sum _{n=0}^{M-1}{u}_N^i(n)w(n)e^{-jwn}{|}^2$$
（4）对每一段估计到的功率谱 $S_{per}^i(w)$ 进行求和平均，则得到Welch法的功率谱估计：
$${\stackrel{-}{S}}_{per}(w)=\frac{1}{L}\sum _{i=1}^L{S}_{per}^i(w)$$

二、仿真实例

设随机过程 $u_N(n)$ 由3个实正弦信号加噪声构成，即：
$$u_N(n)=\sum _{k=1}^3{s}_k(n)+v(n)$$
其中，$s_k(n)=A_{k}cos(2\pi f_{k}n+{\phi }_k)$，归一化频率分别为 $f_1=0.1，f_2=0.25，f_1=0.27$，${\phi }_k$ 是相互独立并在[0,2π]上服从均匀分布的随机相位；$v(n)$ 是均值为0、方差为1的实高斯白噪声序列。
3个正弦信号的幅度 $A_k>0$ ，由每个信号的信噪比 $SN{R}_k$ 决定。它们的关系为：
$$SN{R}_k=10lo{g}_{10}(\frac{A_k^2}{2{\sigma }^2})$$

1、性能比较

2、总结

（1）经典功率谱估计方法可以用FFT进行快速计算，且计算量较小，是目前较常用的谱估计方法；
（2）功率谱的分辨率较低，它正比于 $2\pi /M$，$M$ 是所使用信号的长度。
（3）由于加窗的影响，使得估计的功率谱主瓣展宽，降低了分辨率；
（4）方差性能不好，不是 $S(w)$ 的一致估计，且 $N$ 增大时，谱曲线的起伏加剧；
（5）周期图的平滑和平均与窗函数的使用密切相关。平滑和平均主要是改善了周期图的方差，但会使得分辨率降低和偏差增大。因此，在实际应用时，需要在方差、分辨率和偏差之间进行折中选择。

3、MATLAB代码

clc;
clear all;
close all;

%% 仿真信号
N = 256;     % 观测样本数(周期图法设为64，其他方法用256)
SNR = [30,30,27];
A = sqrt(2)*10.^(SNR./20);
freq = [10,25,27];     % 三个信号的频率
Fs = 100;   % 采样频率
t = (0:N-1)/Fs;   % 时间
fai = 2*pi*rand(3,length(t));
x = zeros(size(t));
for k = 1:3
    x = x+A(k)*cos(2*pi*freq(k)*t+fai(k,:));    % 3个正弦信号的叠加
end
vn = randn(size(t));
y = x+vn;   % 含有高斯白噪声的信号

%% 周期图法、BT法与Welch法估计随机信号的功率谱及性能对比
% 1、周期图法
N_fft = 256;
figure;
periodogram(y,[],'centered',N_fft,Fs);
xlabel('频率/Hz');
ylabel('功率谱密度(dB/Hz)');
title(['周期图法',',N=',num2str(N)]);

% 2、BT法
N_fft = 256;
rm = xcorr(y,'unbiased');
yn = fft(rm,N_fft);     % 对信号进行快速Fourier变换
Ys = fftshift(yn);
figure;
fshift = (-N/2:N/2-1)*(Fs/N);   % zero-centered frequency range
powershift = 20*log10(abs(Ys));       % zero-centered power

plot(fshift,powershift);
xlabel('频率/Hz');
ylabel('功率谱密度(dB/Hz)');
grid on;
title('BT法');

% 3、Welch法
N_fft = 256;
M = 64;         % 窗长
Overlap = 32;   % 重叠点的数目
figure;
pwelch(y,M,Overlap,N_fft,Fs,'centered');
xlabel('频率/Hz');
ylabel('功率谱密度(dB/Hz)');
title('Welch法');

三、参考文献

[1] 何子述, 夏威, 等. 现代数字信号处理及其应用[M]. 北京: 清华大学出版社，2009.

最后

以上就是老迟到鸭子为你收集整理的【信号频率估计】经典功率谱估计及其MATLAB仿真一、经典功率谱估计方法二、仿真实例三、参考文献的全部内容，希望文章能够帮你解决【信号频率估计】经典功率谱估计及其MATLAB仿真一、经典功率谱估计方法二、仿真实例三、参考文献所遇到的程序开发问题。

如果觉得靠谱客网站的内容还不错，欢迎将靠谱客网站推荐给程序员好友。