基于python的ROOT-MUSIC算法

       相比于ESPRIT和MUSIC算法的改写，ROOT-MUSIC花费时间最长，大概一天半的样子（又是一度想放弃）。

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

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 RMUSIC(d,iwave,Rxx):
    K = d.size
    D, EV = LA.eig(Rxx)
    index = np.argsort(D)  # ascending order index
    EN = EV[:, index][:, 0:K - iwave]  #

    z = sympy.symbols('z')
    k1 = np.arange(0,K).reshape(1,-1)
    pz = z**k1
    pz1 = (z**(-1))**k1
    fz = (pz1@EN@(EN.conj().T)@(pz.T)).item()

    fz = sympy.simplify(z**(K-1)*fz)

    zx = sympy.roots(fz,z)
    zx0 = np.array(tuple(zx.keys()))
    zx1 = np.abs(np.abs(zx0)-1)
    index = np.argsort(zx1)
    rx_index = np.arange(0,2*iwave,2)
    rx = zx0[index][rx_index]
    rx1 = np.fromiter(rx, dtype=complex)
    DOAs = np.arcsin(-np.angle(rx1)/np.pi)*radeg
    return DOAs

RMUSCI_DOA = []
RMUSCI_DOA = RMUSIC(d=d,iwave=iwave,Rxx=Rxx)
for DOA in RMUSCI_DOA:
    print(DOA)

因为其中涉及到多项式表达和求解，只采用sympy和numpy中的一种很难完成（当然也可能是我对两种库了解不够），而两种库的数据输出格式并不相同，经常要转换，可以看到代码中有些地方为了匹配格式，显得很生硬。这里sympy主要负责多项式的表达和求解，numpy负责排序，求角度等其他计算。并且该代码执行计算时间较长。