基于MUSIC算法的信号参量估计

空间谱估计基础算法详解——MUSIC算法

本文叙述在均匀线阵接收信号模型下通过MUSIC算法对多个信号进行DOA（Direction of Arrival，波达方向）估计的推导过程。

全文内容来自于何子述老师的现代数字信号处理及

一 、 信号DOA估计

1.1 信号模型

考虑$N$个远场的窄带信号入射到空间阵列中，其中天线阵列为由$M$个阵元组成的均匀线性阵列(ULA，uniform linear array)，通常情况下可以假设$M>N$。窄带信号$\tilde{s}(t)$可以表示成如下的复数形式：
$$\tilde{s}(t)=s(t)e^{j{\omega }_{0}t}$$
其中$s(t)=u(t)e^{j\psi (t)}$是信号的复包络，$u(t)$和$e^{j\psi (t)}$分别是接收信号的幅度和相位，${\omega }_0$是接收信号的载波角频率。

在均匀阵列中，阵列之间等间距分布，相邻间隔取为$d$，当信号入射角方向为$\theta$时，由于远场信号的入射波到达阵列时可以视为平面波，因此相邻阵元之间的波程差为$dsin\theta$，以线阵0阵元为参考阵元，则不难得出阵元$m$达参考阵元的波程差为$mdsin\theta$。对于平面波，$m$阵元相对于0阵元的传播时延为
$${\tau }_m=\frac{mdsin\theta }{c}=m\tau$$
其中，$\tau =\frac{dsin\theta }{c}$为相邻阵元之间的时延差，$c$是光速。

阵元0的接收信号可以表示为：
$$\tilde{x}(t)=\tilde{s}(t)=s(t)e^{j{\omega }_{0}t}$$
则阵元m的接收信号可以表示为：
$${\tilde{x}}_m(t)=\tilde{x}(t-{\tau }_m)=\tilde{s}(t-{\tau }_m)$$
因为信号$\tilde{s}(t)$是窄带信号，即，$s(t)$是慢变化的，因此有如下简化
$$s(t)\approx s(t-{\tau }_m)，m=0,1,...,M-1$$
所以有
$${\tilde{x}}_m(t)=\tilde{s}(t-{\tau }_m)=s(t-{\tau }_m)e^{j{\omega }_0(t-{\tau }_m)}\approx s(t)e^{j{\omega }_{0}t}{e}^{-j{\omega }_0\tau }=\tilde{s}(t)e^{-j{\omega }_0\tau }$$
定义相邻阵元间的相位差为：
$$\varphi ={\omega }_0\tau =2\pi f_{0}dsin\theta /c=2\pi dsin\theta /\lambda$$
其中$f_0$是入射信号的载波频率，$\lambda$是载波的波长。现在定义列向量：
$$\tilde{x}=[{\tilde{x}}_0(t){\tilde{x}}_1(t)...{\tilde{x}}_{M-1}(t){]}^T$$
以及
$$a(\theta )=[1e^{-j\varphi }...e^{-j(M-1)\varphi }{]}^T$$
则各阵元接收到的信号可以表示成
$$\tilde{x}(t)=a(\theta )\tilde{s}(t)=a(\theta )s(t)e^{j{\omega }_{0}t}$$
其中$a(\theta )$称为导向矢量（steering vector）或者方向矢量，也称为阵列的阵列流形。在阵列信号处理中，复载波$e^{j{w}_{0}t}$不含有信号的有用信息，故通常可以只考虑复基带信号。因此公式（10）对应的离散时间复基带信号可以表示为
$$x(n)=a(\theta )s(n)$$
其中时间变量$n$称为快拍(snapshot)，表示n时刻对所有阵元同时采样。

现在考虑当$K$个信号同时从${\theta }_1，{\theta }_2，...，{\theta }_K$方向入射到阵列中，则天线阵列接收到的离散时间基带信号是个入射信号源的和，即可以表示成
$$\mathbf{x}(n)=\mathbf{a}\left({\theta }_1\right)s_1(n)+\mathbf{a}\left({\theta }_2\right)s_2(n)+\cdots +\mathbf{a}\left({\theta }_K\right)s_K(n)$$
其中第$k$个信号源对应的导向矢量为
$$\mathbf{a}\left({\theta }_k\right)={\left[\begin{array}{llll}1& {\mathrm{e}}^{-\mathrm{j}{\varphi }_k}& \cdots & {\mathrm{e}}^{-\mathrm{j}(M-1){\varphi }_k}\end{array}\right]}^{\mathrm{T}},{\varphi }_k=2\pi d\sin {\theta }_k/\lambda ,k=1,2,\cdots ,K$$
分别定义信号向量$s(n)$和方向矩阵(direction matrix)$A$为
$$\begin{array}{l}\mathbf{s}(n)={\left[\begin{array}{cccc}{s}_1(n)& s_2(n)& \cdots & s_K(n)\end{array}\right]}^{\mathrm{T}}\in {\mathbb{C}}^{K\times 1}\\ \mathbf{A}=\left[\begin{array}{llll}\mathbf{a}\left({\theta }_1\right)& \mathbf{a}\left({\theta }_2\right)& \cdots & \mathbf{a}\left({\theta }_K\right)\end{array}\right]\\ =\left[\begin{array}{cccc}1& 1& \cdots & 1\\ {\mathrm{e}}^{-\mathrm{j}{\varphi }_1}& {\mathrm{e}}^{-\mathrm{j}{\varphi }_2}& \cdots & {\mathrm{e}}^{-\mathrm{j}{\varphi }_K}\\ ⋮& ⋮& \ddots & ⋮\\ {\mathrm{e}}^{-\mathrm{j}(M-1){\varphi }_1}& {\mathrm{e}}^{-\mathrm{j}(M-1){\varphi }_2}& \cdots & {\mathrm{e}}^{-\mathrm{j}(M-1){\varphi }_K}\end{array}\right]\in {\mathbb{C}}^{M\times K}\end{array}$$
则公式(14)可以写成向量形式，即：
$$x(n)=As(n)\in {\mathbb{C}}^{M\times 1}$$
从上式中可以看出，阵列接收信号可以由信号矢量和方向矩阵完全确定，由于接收机不可避免地会引入加性噪声，因此在实际工程中阵列接收信号可以表示成
$$x(n)=As(n)+v(n)$$
其中，$v(n)\in {\mathbb{C}}^{M\times 1}$为噪声向量。公式(10)的展开形式为
$$\left[\begin{array}{c}{x}_0(n)\\ x_1(n)\\ ⋮\\ x_{M-1}(n)\end{array}\right]=\left[\begin{array}{cccc}1& 1& \cdots & 1\\ {\mathrm{e}}^{-j{\varphi }_1}& {\mathrm{e}}^{-\mathrm{j}{\varphi }_2}& \cdots & {\mathrm{e}}^{-\mathrm{j}{\varphi }_K}\\ ⋮& ⋮& \ddots & ⋮\\ {\mathrm{e}}^{-\mathrm{j}(M-1){\varphi }_1}& {\mathrm{e}}^{-\mathrm{j}(M-1){\varphi }_2}& \cdots & {\mathrm{e}}^{-\mathrm{j}(M-1){\varphi }_K}\end{array}\right]\left[\begin{array}{c}{s}_1(n)\\ s_2(n)\\ ⋮\\ s_K(n)\end{array}\right]+\left[\begin{array}{c}{v}_0(n)\\ v_1(n)\\ ⋮\\ v_{M-1}(n)\end{array}\right]$$
通常情况下，可以假设各个阵元接收到的噪声是均值为0、方差为${\sigma }^2$的高斯白噪声，并且不同阵元间的接收噪声相互独立，信号与噪声之间也噪声独立，即$\forall n,l$都有
E { v i ( n ) } = 0 , E { v i ( n ) v i ∗ ( l ) } = σ 2 δ ( n − l ) , E { v i ( n ) v m ∗ ( l ) } = 0 , i ≠ m E { s i ( n ) v m ∗ ( l ) } = 0 , \begin{array}{l} \mathrm{E}\left\{v_{i}(n)\right\}=0, \\ \mathrm{E}\left\{v_{i}(n) v_{i}^{*}(l)\right\}=\sigma^{2} \delta(n-l), \\ \mathrm{E}\left\{v_{i}(n) v_{m}^{*}(l)\right\}=0, \quad i \neq m \\ \mathrm{E}\left\{s_{i}(n) v_{m}^{*}(l)\right\}=0, \end{array} E{vi​(n)}=0,E{vi​(n)vi∗​(l)}=σ2δ(n−l),E{vi​(n)vm∗​(l)}=0,i=mE{si​(n)vm∗​(l)}=0,​
或者等价的向量形式
E { v ( n ) } = 0 E { v ( n ) v H ( l ) } = σ 2 I δ ( n − l ) E { s ( n ) v H ( l ) } = 0 \begin{array}{l} \mathrm{E}\{\boldsymbol{v}(n)\}=\mathbf{0} \\ \mathrm{E}\left\{\boldsymbol{v}(n) \boldsymbol{v}^{\mathrm{H}}(l)\right\}=\sigma^{2} \boldsymbol{I} \delta(n-l) \\ \mathrm{E}\left\{\boldsymbol{s}(n) \boldsymbol{v}^{\mathrm{H}}(l)\right\}=\mathbf{0} \end{array} E{v(n)}=0E{v(n)vH(l)}=σ2Iδ(n−l)E{s(n)vH(l)}=0​

1.2 基于MUSIC算法的多个信号的DOA估计

公式实在是太难敲了，直接把书上的截图放上来了

二、信号频率估计

先假设信号$x(n)$是复正弦信号加上高斯白噪声，则可以定义信号模型为
$$x(n)=\sum _{k=1}^K{\alpha }_k{\mathrm{e}}^{\mathrm{j}{\omega }_{k}n}+v(n)$$
其中，${\alpha }_k=|{\alpha }_k|e^{j{\phi }_k}$和${\omega }_k$分别表示信号的复幅度和角频率。初始相位${\phi }_k$是在$[0,2\pi ]$范围内均匀分布的随机变量，同时当$i\ne j$时，${\phi }_i$和${\phi }_k$相互独立；$v(n)$是零均值、方差为${\sigma }^2$的白噪声，同时与信号相互独立。

其向量形式可以写成：
$$x(n)=As(n)+v(n)\in {\mathbb{C}}^{M\times 1}$$
此时方向矩阵$A$、$a(\omega )$、$s(n)$、$v(n)$分别定义为
$$\mathbf{A}=\left[\begin{array}{llll}\mathbf{a}\left({\omega }_1\right)& \mathbf{a}\left({\omega }_2\right)& \cdots & \mathbf{a}\left({\omega }_K\right)\end{array}\right]=\left[\begin{array}{cccc}1& 1& \cdots & 1\\ {\mathrm{e}}^{-j{\omega }_1}& {\mathrm{e}}^{-j{\omega }_2}& \cdots & {\mathrm{e}}^{-j{\omega }_K}\\ ⋮& ⋮& \ddots & ⋮\\ {\mathrm{e}}^{-j(M-1){\omega }_1}& {\mathrm{e}}^{-\mathrm{j}(M-1){\omega }_2}& \cdots & {\mathrm{e}}^{\mathrm{j}(M-1){\omega }_K}\end{array}\right]\in {\mathbb{C}}^{M\times K}$$

$$\mathbf{a}(\omega )=\left[\begin{array}{c}1\\ {\mathrm{e}}^{-\mathrm{j}\omega }\\ ⋮\\ {\mathrm{e}}^{-\mathrm{j}(M-1)\omega }\end{array}\right],\mathbf{s}(n)=\left[\begin{array}{c}{\alpha }_1{\mathrm{e}}^{\mathrm{j}{\omega }_{1}n}\\ {\alpha }_2{\mathrm{e}}^{\mathrm{j}{\omega }_{2}n}\\ ⋮\\ {\alpha }_K{\mathrm{e}}^{\mathrm{j}{\omega }_{K^n}}\end{array}\right],\mathbf{v}(n)=\left[\begin{array}{c}v(n)\\ v(n-1)\\ ⋮\\ v(n-M+1)\end{array}\right]$$

2.2 基于MUSIC算法的信号频率估计