python实现阵列信号处理（三）：多重信号分类Music算法

一、概述

  MUSIC算法是学者 Schmidt 等人 1979 年提出的, 该算法是空间谱估计理论体系中的标志性算法, 它开创了空间谱估计算法研究的新时代, 促进了特征结构类算法的兴起和发展。MUSIC 算法的基本思想是将阵列输出数据的协方差矩阵进行特征分解, 得到与信号分量相对应的信号子空间和与信号分量正交的噪声子空间, 然后利用这 2 个空间的正交性来估计信号的人射方向、极化信息及信号强度等参数。

二、Music算法原理

  在矩阵理论中关于特征值与特征矢量描述如下:
$${\mathbf{R}}_x{\mathbf{u}}_i={\lambda }_i{\mathbf{u}}_i(i=1,2,\cdots ,M)\text{\hspace{1em}(2.1)}$$
  假定特征值按照降序排列:
$${\lambda }_1>{\lambda }_2>{\lambda }_3>\cdots >{\lambda }_M\text{\hspace{1em}(2.2)}$$
  那么：
$${\lambda }_1={\lambda }_{\max }\text{\hspace{1em}(2.3)}$$
  特征矢量组成矩阵:
$$\mathbf{U}=\left[\begin{array}{llll}{\mathbf{u}}_1& {\mathbf{u}}_2& \cdots & {\mathbf{u}}_M\end{array}\right]\text{\hspace{1em}(2.4)}$$
  那么:
$$\begin{array}{r}{\mathbf{R}}_x\mathbf{U}={\mathbf{R}}_x\left[\begin{array}{llll}{u}_1& u_2& \cdots & u_M\end{array}\right]=\left[\begin{array}{llll}{\lambda }_1{u}_1& {\lambda }_2{u}_2& \cdots & {\lambda }_M{u}_M\end{array}\right]=\\ \left[\begin{array}{llll}{u}_1& u_2& \cdots & u_M\end{array}\right]\left[\begin{array}{cccc}{\lambda }_1& 0& \cdots & 0\\ 0& {\lambda }_2& \cdots & 0\\ ⋮& ⋮& & ⋮\\ 0& 0& \cdots & {\lambda }_M\end{array}\right]=\mathbf{U}\end{array}\text{\hspace{1em}(2.5)}$$

  将阵列快拍信号 $\mathbf{x}(t)$ 的协方差矩阵特征矢量矩阵分解为前 $P$ 个特征矢量矩阵 ${\mathbf{U}}_{s(M\times P)}$ 与后 $Q$ 个特征矢量矩阵 ${\mathbf{U}}_{n(M\times Q)}$, 即：

$${\mathbf{U}}_{(M\times M)}=\left[\begin{array}{ll}{\mathbf{U}}_{s(M\times P)}& {\mathbf{U}}_{n(M\times Q)}\end{array}\right]\text{\hspace{1em}(2.6)}$$

  前 $P$ 个大特征值对应特征矢量 ${\mathbf{U}}_s$, 由信号特征矢量组成; 后 $M-P$ 个小特征值对应特征矢量 $U_n$, 由噪声特征矢量组成。

  考虑到实际接收数据矩阵是时间有限长的, 数据协方差矩阵的估计为：

$${\mathbf{R}}_x=\frac{1}{L}\sum _{i=1}^L\mathbf{X}{\mathbf{X}}^{\mathrm{H}}\text{\hspace{1em}(2.7)}$$

  对 ${\mathbf{R}}_x$ 进行特征值分解可以计算得到信号子空间 $U_s$ 噪声子空间 $U_n$ 和由特征值构成的对角矩阵 $\mathbf{\Sigma }$ ，即：

$${\mathbf{a}}_j^{\mathrm{H}}{\mathbf{U}}_n{\mathbf{U}}_n^{\mathrm{H}}{\mathbf{a}}_j=0(j=1,2,\cdots ,P)\text{\hspace{1em}(2.8)}$$

  实际上求 DOA 是以最小优化搜索实现的, 即：

$${\theta }_{\text{Music }}=\underset{\theta }{\arg }\left\{\min \left({\mathbf{a}}^{\mathrm{H}}(\theta ){\mathbf{U}}_n{\mathbf{U}}_n^{\mathrm{H}}\mathbf{a}(\theta )\right)\right\}\text{\hspace{1em}(2.9)}$$

  习惯上, 我们构造类似功率谱函数, 搜索最大值替代搜索最小值, 即：

$$P_{\text{MusIC }}=\frac{1}{{\mathbf{a}}^{\mathrm{H}}(\theta ){\mathbf{U}}_n{\mathbf{U}}_n^{\mathrm{H}}\mathbf{a}(\theta )}\text{\hspace{1em}(2.10)}$$

  按方位角度 $\theta$ 进行搜索,取得峰值的 ${\theta }_1,{\theta }_2,\cdots ,{\theta }_N$ 就是 $N$ 个信 号源的波达方向估计值。

  也可以采用信号子空间方法构造空间谱函数：

$$P_{\mathrm{M}\mathrm{U}\mathrm{S}\mathrm{I}\mathrm{C}}=\frac{1}{{\mathbf{a}}^{\mathrm{H}}(\theta )\left(I-{\mathbf{U}}_S{\mathbf{U}}_{\mathrm{S}}^{\mathrm{H}}\right)\mathbf{a}(\theta )}\text{\hspace{1em}(2.11)}$$

三、python语言实现Music算法

# 导入模块
import numpy as np
from scipy import signal
import matplotlib.pyplot as plt

# 生成快拍数据
def gen_signal(fre, t_0, theta, speed, numbers, space):
    res = []
    for i in range(numbers):
        res.append(np.exp(2j*np.pi*fre*t_0 - 2j*np.pi*fre*i*space*np.cos(theta)/speed))
    return np.array(res)

# 生成方向矢量
def steer_vector(fre, theta, speed, numbers, sapce):
    alphas = []
    for i in range(numbers):
        alphas.append(np.exp(-2j*np.pi*fre*i*space*np.cos(theta)/speed))   
    return np.array(alphas).reshape(-1, 1)
 
# Music算法
def cal_music(fre, speed, numbers, space, signals, method='signal'):
    R_x = np.matmul(signals, np.conjugate(signals.T)) / signals.shape[1]
    lamda, u = np.linalg.eig(R_x)
    u_s = u[:, np.argmax(lamda)].reshape(-1, 1)
    u_n = np.delete(u, np.argmax(lamda), axis=1)
    P = []
    thetas = np.linspace(-np.pi/2, np.pi/2, 180)
    for _theta in thetas:
        _alphas = steer_vector(fre, _theta, speed, numbers, space).reshape(-1, 1)
        if method == 'signal':
            P_x = 1 / np.matmul(np.matmul(np.conjugate(_alphas).T, np.eye(len(u_s)) - np.matmul(u_s, np.conjugate(u_s.T))), _alphas)
        elif method == 'noise':
            P_x = 1 / (np.matmul(np.matmul(np.matmul(np.conjugate(_alphas).T, u_n), np.conjugate(u_n.T)), _alphas))
        else:
            print('there is no ' + method)
            break
        P.append(P_x)
    P = np.array(P).flatten()
    return thetas/np.pi*180, P

# 初始化数据
fs = 20000
# 定义源信号
fre = 200
t = np.arange(0, 0.01, 1/fs)
theta1 = np.pi / 3
theta2 = 2 * np.pi / 3
# 传播速度
speed = 340
# 阵元数量
numbers = 32
# 阵元之间距离
space = 1
# 生成模拟快拍数据
signals = []
for t_0 in t:
    signal1 = gen_signal(fre, t_0, theta1, speed, numbers, space)
    signal2 = gen_signal(fre, t_0, theta2, speed, numbers, space)
    signal = signal1 + signal2
    signals.append(signal.tolist())
signals = np.array(signals)

# Music算法处理结果
thetas, P = cal_music(fre, speed, numbers, space, signals, method='noise')
plt.figure(figsize=(10, 2))
plt.plot(thetas, abs(P))
plt.xlim(-90, 90)
plt.xlabel('degree')
plt.ylabel('mag')
plt.show()

图3.1 Music算法处理结果

四、Tips

  文中如有错误或不清晰的地方，敬请谅解，可留言指出，笔者将及时答复、更改，希望共同进步，谢谢支持。