改进型MUSIC算法（python）

最新推荐文章于 2024-01-05 15:02:54 发布

dongming2017

最新推荐文章于 2024-01-05 15:02:54 发布

阅读量2.3k

点赞数

文章标签： python 算法矩阵

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

版权

改进型MUSIC算法（IMUSIC）用来估计相干信号，通过对阵列输出的协方差矩阵进行处理，转化为Toeplitz矩阵，实现信号的DOA估计。

该算法程序比较简单，核心内容就是一个矩阵转换公式，

$R_{X}=R+I_{v}R^{*}I_{v}$

其中 $R$ 代表原协方差矩阵， $R_{X}$ 代表重构的Toeplitz矩阵， $I_{v}$ 表示反单位矩阵。

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
snr = 10
iwave = theta.size

def coherent_signal_Rxx(d, theta, n, snr):
    A = np.exp(-1j * 2 * np.pi * d.reshape(-1, 1) @ np.sin(theta * derad))
    S = np.random.randn(iwave-1, n)
    S = np.concatenate((S[0,:].reshape(1,-1),S), axis=0)
    X = A @ S
    X = awgn(X, snr)
    Rxx = X @ (X.conj().T) / n
    return Rxx

def IMUSIC(Rxx, iwave, d):
    K = Rxx.shape[0]
    J = np.flip(np.eye(K),axis=0)
    Rxx = Rxx + J@Rxx.conj()@J


    D, EV = LA.eig(Rxx)
    index = np.argsort(D)  # ascending order index
    EN = EV[:, index][:, 0:K - iwave]  #

    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)[0:K] * np.sin(Angles[i]))
        SP[i] = ((a.conj().T @ a) / (a.conj().T @ EN @ EN.conj().T @ a))[0, 0]

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

    return x, SP


co_Rxx = coherent_signal_Rxx(d=d,theta=theta,n=n,snr=snr)
x, SP = IMUSIC(Rxx=co_Rxx,iwave=iwave,d=d)
plt.plot(x, SP)
plt.show()