现代信号处理-经典功率谱估计＜改进前＞

样本函数的总能量是无限的，但是其平均功率是有限值。估计的意义是取一段有限的信号值，去反映原信号的波形或者统计特性。然而对于随机信号而言，它没有确切信号，因此直接研究其频谱没有意义，所以才回去研究平均功率随频率的分布，也就是功率谱估计。
功率谱密度有很多方法，本篇是从经典谱估计出发。

经典功率谱估计共分为两大部分：第一部分是普遍用的几种方法：周期图法（Periodogram，直接法）、BT法（Blankman-Tukey，自相关法，间接法）以及BT法的两种情况。第二部分是利用两种分段方法（Bartlett法和Welch法）改进功率谱估计。

---

前言

在周期图法和BT法中，功率谱估计的原信号设为高斯白噪声，因为BT法必须用在广义平稳随机过程，所以现有的函数只有周期图法的函数，所以前半部分中，两种方法我都用单一样本表示，改进前的matlab程序是从原始公式开始写，这样可以更方便地了解两种方法的内部过程。最后也会给出功率谱估计的质量：偏差和方差。

而在改进方法中，我们直接用Pwelch函数来实现分段的周期图法；而由于没有直接的分段的BT法公式（可能是因为BT法只能用在广义平稳过程中），所以在单一样本的基础上，先进行分段处理，每段再各自利用BT法处理，最后再求平均，来实现分段的BT法。第三个方法是利用MTM法来实现分辨率和方差的平衡。

---

提示：以下是本篇文章正文内容，下面案例可供参考，全文的几个主程序是一个整体，请勿分开复制直接运行。

一、功率谱估计方法（改进前）

单一样本且不分段

由于高斯白噪声具有各态历经性，所以可以用足够长的单一样本来估计，用时间平均来代替集总平均，所以下面两种方法都是用单一样本来估计。两种方法的表达式如下：

第一个式子是BT法，也就是利用维纳-辛钦定理，先求出原信号的自相关函数，再求傅里叶得到功率谱密度，所以也称这种方法为直接法。但由于此定理只能用在广义平稳随机过程中，所以BT法具有很大的局限性。
第二个式子是周期图法，是从能量和的角度出发，用傅里叶系数的平方来表示，这是最早提出的功率谱估计的方法，沿用至今，也被称为直接法。

1.周期图法

虽然这种方法最早出现，但现在相比以前，多了FFT来做傅里叶变换，加快了运算速度。公式如下表示：
下面是此方法实现的代码和结果，由于理论上周期图法和BT法是等价的，在表达式中也可以相互推导，只有点数不同：周期法补零为2*N-1时，与BT法M=N-1完全等效。所以在周期图法中给出了补零值的情况，以便后面与BT法M=N-1的情况作比较。

clc;clear all;close all;
N=200;
rx_BT=zeros(1,N);
x=randn(1,N);
%% 周期图解法PER
xk=abs(fft(x));       %%做的是N1点fft
for ki=1:N
    sf_per(ki)=((xk(ki))^2)/N;
end
rx_per=fftshift(abs(ifft(sf_per)));
figure;
subplot(2,2,1);
plot(rx_per);title('per的自相关函数');
subplot(2,2,2);
plot(sf_per);title('per的功率谱密度');
axis([-inf +inf -5 5]);

%% 周期图解法PER,XN补零成X2N-1
xn_2n=[x,zeros(1,N-1)];  %%补零
xk_2n=abs(fft(xn_2n));   %%做的是2*N1点fft
    for ki=1:(2*N-1)
        sf_per_2n(ki)=((xk_2n(ki))^2)/(2*N-1);
    end
rx_per_2n=fftshift(abs(ifft(sf_per_2n)));
subplot(2,2,3);
plot(rx_per_2n);title('per的自相关函数,XN插零成X2N');
subplot(2,2,4);
plot(sf_per_2n);title('per的功率谱密度');
axis([-inf +inf -5 5]);

上面周期图法的程序是用FFT写的，也可以用periodogram函数来写：

[Pxx,f] = periodogram(x,window,nfft,fs)

至此，也可以证明一下，BT法M=2N-1等效于周期图法的2N-1点，由于周期图法补完零，功率减半，这里为了与BT法对比，将其结果乘2。

plot(lags,sf_BT,lags,sf_per_2n*2);title('证明BT法等效周期图法');
legend('BT法M=2N-1','周期图法2N-1点(结果*2）');

2.BT法

BT法虽然仅用于广义平稳过程，但是它的优点是可以在自相关函数上二次截短加窗，使功率谱更加平滑。
例如在高斯白噪声中，自相关函数理论上是冲激函数，而加窗就相当于使其逼近理论函数，功率谱也更加逼近常数1。
加窗的实现是下式中的M取值，当M=N-1时等效于周期图法，而当M<N-1时，也就意味使自相关函数加窗后再做傅里叶变换。

3.周期图法和BT法的关系

根据上述所述，需要考虑BT法下M的取值：

下面讨论了M的三种方法：M=N-1;M=N/4;M=N/8

%% 间接法BT法M=N-1
[rx_BT,lags]=xcorr(x,'coeff');
sf_BT=(abs(fft(rx_BT)));
figure;
subplot(3,2,1);
plot(lags,rx_BT);title('BT的自相关函数,M=N-1');
subplot(3,2,2);
plot(lags,sf_BT);title('BT的功率谱密度');
axis([-inf +inf -5 5]);

%% 间接法BT法M<N-1
M1=N/4;
rx_BT_M1=zeros(1,2*N-1);
rx_BT_M1(N)=rx_BT(N);
for ii=1:M1
    rx_BT_M1(N+ii)=rx_BT(N+ii);
    rx_BT_M1(N-ii)=rx_BT(N-ii);
end
sf_BT_M1=(abs(fft(rx_BT_M1)));
subplot(3,2,3);
plot(lags,rx_BT_M1);title('BT的自相关函数,M=N/4<N-1');
subplot(3,2,4);
plot(lags,sf_BT_M1);title('BT的功率谱密度');
axis([-inf +inf -5 5]);

M2=N/8;
rx_BT_M2=zeros(1,2*N-1);
rx_BT_M2(N)=rx_BT(N);
for ii=1:M2
    rx_BT_M2(N+ii)=rx_BT(N+ii);
    rx_BT_M2(N-ii)=rx_BT(N-ii);
end
sf_BT_M2=(abs(fft(rx_BT_M2)));
subplot(3,2,5);
plot(lags,rx_BT_M2);title('BT的自相关函数,M=N/8<N-1');
subplot(3,2,6);
plot(lags,sf_BT_M2);title('BT的功率谱密度');
axis([-inf +inf -5 5]);

值得注意的是：为了对比实验组，加完窗做的FFT还是加窗前的点数N。如果想要去掉加完窗两边的零点，再做FFT，就需要改变自相关函数[rx,lags]=xcorr(x,y,maxlag,‘coeff’)里的maxlag参数，这里的maxlag（最大滞后）的意思是将rx加上矩形窗，再将加窗后的自相关函数做FFT得到功率谱密度。
补零和不补零得到的功率谱密度函数是相同的，只是采样率有所不同，这里对比了补零和不补零做FFT的区别：

4.功率谱估计的质量

一般用偏差和方差来估计功率谱估计的好坏，二者的求解如下所示。本文依照此公式，求得上述五种角度的功率谱估计得出各自的估计质量。

按照上式，求五种角度的偏差和方差

L1=length(sf_BT);
L2=length(sf_BT_M1);
L3=length(sf_BT_M2);
L4=length(sf_per);
L5=length(sf_per_2n);
L=[L1 L2 L3 L4 L5];
LMAX=max(L);
SF=zeros(length(L),LMAX);
for k=1:L(1)
    SF(1,k)=sf_BT(k);
end
for k=1:L(2)
    SF(2,k)=sf_BT_M1(k);
end
for k=1:L(3)
    SF(3,k)=sf_BT_M2(k);
end
for k=1:L(4)
    SF(4,k)=sf_per(k);
end
for k=1:L(5)
    SF(5,k)=sf_per_2n(k);
end
%% 求几种方法的偏差
for i=1:length(L)
    SF_E(i)=sum(SF(i,:))/L(i);
    bia(i)=SF_E(i)-1;   %%E(P(K))-P(K)  功率谱密度理论值是常数1，因为是高斯白噪声
end   
figure;
subplot(2,1,1);
bar(abs(bia),0.5);axis([0 length(L)+1 -inf inf]); 
for i=1:length(bia)
    text(i,bia(i),[num2str(bia(i))],'HorizontalAlignment','center');
end
set(gca,'XTick',1:length(L));
set(gca,'XTickLabel',{'BT法不取窗','BT法取M=N/4窗','BT法取M=N/8窗','周期图法N点','周期图法2N点'});
title([num2str(length(L)),'种方法的偏差']);

%% 求几种方法的方差
for ii=1:length(L)
    for kk=1:L(ii)
        SF_2(ii,kk)=(SF(ii,kk)-SF_E(ii))^2;
    end
    var(ii)=sum(SF_2(ii,:))/L(ii);
end
subplot(2,1,2);
bar(var,0.5);axis([0 length(L)+1 -2 2]); 
for i=1:length(var)
    text(i,var(i),[num2str(var(i))],'HorizontalAlignment','center');
end
set(gca,'XTick',1:length(L));
set(gca,'XTickLabel',{'BT法不取窗','BT法取M=N/4窗','BT法取M=N/8窗','周期图法N点','周期图法2N点'});
title([num2str(length(L)),'种方法的方差']);

可以得出结论：
1.方差和分辨率成反比
2.增大数据N，方差会减小
3.取窗越小，方差也越小
4.偏差和方差没有必然关系

二、功率谱估计的改进

1.Pwelch法（周期图法下的分段方法）

2.BT法下的分段方法

3.MTM法

由于篇幅问题，改进后的功率谱估计留着下一节

---

总结

今天从基础公式出发，介绍和仿真了两种功率谱估计的方法。但是由于这种情况的方差很大，功率谱密度的波形也很粗糙，所以改进主要是从分段出发，分段有两种方式：也就是不重叠分段（Bartlett 平均）和重叠分段（Welch 平均）。又由于功率谱密度的分辨率和方差是一对矛盾，所以MTM法主要是解决二者的平衡关系。