Python-用Burg算法估计随机信号的频率并求估计的均方误差

估计频率：

import math
import numpy as np
import matplotlib.pyplot as plt
import scipy.signal
import statsmodels.tsa.api as smt
pi = math.pi


#添加噪声的函数
def awgn(x,snr):
    snr=10**(snr/10)  #转化为单位为1
    xp=np.sum(x**2)
    nop=xp/snr   #计算噪声能量
    noise=np.random.randn(len(x)) #得出噪声序列
    pnoise=sum(noise**2)  #计算噪声能量
    n_noise=math.sqrt(pnoise)  #用于计算噪声值
    noise=noise/n_noise   #将噪声的能量变为1
    noise=noise*(math.sqrt(nop))   #得出噪声序列
    xr=x+noise
    return xr


N=128
n = np.arange(0,N,1)
x=10*np.cos(2*pi*0.1*n+pi/3)+2*np.cos(2*pi*0.11*n+pi/5)
xr = awgn(x,10)
p=40      #阶次不一定是40，由实际情况决定

ef=np.zeros((p+1,N))  #前向预测误差
eb=np.zeros((p+1,N))  #后向预测误差
k=np.zeros((1,p+1))   #反射系数,1
a=np.zeros((p+1,p+1))  #第m阶的m个参数,1,1
pp=np.zeros((1,p+1))  #预测误差功率,0
pp[0,0]=sum(xr**2)

ef[0,:]=xr
eb[0,:]=xr


for m in range(1,p+1):

    fenzi=0
    fenmu=0
    for n in range(m,N):
        fenzi=fenzi+ef[m-1,n]*eb[m-1,n-1]
        fenmu=fenmu+ef[m-1,n]**2+eb[m-1,n-1]**2

    k[0,m]=-2*fenzi/fenmu
    a[m,m]=k[0,m]

    ebyou=list(eb[m-1])
    ebyou=[0]+ebyou[0:N-1]
    ebyou=np.array(ebyou)
    ef[m,:]=ef[m-1,:]+k[0,m]*ebyou
    eb[m,:]=ebyou+k[0,m]*ef[m-1,:]

    #求a(m,)
    for i in range(1,m):
        a[m,i]=a[m-1,i]+k[0,m]*a[m-1,m-i]

    a[m,m]=k[0,m]

    #求pp()
    pp[0,m]=pp[0,m-1]*(1-(k[0,m])**2)
    if abs(pp[0,m]-pp[0,m-1])<0.001:
        break



G=math.sqrt(pp[0,m])
A1=a[m,1:m+1]
A=list(A1)
A.insert(0,1)

B=G
w,H=scipy.signal.freqz(B,A,worN=512)
H=abs(H)




ppn=np.arange(0,m+1,1)
w=w/(2*pi)
#plt画图是找不到字体，需要手动设置：
plt.rcParams['font.sans-serif']=['SimHei'] #显示中文标签
plt.rcParams['axes.unicode_minus']=False

plt.figure(1)
ax1=plt.subplot(2,1,1)
ax2=plt.subplot(2,1,2)

plt.sca(ax1)
plt.plot(ppn,pp[0,0:m+1])
plt.xlabel('阶次')
plt.ylabel('pp')
plt.title('误差功率')

plt.sca(ax2)
plt.plot(w,H)
plt.xlabel('f/Hz')
plt.ylabel('p')
plt.title('功率谱')
plt.show()




#估计频率
ff1=np.argmax(H)
H[ff1]=0
f1=ff1/(2*len(H))
for i in range(0,len(H)):
    ff2=np.argmax(H)
    f2=ff2/(2*len(H))
    if abs(f2-f1)<0.01:
        H[ff2] = 0
    else:
        break


print('f1的估计为：',f1)
print('f2的估计为：',f2)

均方误差：

import math
import numpy as np
import matplotlib.pyplot as plt
import scipy.signal
import statsmodels.tsa.api as smt
pi = math.pi

#添加噪声的函数
def awgn(x,snr):
    snr=10**(snr/10)  #转化为单位为1
    xp=np.sum(x**2)
    nop=xp/snr   #计算噪声能量
    noise=np.random.randn(len(x)) #得出噪声序列
    pnoise=sum(noise**2)  #计算噪声能量
    n_noise=math.sqrt(pnoise)  #用于计算噪声值
    noise=noise/n_noise   #将噪声的能量变为1
    noise=noise*(math.sqrt(nop))   #得出噪声序列
    xr=x+noise
    return xr

N=128
n = np.arange(0,N,1)
x=10*np.cos(2*pi*0.1*n+pi/3)+2*np.cos(2*pi*0.13*n+pi/5)

flag=0
p=40  #阶次不一定是40，由实际情况决定
num=512
mse1=[]    #用于记录f1估计的均方差
mse2=[]    #用于记录f2估计的均方差
c=np.zeros((1,100))    #存储每一个信噪比中100次f1的估计值
d=np.zeros((1,100))    #存储每一个信噪比中100次f1的估计值
e=np.zeros((1,100))    #为了求f1估计的方差
f=np.zeros((1,100))    #为了求f2估计的方差

for snr in range(-30,30,2):
    for i in range(0,100):
        xr=awgn(x,snr)
        ef = np.zeros((p + 1, N))  # 前向预测误差
        eb = np.zeros((p + 1, N))  # 后向预测误差
        k = np.zeros((1, p + 1))  # 反射系数
        a = np.zeros((p + 1, p + 1))  # 第m阶的m个参数
        pp = np.zeros((1, p + 1))  # 预测误差功率
        pp[0, 0] = sum(xr ** 2)

        ef[0, :] = xr
        eb[0, :] = xr

        for m in range(1, p + 1):

            fenzi = 0
            fenmu = 0
            for nn in range(m, N):
                fenzi=fenzi+ef[m-1,nn]*eb[m-1,nn-1]
                fenmu=fenmu+ef[m-1,nn]**2+eb[m-1,nn-1]**2

            k[0, m] = -2 * fenzi / fenmu
            a[m, m] = k[0, m]

            ebyou = list(eb[m - 1])
            ebyou = [0] + ebyou[0:N - 1]
            ebyou = np.array(ebyou)
            ef[m, :] = ef[m - 1, :] + k[0, m] * ebyou
            eb[m, :] = ebyou + k[0, m] * ef[m - 1, :]

            # 求a(m,)
            for index in range(1, m):
                a[m,index]=a[m-1,index]+k[0,m]*a[m-1,m-index]

            a[m, m] = k[0, m]

            # 求pp()
            pp[0, m] = pp[0, m - 1] * (1 - (k[0, m]) ** 2)
            if abs(pp[0, m] - pp[0, m - 1]) < 0.001:
                break

        G = math.sqrt(pp[0, m])
        A1 = a[m, 1:m + 1]
        A = list(A1)
        A.insert(0, 1)
        B = G
        w, H = scipy.signal.freqz(B, A, worN=512)
        H = abs(H)

        #下面估计频率f1,f2
        pxr = abs(H)  #执行后pxr是信号的功率谱
        ff1 = np.argmax(pxr)
        pxr[ff1] = 0
        if ff1==0:
            pxr[num-1]=0
        else:
            pxr[num - ff1] = 0
        f1 = ff1 / (2*(len(pxr)))
        for ii in range(0, len(pxr)):
            ff2 = np.argmax(pxr)
            f2 = ff2 / (2*(len(pxr)))
            if abs(f2 - f1) < 0.01:
                pxr[ff2] = 0
                if ff2==0:
                    pxr[num-1]=0
                else:
                    pxr[num-ff2]=0
            else:
                break
        c[0,i]=f1
        d[0,i]=f2

    ave1=sum(c[0,:])/100    #此信噪比下频率f1的均值
    ave2=sum(d[0,:])/100  #此信噪比下频率f2的均值
    for j in range(0,100):
        e[0,j]=(c[0,j]-ave1)**2
        f[0,j]=(d[0,j]-ave2)**2

    var1=sum(e[0,:])/100    #此信噪比下f1估计的方差
    var2=sum(f[0,:])/100  #此信噪比下f2估计的方差
    bia1=0.1-ave1   #f1估计的偏
    bia2=0.13-ave2   #f2估计的偏
    tem1=var1+bia1**2
    mse1.append(tem1)   #python索引从0开始
    tem2=var2+bia2**2
    mse2.append(tem2)
    flag=flag+1
    print(flag)

burg_mse1=mse1
burg_mse2=mse2

#下面画图
nf1=np.arange(0,2*len(mse1),2)
nf1=nf1-30
nf2=np.arange(0,2*len(mse2),2)
nf2=nf2-30

#plt画图是找不到字体，需要手动设置：
plt.rcParams['font.sans-serif']=['SimHei'] #显示中文标签
plt.rcParams['axes.unicode_minus']=False

plt.figure(1)
ax1=plt.subplot(2,1,1)
ax2=plt.subplot(2,1,2)

plt.sca(ax1)
plt.plot(nf1,mse1[0:])
plt.xlabel('snr/dB')
plt.ylabel('mse')
plt.title('burg_f1的均方差')

plt.sca(ax2)
plt.plot(nf2,mse2[0:])
plt.xlabel('snr/dB')
plt.ylabel('mse')
plt.title('burg_f2的均方差')
plt.show()

若有不当之处请多多指正