用Python编程实现功率谱估计的平滑改进

周期图法求解信号的功率谱

周期图法求解信号功率谱的基本思想是：将随机信号x(n)的N点观测数据视为能量有限信号，对该信号求解傅里叶变换，求解变换结果幅值的平方，并且除以N。
这一部分的代码如下：

#sig为所给信号，N为信号的长度
Pper=np.zeros(N)
xfft=np.fft.fft(sig)
for i in np.arange(N):
    Pper[i]=np.abs(xfft[i])*np.abs(xfft[i])/N

试验时所用的信号为：

N=200
sig=np.random.randn(N,)

程序运行结果为：

对周期图法计算出的功率谱进行平滑改进

在本实验中，利用了Welch法对周期图法进行了改进，其基本思想是：对原始信号x(n)进行分段，并且允许每一段的数据有部分交叠；此外Welch法还允许选择不同的窗函数来进行分段。本文利用了最简单的窗函数来进行分段。
在利用Welch法时，首先要计算信号的分段数目，计算信号分段数的程序为：

def segL(N,M, overlap):     #N：信号长度；M:分段长度；overlap:重叠长度；返回值为分段数目sn。
    d=M-overlap
    n=np.int32(np.ceil((N-overlap)/d))
    return n

在计算完分段数后，就可以构造一个由所有段信号组成的sn*M的xi矩阵，对该矩阵的每一行再进行加窗操作，窗函数的宽度为该矩阵的列数。
加窗函数的程序为：

for i in np.arange(sn):
    xi=xi_snM[i]                     #xi_snM为求出的sn*M的xi矩阵
    for j in np.arange(M):
        xid[j]=xi[j]*d[j]            #d为所用的窗函数
    xifft[i]=np.fft.fft(xid)         #计算加窗后信号的傅里叶变换

到此为止，计算出了snM矩阵的每一行信号的傅里叶变换，其中xifft也是一个snM的矩阵。利用Welch法计算功率谱时，需要对xifft进行变换，变换的目的主要是让xifft的各行之间按照分段时的原则进行每段功率谱的相加求和。本文用的M=50，overlap=25。计算出的每小段开始时的索引号为：

计算完成的xifft和加窗后并且进行索引号变换的的xifftwin的数据如下图所示：


原始信号和各小段信号组合后的新信号的图形如下图所示：

最后就算xifftwin每一列幅值的平方和，并且除以MUsn，就可以得到改进后的功率谱。其中U为归一化因子，计算的程序为：

UWinwd=0
for n in np.arange(Winwd):      #Winwd为窗函数的宽度，此处Winwd=M。
    UWinwd=d[n]*d[n]+UWinwd     #d为所选用的窗函数，本文选用的为矩形窗，计算出的归一化因子为1。 
U=UWinwd/Winwd    

计算改进后的功率谱的程序为：

Pperwin=np.zeros(Nex)
for j in  np.arange(Nex):
    for i in np.arange(sn):
        Pperwin[j]=np.abs(xifftwin[i,j])*np.abs(xifftwin[i,j])+Pperwin[j]
    Pperwin[j]=Pperwin[j]/(M*U*sn)

最终得到的结果为：

如果本文有什么错误，还请读者谅解，并且能够帮忙指正，感激不尽！！！