基于Python的TLS-ESPRIT算法

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([10,30,60]).reshape(1,-1)
n=500
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,10)
Rxx = X@(X.conj().T)/n

def TLS_ESPRIT(K,d,iwave,Rxx):
    dd = d[1]-d[0]
    D,EV = LA.eig(Rxx)
    index = np.argsort(D)
    EN = EV[:,index][:,K-iwave:K]

    Exy = np.concatenate((EN[0:K-1,:],EN[1:K,:]),axis=1)
    E_xy = Exy.conj().T@Exy
    D_xy,EV_xy = LA.eig(E_xy)
    index1 = np.argsort(D_xy)
    EV_xy = EV_xy[:,index1]



    E12 = EV_xy[0:iwave,0:iwave]
    E22 = EV_xy[iwave:2*iwave,0:iwave]
    Psi = -E12@LA.pinv(E22)

    D_si,V_si = LA.eig(Psi)
    ephi = np.angle(D_si)
    angles = -np.arcsin(ephi/(2*np.pi)/dd)*radeg
    return angles

ESPRIT_DOA = []
ESPRIT_DOA = TLS_ESPRIT(K=8,d=d,iwave=iwave,Rxx=Rxx)
for DOA in ESPRIT_DOA:
    print(DOA)

ESPRIT算法建立在子空间旋转不变的基础上，相比于MUSIC算法，不需要全局搜索，减少了运算量。

编写的过程中，易出现的问题：

TLS-ESPRIT算法一共需要进行3次特征分解

1、第一次，对接收信号取K个大特征值对应的特征向量，构成信号子空间。

2、第二次，取K个小特征值对应的特征空间。

一般公式中，矩阵乘法省略两个矩阵之间的乘号，矩阵乘法对应MATLAB中（*），而在Python中矩阵乘法用np.dot()或者@表示，（*）表示对应位相乘（也就是Hadamard积）。