利用经典谱估计法估计信号的功率谱(随机信号)

随机信号

利用经典谱估计法估计信号的功率谱

作业综述：

给出一段信号“asd.wav”,利用经典谱估计法的原理，通过不同的谱估计方法，求出信号的功率谱密度函数。采用MATLAB语言，利用MATLAB语言强大的数据处理和数据可视化能力，通过GUI的对话框模板，使操作更为简便！在一个GUI界面中，同时呈现出不同方法产生出的功率谱。

这里给出了几种不同的方法：BT法，周期图法，平均法以及Welch法。把几种不同方法所得到的功率谱都呈现在一个界面中，便于对几种不同方法得到的功率谱作对比。

一．题目要求

给出一段信号及采样率，利用经典谱估计法估计出信号的功率谱。

二．基本原理及方法

经典谱估计的方法，实质上依赖于传统的傅里叶变换法。它是将数据工作区外的未知数据假设为零，相当于数据加窗，主要方法有BT法，周期图法，平均法以及Welch法。

1. BT法(Blackman-Tukey)

●理论基础：

(1)随机序列的维纳-辛钦定理

由于随机序列{X(n)}的自相关函数Rx(m)=E[X(n)X(n+m)]定义在离散点m上，设取样间隔为，则可将随机序列的自相关函数用连续时间函数表示为

等式两边取傅里叶变换，则随机序列的功率谱密度

(2)谱估计

BT法是先估计自相关函数Rx(m)(|m|=0,1,2…,N-1),然后再经过离散傅里叶变换求的功率谱密度的估值。即

其中可有式得到。

2. 周期图法

●理论基础：

周期图法是根据各态历经随机过程功率谱的定义来进行谱估计的。在前面我们已知，各态历经的连续随机过程的功率谱密度满足

式中 是连续随机过程第i 个样本的截取函数 的频谱。对应在随机序列中则有

由于随机序列中观测数据 仅在 的点上存在，则 的N 点离散傅里叶变换为：

因此有随机信号的观测数据 的功率谱估计值（称“周期图”）如下：

由于上式中的离散傅里叶变换可以用快速傅里叶变换计算，因此就可以估计出功率

谱。

3.平均法：

理论基础：

平均法可视为周期图法的改进。周期图经过平均后会使它的方差减少，达到一致估计的目的，有一个定理：如果 ， ， ， 是不相关的随机变量，且都有个均值 及其方差 ，则可以证明它们的算术平均的均值为 ，方差为

。

由定理可见：具有 个独立同分布随机变量平均的方差，是单个随机变量方差的 ，

当 时，方差

，可以达到一致估计的目的。因此，将 个独立的估计量经过算术

平均后得到的估计量的方差也是原估计量方差的 。

平均图法即是将数据 ， ， 分段求周期图法后再平均。例如，给定N=1000个数据样本（平均法适用于数据量大的场合），则可以将它分成10个长度为100的小段，分别计算每一段的周期图

()()2

1001100,100(1)

1

,1,2,```,10100

l j l n l G w X e l ω-=-=

=∑

然后将这10个周期图加以平均得谱估计值：

()()

10

100100,1

110l l G w G w ==∑

由于这10小段的周期图取决于同一个过程，因而其均值相同。若这10个小段的周期图是统计独立的，则这10个小段平均之后的方差却是单段方差的 。

D[ ]= D[ ]

即：平均法将 ， ， 的N 个数据分成L 段（N=ML ），若各数据段相互独立，则平方后估计量的方差是原来不分段估计量方差的 。所以当 时，估计量的方差趋于0，达到一致估计的目的。但是，随着分段数L 的增加，M 点数减少，分辨率减少，使估计变成有偏估计。相反，若L 减少，M 增加，虽偏差减少，但方差增大。所以，在应用中，必须兼顾分辨率和方差的要求来适当选择M 和L 的值。

4.Welch 法：

理论基础：

Welch 法又称修正周期法，其步骤为：

○

1先将N 个数据分成L 段，每数据段M 个数据，N=ML ； ○

2选择适当的窗函数w （n ），并用该w （n ）依次对每段数据做相应的加权，然后确定每段的周期图

()2

1,(1)

1

(),1,2,```,Ml j n M l n n M l G w x w n e l L MU

ω--=-=

=∑

U 为归一化因子

1

20

1(),M n U w n M

-==

∑

对每段周期图进行平均得到功率谱估计：

()(),1

1L

M M l l G w G w L ==∑

当数据量一定时，若分段数L 增加，M 点数减少，则分辨率下降；若L 减少，虽M 增加，但方差增大。解决这一矛盾的方法是，让数据段间适当“重叠”。

三．算法设计与实现

1程序流程图

（1） BT 法的流程图

（2）周期图法的流程图

（3）平均法的流程图

（4）Welch法的流程图

2.主要模块的设计:

(1)产生原始信号

Fs=600;

nfft=1024;

n=0:1/Fs:1;

y=cos(2*pi*30*n)+3*cos(2*pi*100*n)+randn(size(n));

axes(handles.axes1);

plot(n,y);

title(‘原始信号波形‘);

原始信号波形如下：