Python 功率谱估计直接法

xn=10sin(2pi0.1+pi/3)+4sin(2pi0.4+pi/4)，噪声为均值为零方差为一的高斯白噪声，用经典谱估计直接法来估计f1（0.1）和f2（0.4），并计算均方误差（MSE）。

原理

对N点DFT取绝对值后得到幅频特性，取绝对值平方后得到功率谱。
在估计f1、f2时，采取对幅频特性拟合后根据波峰处值最大来得到所要估计的f对应的k。
将整个程序循环多次即可计算估计值的期望和均方误差。
如何调整拟合效果见上篇博文。

结果

代码

import matplotlib.pyplot as plt
import numpy as np
from math import *
import random as ran
#定义Wn
def WN(N,k,n):
    return complex(cos(-2*pi/N*k*n), sin(-2*pi/N*k*n))
#整个程序循环多次来计算MSE
f_1 = 0
f_2 = 0
Error_f1 = 0
Error_f2 = 0
MSE_f1 = 0
MSE_f2 = 0
times = 5
for time in range(0, int(times)):
    # 设定初始参数
    f1 = 0.1
    f2 = 0.4
    N = 1024
    # 生成空列表存储结果
    Xk = []
    power = []
    #计算DFT和功率谱
    sum = 0
    for k in range(0, N):
        for n in range(0, N):
            xn = 10 * sin(2*pi*f1*n+3/pi)+4*sin(2*pi*f2*n+4/pi)+ran.gauss(0, 1)
            sum += xn * WN(N, k, n)
        Xk.append(sum)
        power.append(sum ** 2)
        sum = 0
    # 画图：频谱、功率谱（实际上为散点图，画为连续曲线）

    plt.title('spectrum map')
    f = [k / N for k in range(0, N)]
    dB_1 = [10*log10(i/N) for i in np.abs(Xk)]
    plt.plot(f, dB_1)
    plt.show()
    plt.title('power spectrum')
    dB_2 = [10*log10(i/N) for i in np.abs(power)]
    plt.plot(f, dB_2)
    plt.show()
    # 取DFT前半部分结果进行拟合
    for i in range(0, int(N / 2)):
        Xk.pop()
    Xk = list(np.abs(Xk))
    x = [i for i in range(0, int(N / 2))]
    y = np.polyfit(x, Xk, 6)
    Xk_smooth = np.polyval(y, x)
    # 对f1、f2进行估计得到对应的k值并存入列表F：波峰处值最大=对于任意的k，X(k)比两侧大即为波峰
    F = []
    for k in range(1, int(N / 2) - 1):
        if Xk_smooth[k] > Xk_smooth[k - 1] and Xk_smooth[k] > Xk_smooth[k + 1]:
            F.append(k)
    # 计算MSE
    f = [k/N for k in F]
    f_1 += f[0]
    f_2 += f[1]
    Error_f1 += (f[0]-f1)**2
    Error_f2 += (f[1]-f2)**2
f_1 = f_1/int(times)
f_2 = f_2/int(times)
print('f_1 = '+str(f_1))
print('f_2 = '+str(f_2))
MSE_f1 = Error_f1/int(times)
MSE_f2 = Error_f2/int(times)
print('MSE_f1 = '+str(MSE_f1))
print('MSE_f2 = '+str(MSE_f2))