能量域特征提取方法

目录

1.经典谱估计方法

1.1直接法

1.2 间接法

1.3 Welch法

2. 现代谱估计方法

2.1 AR模型

（1）Burg算法

（2）MTM(Multiaper)方法

---

        所谓的功率谱，定义1：描述信号功率随频率变化的函数；定义2：以频率函数形式表达的信号或噪声频谱分量的幅度平方之半的分布。主要研究信号在频域中的各种特征，目的是根据有限数据在频域内提取被淹没在噪声中的有用信号。

        所谓谱估计或功率谱估计，就是用已观测到的一定数量的样本数据估计一个平稳随机信号的功率谱。这是因为功率谱在随机信号的分析和变换中起着类似于频谱在确定性信号分析中所起的作用，而且平稳随机过程统计特征的计算要求信号无限长，而实际上只能用一个样本，即有限长序列来计算，因此所得的计算值不是随机信号真正的统计值，而仅仅是一种估计。

        功率谱估计方法可以分为经典谱估计法与现代谱估计法。经典谱估计方法的基础是快速傅里叶变换(FFT)，因此该方法的计算效率很高。在一些实际问题中，当分辨率要求不是很髙时，经典谱估计方法就能很好地解决该问题。但是，该方法有着不可避免的缺点。当数据长度足够长时，经典谱估计方法的分辨率高，估计性能却较差。当数据很短时，经典谱估计方法的频率分辨率很低。

1.经典谱估计方法

        经典谱估计法又可分为直接法与间接法，直接法是利用快速傅里叶变换FFT算法对有限个样本数据进行傅里叶变换得到功率谱的方法，又称为周期图法；间接法是先得到样本数据的自相关函数估计，然后进行傅里叶变换得到功率谱的方法。由于直接法得到的功率谱估计存在方差大，而且不随点数N的增大而减小，谱曲线起伏大，或谱分辨率不高等缺点，于是又提出了几种改进算法：Bartlett法、Welch法及Nuttall法等。

1.1直接法

        直接法又称为周期图法，它是把随机信号的N点观察数据视为一能量有限信号，直接取的傅里叶变换，得，然后再取其幅值的平方，并除以N，作为对真实功率谱的估计。以表示用周期图法估计出的功率谱，则

temp=load('data');%导入自己数据
N=length(temp);   %数据点数
nfft=1024; %FFT所用数据长度
fs=150;%采样率

window=hamming(N);%汉明窗
[Pxx,f]=periodogram(temp,window,nfft,fs);

dataP=10*log10(Pxx);%功率谱

figure,plot(dataP)

 图1 周期图法功率谱

1.2 间接法

        直接法估计出的功率谱性能不好，当数据长度N太大时，谱线起伏加剧，N太小时，谱的分辨率又不好，因此需要加以改进。此处说的改进，主要是改进其方差特性。间接法是对直接法的一种改进，又称之为周期图法的平滑，此方法的基础是维纳-辛钦定理。

        间接法，又称自相关法，即BT法，1958年Blackman和Tukey给出了这一方法的具体实现，先由随机信号估计出其自相关函数，然后对求傅里叶变换得到的功率谱，记之为。

temp=load('data');%导入自己数据

fs=150;%采样率
nfft=1024;
corr=xcorr(temp,'unbiased');%计算序列的自相关函数
P=fft(corr,nfft);
Pxx=abs(P);
dataP=10*log10(Pxx);

figure,plot(dataP)

  图2 间接法功率谱

1.3 Welch法

        Welch功率谱算法是一种改进的周期图法，他在Bartlett法的基础上进行了两方面改进。首先，当信号序列分段时，该算法可将相邻两段重叠以提高方差性能，同时使用非矩形窗函数提升谱估计的分辨率，达到进一步减小方差的目的。Welch功率谱算法如下：

        一个随机信号序列,将其分割为L段，每段长度为M，相邻两段重合长度为M-k，则第i段的信号序列为：

其中，，L和M满足：。

        第i段信号的功率谱为：

式中U为归一化因子，为非矩形窗函数，则信号的Welch功率谱定义为：

temp=load('data');%导入数据

nfft=1024; 
L=256;   %数据段长度取值(64/128/256)
for i=1:1:length(L)
   sum_Pxx=0;
   for j=1:1:N/L(i)
       window=hamming(L(i));%汉明窗
       noverlap=length(window)*0.5;%重叠数
       [Pxx,f]=pwelch(temp,window,noverlap,nfft);
       sum_Pxx=sum_Pxx+abs(Pxx);
    end
 dataW=10*log10(sum_Pxx/N/L(i));
end

figure,plot(dataW)

   图3 Welch法功率谱

2. 现代谱估计方法

        现代谱估计方法，即参数化谱估计方法，是以1967年Burg提出的最大熵谱分析法为代表的谱估计方法。该方法是通过样本数据建立参数模型再按照求参数模型输出功率的方法估计该随机信号的功率谱，从而将功率谱估计问题转化为求参数问题。它并没有像经典谱估计方法假设的那样，不认为在观察到的N个数据以外的数据全为零，因此克服了经典谱估计分辨率低、方差性能不好以及不适合短样本等缺点，并且通过计算机模拟仿真和理论证实该谱估计方法是功率谱的有偏估计，谱估计的分辨率随样本数据序列长度的増加而提高。

        现代谱估计的线性参数模型主要包括三种，即AR模型（自回归模型），MA模型（滑动平均模型）和ARMA模型（自回归滑动平均模型），其中基于AR模型的功率谱估计是现代谱估计中最常用的一种方法，这是因为AR模型的参数可以通过解一组线性方程求得，计算简单。而对于MA模型和ARMA模型功率谱估计来说，其参数的计算需要解一组高阶的非线性方程，计算比较复杂。

2.1 AR模型

        AR模型又称为自回归模型，该随机信号是由本身的若干次过去值和激励的现实值的线性组合。

        对于p阶AR模型，有

其中为均值为零、方差为的白噪声，为p阶AR模型的参数。

        该AR模型系统的传递函数为

由上式可以看出，该系统的传递函数中只有极点，没有零点，所以AR模型是一个全极点的模型。

        因此，的功率谱为

        如果我们能得到AR模型的参数，代入上式，就能得到功率谱估计，从而将功率谱估计问题变成了AR模型阶数p和AR模型p+1个参数确定的问题。

        经过推导，可得

其中为的自相关函数。上式称为AR模型的正则化方程(Yule-Walker方程)。很容易得到，AR模型的正则方程是一组线性方程。

（1）Burg算法

        Burg算法是建立在样本数据的基础上的，是对AR模型参数求解的有效算法。该算法是自相关法的改进，使前、后向预测误差的平均功率最小来计算反射系数，以Levinson-Dubin递推为限制条件，求出随机信号的功率谱。

        Burg算法是一种迭代算法，通过不断迭代寻找自回归模型的系数，以逼近最小均方误差。它是一种无记忆AR模型，因此可以适用于实时处理，且计算效率高。Burg算法的优点是具有较高的精度和严格的数学基础，对于多变量预测和滤波具有很好的效果。但是缺点是对于高阶自回归模型，计算复杂度会很高，处理时需要考虑过程中的不稳定性问题。​

temp=load('data');%导入数据

p=8;%阶数的取值
nfft=1024;
for i=1:1:length(p)
   sum_Pxx=0;
   [Pxx,f]=pburg(temp,p(i),nfft);
   sum_Pxx=sum_Pxx+abs(Pxx);
   dataB=10*log10(sum_Pxx);
end

figure,plot(dataB)

   图4 Burg法功率谱

（2）MTM(Multiaper)方法

        多窗谱估计法（MTM）利用正交的窗口截取数据段并获取相互独立的功率谱估计，然后再结合这些结果以得到最终的估计。多窗谱估计法最重要的参数是时间与带宽乘积（设为）。这个参数将会直接决定谱估计窗的个数，这里窗的个数为个。由此可以得知，NW的增大，窗的个数也会迅速增多，这样就会有更多的谱估计，那么谱估计的方差便得到了减小。但是问题也会随之而来，那就是谱泄漏的增大，并且正的谱估计的结果之偏差更为严重。因此，NW的选择将会决定在使用本方法进行功率谱估计时结果的好坏，这就要仔细考虑偏差和方差之间的平衡性问题了。

temp=load('data');%导入数据

fs=150;%采样率
[Pxx,f]=pmtm(temp,4,nfft,fs);
f=f/max(f)/2;%归一化
dataM=10*log10(Pxx);

figure,plot(dataM)

   图5 MTM法功率谱