代码：

import numpy as np
import math
import matplotlib.pyplot as plt
from scipy.fftpack import fft
pi = np.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

#求自相关函数的函数
def zxg(x):
    x=list(x)
    l=len(x)
    r=[0]*(l)
    xnp=np.array(x)
    xx=x
    r[0]=sum(xnp*xnp)
    for i in range(0,l-1):
        xx.insert(0,0)
        xxnp=np.array(xx[0:l])    #时延为几就在列表前加几个0，然后取前128个
        r[i+1]=sum(xnp*xxnp)
    tem=list(reversed(r))
    tem1=tem[0:l-1]
    rx=tem1+r     #rx为自相关函数
    return rx

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
num=255
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)
        rx = np.array((zxg(xr))) / N

        for index in range(0, N - 1):
            t = rx[1:255]
            tt = rx[0]
            rx[0:254] = t
            rx[254] = tt

        pxr = fft(rx, num)
        #下面估计频率f1,f2
        pxr = abs(pxr)  #执行后pxr是信号的功率谱
        ff1 = np.argmax(pxr)
        pxr[ff1] = 0
        if ff1==0:
            pxr[num-1]=0
        else:
            pxr[num - ff1] = 0
        f1 = ff1 / (len(pxr))
        for ii in range(0, len(pxr)):
            ff2 = np.argmax(pxr)
            f2 = ff2 / (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)


zxg_mse1=mse1
zxg_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('自相关的fft_f1的均方差')

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

如有不当之处请多多指正