谱峰搜索传播算子（python）

dongming2017

于 2021-12-23 17:20:18 发布

阅读量1.3k

点赞数 1

文章标签： python 线性代数机器学习

本文链接：https://blog.csdn.net/wudongming2017/article/details/122112470

版权

使用谱峰搜索传播算子的方法进行DOA估计

关键在于构造矩阵Q

import numpy as np
import scipy.signal as ss
import scipy.linalg as LA
import matplotlib.pyplot as plt

def awgn(x, snr):
    spower = np.sum((np.abs(x) ** 2)) / x.size
    x = x + np.sqrt(spower / snr) * (np.random.randn(x.shape[0], x.shape[1]) + 1j * np.random.randn(x.shape[0], x.shape[1]))
    return x

derad = np.pi / 180
radeg = 180 / np.pi
d = np.arange(0,4,0.5)
theta = np.array([0,30,60]).reshape(1,-1)
n=500
snr = 10
iwave = theta.size

A = np.exp(-1j*2*np.pi*d.reshape(-1,1)@np.sin(theta*derad))
S = np.random.randn(iwave,n)
X = A@S
X = awgn(X,snr)
Rxx = X@(X.conj().T)/n

K = Rxx.shape[0]
G = Rxx[:,0:iwave]
H = Rxx[:,iwave:K]
P = LA.pinv((G.conj().T)@G)@(G.conj().T)@H

Im = -np.diag(np.ones(K-iwave))
Q = np.concatenate((P.conj().T,Im),axis=1)

Angles = np.linspace(-np.pi / 2, np.pi / 2, 360)
numAngles = Angles.size
SP = np.empty(numAngles, dtype=complex)

for i in range(numAngles):
    a = np.exp(-1j*2*np.pi*d.reshape(-1,1)*np.sin(Angles[i]))
    SP[i] = (1/(a.conj().T@Q.conj().T@Q@a))[0,0]

SP = np.abs(SP)
SPmax = np.max(SP)
SP = 10 * np.log10(SP / SPmax)
x = Angles * radeg


plt.plot(x, SP)
plt.show()