Python 功率谱估计直接法（二）

做了什么

对上文的代码进行了改动，使其能在不同幅度的噪声条件下多次运行并计算同一噪声条件下均方误差的期望，而且在最后一次运行时输出最大噪声条件下的振幅谱和功率谱，并最终输出能反映均方误差（MSE）随信噪比（SNR）变化的曲线。其中f1、f2的估计也是在全部循环完成后取所有估计值的期望，因此估计效果较上文差。

结果

代码

import matplotlib.pyplot as plt
import numpy as np
from math import *
import random as ran
def WN(N,k,n):
    return complex(cos(-2*pi/N*k*n), sin(-2*pi/N*k*n))
a = 100
b = 90
times = 10
MSE_1 = []
MSE_2 = []
SNR = []
F1 = 0
F2 = 0
f_1 = 0
f_2 = 0
mse_1 = 0
mse_2 = 0
for c in range(1, int(times)+1):
    for theone in range(0,10):
        f1 = 0.1
        f2 = 0.4
        N = 512
        Xk = []
        power = []
        F = []
        P = []
        sum = 0
        for k in range(0, N):
            for n in range(0, N):
                xn = 100 * sin(2 * pi * f1 * n + 3 / pi) + 90 * sin(2 * pi * f2 * n + 4 / pi) + c * ran.gauss(0, 1)
                sum += xn * WN(N, k, n)
            Xk.append(sum)
            power.append(sum ** 2)
            sum = 0
        f = [k / N for k in range(0, N)]
        dB_1 = [10 * log10(i / N) for i in np.abs(Xk)]
        dB_2 = [10 * log10(i / N) for i in np.abs(power)]
        if c == int(times) - 1 and theone ==9:
            plt.subplot(2, 1, 1)
            plt.title('spectrum map')
            plt.plot(f, dB_1)
            plt.subplot(2, 1, 2)
            plt.title('power spectrum')
            plt.plot(f, dB_2)
            plt.show()
        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, 7)
        Xk_smooth = np.polyval(y, x)
        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)
                P.append(power[k])
        f = [k / N for k in F]
        f_1 += f[0]
        f_2 += f[1]
        mse_1 += (f_1-f1)**2
        mse_2 += (f_2-f2)**2
    f_1 = f_1 / 10
    F1 += f_1
    f_2 = f_2 / 10
    F2 += f_2
    MSE_1.append(mse_1/(10*times))
    MSE_2.append(mse_2/(10*times))
    SNR.append(10*log10((a**2+b**2)/2/c**2))
F1 = F1/times
F2 = F2/times
print('f1 = '+str(F1))
print('f2 = '+str(F2))
plt.suptitle('MSE_SNR')
plt.subplot(1,2,1)
plt.title('f1')
plt.plot(SNR,MSE_1)
plt.subplot(1,2,2)
plt.title('f2')
plt.plot(SNR, MSE_2)
plt.show()