开发: C++知识库 Java知识库 JavaScript Python PHP知识库人工智能区块链大数据移动开发嵌入式开发工具数据结构与算法开发测试游戏开发网络协议系统运维
教程: HTML教程 CSS教程 JavaScript教程 Go语言教程 JQuery教程 VUE教程 VUE3教程 Bootstrap教程 SQL数据库教程 C语言教程 C++教程 Java教程 Python教程 Python3教程 C#教程
数码: 电脑笔记本显卡显示器固态硬盘硬盘耳机手机 iphone vivo oppo 小米华为单反装机图拉丁

-> 人工智能 -> 现代信号处理-经典功率谱估计＜改进前＞ -> 正文阅读

[人工智能]现代信号处理-经典功率谱估计＜改进前＞

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

经典功率谱估计共分为两大部分：第一部分是普遍用的几种方法：周期图法（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]);

BT法
值得注意的是：为了对比实验组，加完窗做的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法主要是解决二者的平衡关系。