估计频率：
import math
import numpy as np
import matplotlib.pyplot as plt
from scipy.linalg import toeplitz
from numpy.linalg import eig
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
#求自相关函数的函数
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=127
n = np.arange(0, N+1, 1)
x=10*np.cos(2*pi*0.1*n+pi/3)+10*np.cos(2*pi*0.11*n+pi/5)
xr = awgn(x, 10)
rx=zxg(xr)
rx=rx[127:255]
Rx=toeplitz(rx) #求出特征矩阵，即托步列兹阵
w, v = eig(Rx) #求特征分解，w为特征值，v为特征向量
w = list(sorted(w, reverse=True))
#下面求出信号的个数m
m = 0
for i in range(0, len(w)):
if w[i]/w[0] > 0.5:
m = m+1
else:
break
buchang=0.01
ww=np.arange(0,pi,buchang)
leng=len(ww)
p=[]
ll=1 #用于给p赋值
v=np.mat(v) #将特征向量v变为矩阵形式
vh=v.H #特征向量v的共轭转置
for i in range(0,leng):
ww=i*buchang #ww为频率
ei=[]
for k in range(0,N+1):
temp=complex(np.cos(k*ww),np.sin(k*ww))
ei.append(temp)
ei=np.mat(ei)
ei=ei.T #信号向量的转置
ei1=ei.H #信号向量的共轭转置
ei2=ei
summ=0
for n in range(m,N+1): #注意：python中的小特征值的特征向量从第3个开始，索引值为2
summ=summ+v[:,n]*vh[n,:]
pp=1/(ei1*summ*ei2)
pp=pp[0,0]
p.append(abs(pp))
#plt画图是找不到字体，需要手动设置：
plt.rcParams['font.sans-serif']=['SimHei'] #显示中文标签
plt.rcParams['axes.unicode_minus']=False
f=list(range(0,len(p)))
f=np.array(f)
f=f/(2*len(p))
f=list(f)
plt.plot(f,p)
plt.xlabel('f/Hz')
plt.ylabel('p')
plt.title('似功率谱')
plt.show()
p=np.array(p)
ff1=np.argmax(p)
p[ff1]=0
f1=ff1/(2*len(p))
ff2=np.argmax(p)
f2=ff2/(2*len(p))
print('f1的估计为：',f1)
print('f2的估计为：',f2)
均方误差：
import math
import numpy as np
import matplotlib.pyplot as plt
from scipy.linalg import toeplitz
from numpy.linalg import eig
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
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=127
n = np.arange(0,N+1,1)
x=np.cos(2*pi*0.1*n+pi/3)+np.cos(2*pi*0.11*n+pi/5)
flag=0
mse1=[] #用于记录f1估计的均方差
mse2=[] #用于记录f2估计的均方差
c=np.zeros((1,100)) #存储每一个信噪比中10次f1的估计值
d=np.zeros((1,100)) #存储每一个信噪比中10次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 = zxg(xr)
rx = rx[127:255]
Rx = toeplitz(rx) # 求出特征矩阵，即托步列兹阵
w, v = eig(Rx) # 求特征分解，w为特征值，v为特征向量
w = list(sorted(w, reverse=True))
# 下面求出信号的个数m
m = 0
for ii in range(0, len(w)):
if w[ii] / w[0] > 0.5:
m = m + 1
else:
break
buchang = 0.01
ww = np.arange(0, pi, buchang)
leng = len(ww)
p = []
ll = 1 # 用于给p赋值
v = np.mat(v) # 将特征向量v变为矩阵形式
vh = v.H # 特征向量v的共轭转置
for ii in range(0, leng):
ww = ii * buchang # ww为频率
ei = []
for k in range(0, N + 1):
temp = complex(np.cos(k * ww), np.sin(k * ww))
ei.append(temp)
ei = np.mat(ei)
ei = ei.T # 信号向量的转置
ei1 = ei.H # 信号向量的共轭转置
ei2 = ei
summ = 0
# 注意：python中的小特征值的特征向量从第3个开始，索引值为2
for n in range(m, N + 1):
summ = summ + v[:, n] * vh[n, :]
pp = 1 / (ei1 * summ * ei2)
pp = pp[0, 0]
p.append(abs(pp))
#下面估计频率f1,f2
num=len(p)
p=np.array(p)
pxr = abs(p) #执行后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
if (abs(f1-0.1)abs(f2-0.1)):
c[0,i]=f2
d[0,i]=f1
elif (abs(f1-0.1)>=abs(f1-0.13)) and (abs(f1-0.13)<=abs(f2-0.13)):
c[0,i]=f2
d[0,i]=f1
elif (abs(f1-0.1)>=abs(f1-0.13)) and (abs(f1-0.13)>abs(f2-0.13)):
c[0,i]=f1
d[0,i]=f2
else:
c[0,i]=f2
d[0,i]=f1
print(flag,i)
flag=flag+1
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)
te_mse1=mse1
te_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('特征分解_f1的均方差')
plt.sca(ax2)
plt.plot(nf2,mse2[0:])
plt.xlabel('snr/dB')
plt.ylabel('mse')
plt.title('特征分解_f2的均方差')
plt.show()
若有不当之处请多多指正
