0. 参考

1. 理论基础

1.1 信号

若有正弦信号

$y = \sin{\omega t} =\sin{2\pi ft}

则有

$T = \frac{2\pi}{\omega} = \frac{1}{f}$

1.2 共轭转置

若$\emph{X} = [x{ij}] = [a{ij}+b{ij}\cdot j]$，则其共轭转置$\emph{X}^{H} = [a{ji}-b_{ji} \cdot j]$。

1.3 欧拉公式

对于任意复数

$z = a + b \cdot j = r(\cos\theta + j\cdot \sin\theta) = re^{j\cdot\theta}$

其中

$r = \sqrt{a^2+b^2},\qquad\theta = \arctan\left( \frac{b}{a} \right)$

1.4 拉格朗日乘子法

$\underset{\overset{\rightarrow}{x}}{\min \text{or} \max}\qquad f(\overset{\rightarrow}{x})$

$\text{subject to} \qquad g(\overset{\rightarrow}{x}) = 0$

其中

$\overset{\rightarrow}{x} = [x_1, x_2, ... , x_n]$

令

$\left{ \begin{aligned} &\nabla f(\overset{\rightarrow}{x}) + \lambda \nabla g(\overset{\rightarrow}{x}) &= 0 \ &g(\overset{\rightarrow}{x}) &= 0 \end{aligned} \right.$

解得

$\overset{\rightarrow}{x} = \overset{\rightarrow}{x^*}$

即为极值点，但需要检查其极值性。

1.5 Wirtinger梯度

在复数域里，复函数$f(z)$对$z$严格求导需要满足Cauchy–Riemann条件（即$f(z)$为解析函数）。

但在信号处理里常用到的功率、相关等函数大多不是解析函数，所以无法使用标准复分析。

Wirtinger梯度的解决方法为，给定复数$z = x + jy$，把$ z $和它的共轭$z^∗ $当作两个独立变量，并定义两个偏导：

$\begin{aligned} &\frac{\partial f(z, z^)}{\partial z} = \frac{1}{2}\left(\frac{\partial}{\partial x} - j \frac{\partial}{\partial y}\right) \ &\frac{\partial f(z, z^)}{\partial z^*} = \frac{1}{2}\left(\frac{\partial}{\partial x} + j \frac{\partial}{\partial y}\right) \end{aligned}$

则有：

$\begin{aligned} & \frac{\partial z}{\partial z} = 1, \quad \frac{\partial z^}{\partial z} = 0 \ & \frac{\partial z^}{\partial z} = 0, \quad \frac{\partial z^}{\partial z^} = 1 \ & \frac{\partial (zz^)}{\partial z} = z^, \quad \frac{\partial (zz^)}{\partial z^} = z^* \end{aligned}$

并有：

$\frac{\partial f}{\partial z} = \overline{\frac{\partial f}{\partial z^*}}$

需要注意的是，在对一个复数向量求Wirtinger梯度时，梯度的维度应当和该向量的维度相同，具体可能体现在结果需要做个转置等。

1.6 奈奎斯特定理

要从抽样信号中无失真地恢复原信号，抽样频率应大于2倍信号最高频率。

抽样频率小于2倍频谱最高频率时，信号的频谱有混叠（相关性？）。

1.7 傅里叶变换

傅里叶变换是傅里叶级数在非周期函数上的推广，是时域$t$和频域$\omega$的转换。

1.7.1 连续傅里叶级数

傅里叶级数和变换学习笔记

对于周期为$T$（数字频率$w_0=\frac{2\pi}{T}$）的信号$f(t)$，那么它可以展开为一串正弦和余弦信号的叠加，且正余弦信号的数字频率间隔$\Delta w = w_0=\frac{2\pi}{T}$：

$\begin{aligned} f(t)&=\frac{a0}{2}+\sum{n=1}^{\infty}[a_n\cos(n\omega_0 t)+b_n\sin(n\omega0 t)] \ &= \sum{n=-\infty}^{\infty}c_ne^{jn\omega_0 t} \end{aligned} $

其中：

$\begin{aligned} an &= \frac{2}{T}\int{t_0}^{t_0+T}f(t)\cos(n\omega_0 t)dt \ bn &= \frac{2}{T}\int{t_0}^{t_0+T}f(t)\sin(n\omega_0 t)dt \ cn &= \frac{1}{T}\int{t_0}^{t_0+T}f(t)e^{-jn\omega_0 t}dt \end{aligned}$

1.7.2 连续傅里叶变换FT

在上述周期为$T$的信号$f(t)$，有：

$\begin{aligned} f(t)&= \sum_{n=-\infty}^{\infty}c_ne^{jn\omega0 t} \ &= \sum{n=-\infty}^{\infty}[\frac{1}{T}\int_{t_0}^{t_0+T}f(\tau)e^{-jn\omega_0 \tau}d\tau]e^{jn\omega0 t} \ &= \frac{1}{2\pi}\sum{n=-\infty}^{\infty}[\int_{t_0}^{t_0+T}f(\tau)e^{-jn\omega_0 \tau}d\tau]e^{jn\omega_0 t}\cdot\Delta \omega \end{aligned} $

当$f(t)$不再是周期信号，也就是$T\to \infty$时，有$\Delta w = w_0 \to 0$，

令基信号的数字频率为$\omega$，则$\omega = n\omega_0 \to 0$，

此时$\omega$由离散变量变为连续变量，上述$f(t)$可变形为：

$ f(t) = \frac{1}{2\pi}\int{-\infty}^{\infty}[\int{-\infty}^{\infty}f(\tau)e^{-j\omega \tau}d\tau]e^{j\omega t}\cdot d\omega$

令：

$F(\red{\omega}) = \mathcal{F}[f(\tau)] = \int_{-\infty}^{\infty}f(\tau)e^{-j\red{\omega} \tau}d\tau$

则有：

$f(t)=\mathcal{F}^{-1}[F(\omega)] = \frac{1}{2\pi}\int_{-\infty}^{\infty}F(\omega)e^{j\omega t}d\omega$

或令$f = \frac{\omega}{2\pi} $：

$X(\red{f}) = \mathcal{F}[x(t)] = \int_{-\infty}^{\infty}x(t)e^{-j2\pi \red{f} t}dt$

$x(t)=\mathcal{F}^{-1}[X(f)] = \frac{1}{2\pi}\int_{-\infty}^{\infty}X(f)e^{j2\pi ft}df$

1.7.3 离散傅里叶级数

对于周期为$N$的离散信号$x[n]$：

$x[n] = \sum_{k=0}^{N-1}c_ke^{jk\frac{2\pi}{N}n}$

其中：

$ck = \frac{1}{N}\sum{n=0}^{N-1}x[n]e^{-jk\frac{2\pi}{N}n}$

1.7.4 离散傅立叶变换DFT

由上述离散信号$x[n]$可得：

$x[n] =\frac{1}{N} \sum{k=0}^{N-1} [\sum{\tau=0}^{N-1}x[\tau]e^{-jk\frac{2\pi}{N}\tau}] e^{jk\frac{2\pi}{N}n}$

令：

$X[\red{k}] = \sum_{\tau=0}^{N-1}x[\tau]\cdot e^{-j\red{k}\frac{2\pi}{N}\tau}$

则有：

$x[n] = \frac{1}{N}\sum_{k=0}^{N-1}X[k]\cdot e^{jk\frac{2\pi}{N} n}$

1.7.5 快速傅里叶变换FFT

TODO

2. 名词解释

2.1 窄带信号

指信号的带宽小于信号载频的$10\%$，其数字频率$w_0$一般很大，非窄带信号一般都称为宽带信号。

正余弦信号是典型的窄带信号，通信雷达领域常用的chirp信号是宽带信号。

2.2 波达方向（Direction of Arrival）

根据各个传感器接收到的信号，估计声源/信号是从哪个方向来的。

2.3 均匀线阵（Uniform Linear Array）

一种阵元（传感器）的排列方式：阵元在一条直线上排列，且间隔相同。

2.4 导向矢量（Steering Vector）

用于描述在某种阵元排列方式下，不同阵元接收信号的相位差关系。

3. DOA算法

3.0 前提

假设有$M$个阵元以$ULA$排列，沿$x$轴，阵元间距为$d$，第$m$个阵元位置为$x_m=md$（以$m=0$为参考阵元），平面波的入射角为$\theta$（与阵列的法线夹角）。

则相邻阵元的路径差为：$\Delta d = d\sin \theta$，

因此第$m$个阵元相对参考的时延为：$\tau_m = \frac{md\sin \theta}{c}$，

$c$为波速（电磁波为光速约$3*10^8 m/s$，声波为声速约$340m/s$）。

若信号为窄带单频分量$s(t)=Ae^{j\omega_0 t}$，第$n$个阵元接收到的实际信号（相对参考阵元）为：

$s(t-\tau_m) = Ae^{j\omega_0(t-\tau_m)} = s(t)\cdot e^{-j\omega_0 \tau_m}$

令波数$k=\frac{2\pi}{\lambda}$，

得到第$n$阵元上的相位因子为：$e^{j\omega_0 \tau_m} = e^{j\omega_0 \frac{md\sin \theta}{c}} = e^{j\frac{2\pi}{\lambda}md\sin\theta} = e^{jkmd\sin\theta}$（$\omega_0 = 2\pi f$, $c = \lambda f$）。

因此相邻阵元的相位差为：$\Delta \phi = \omega_0 \frac{d\sin\theta}{c} = \frac{2\pi}{\lambda} d\sin\theta = kd\sin\theta$。

由此可得导向矢量：

$ a(\theta) = \begin{bmatrix} 1 \ e^{-jkd\sin\theta} \ e^{-jk2d\sin\theta} \ \vdots \ e^{-jk(M-1)d\sin\theta} \end{bmatrix} $

若有$K$个同频窄带信号（不考虑噪声），则阵元体在时刻$t$的接受信号为：

$X(t){M×1} = A{M×K}S_{K×1}(t), \qquad A = [\alpha(\theta_1), ..., \alpha(\theta_k)]$

3.1 朴素算法

假设单个窄带信号的入射角为$\alpha$，那么该信号对于给定$ULA$阵列的导向矢量为：

$ a(\alpha) = \begin{bmatrix} 1 \ e^{-jkd\sin\alpha} \ e^{-jk2d\sin\alpha} \ \vdots \ e^{-jk(M-1)d\sin\alpha} \end{bmatrix} $

用该导向矢量$a(\alpha)$和阵列实际接收到的信号$X(n)$做内积，得：

$\begin{aligned} y &= a^H{1×M}(\alpha)\cdot X{M×1}(n) = a^H(\alpha) a(\theta) s(n) \ &= \left[1 + e^{-jkd(\sin\theta-\sin\alpha)} + e^{-jk2d(\sin\theta-\sin\alpha)} + ... + e^{-jk(M-1)d(\sin\theta-\sin\alpha)}\right]s(n) \ &\le M_{max}s(n) \end{aligned}$

等号仅在$\alpha = \theta$处取得。

由此很自然的得到一$DOA$估计的朴素算法：

$\begin{aligned} &for: \ \ \alpha \ from\ -\frac{\pi}{2}\ to \ \frac{\pi}{2} \ &\qquad calculate \ \ y\alpha = a^H(\alpha)\cdot X(n) \ &end \ \ for \ &\theta = \underset{\alpha}{\max} \ y\alpha \end{aligned}$

此外，假设有多个同频窄带信号，真实入射角分别为$\theta_1, \theta_2, ..., \theta_k$，

此时有（假设多个信号的总接收值等于单个信号接收值累加）：

$X{M×1}(n) =\begin{bmatrix} \sum{k'=1}^{K}e^{-jkd\sin k'} \ \sum{k'=1}^{K}e^{-jk2d\sin k'} \ \vdots \ \sum{k'=1}^{K}e^{-jk(M-1)d\sin k'} \end{bmatrix} s(n)
$

$\red{ \begin{aligned} y &= a^H{1×M}(\alpha)\cdot X{M×1}(n) \ &=

\left[ 1

e^{-jkd\sin\alpha}\sum_{k'=1}^{K}e^{-jkd\sin k'}
e^{-jk2d\sin\alpha}\sum_{k'=1}^{K}e^{-jk2d\sin k'}
...
e^{-jk(M-1)d\sin\alpha}\sum_{k'=1}^{K}e^{-jk(M-1)d\sin k'} \right]

s(n) \ &\le M_{max}s(n) \end{aligned} }$

该朴素算法也具有一定作用，但当发射信号的入射角比较接近时，该算法的识别分辨率较低，只能识别出一个入射角。

3.2 MVDR

假设只有一个远场窄带信号从$\theta$方向入射，

对各个阵元接收到的信号$x(n)$进行加权求和，权系数$w_{M\times 1} = \left[ w_0, w1, w{M-1} \right]^T$，

输出$y(n) = w^H x(n) = w^Ha(\theta)s(n)$，

该阵列输出信号的平均功率为：

$P(n) = E\left[ \left| y(n) \right|^2 \right] = E\left[ w^Hx(n)x^H(n)w \right] = w^HRw$

其中$R = E\left[ x(n)x^H(n) \right]$为$x(n)$的协方差矩阵，并且有$R^H=R$，得$R$和$R^{-1}$均为Hermitian矩阵。

猜想$\theta = \theta_0$，则阵列接收信号$x_0(n) = a(\theta_0)s(n)$，为使$x_0(n)$通过空域滤波器后无失真，

权矢量需满足：$w^Ha(\theta_0) = 1$。

如果同时尽可能让$P(n)$的值越小，则该权矢量$w$实现了对干扰信号及噪声抑制和期望信号估计的功能。

构造拉格朗日代价函数如下：

$J(w,\lambda)=w^HRw+\lambda(w^Ha(\theta_0)-1)$

计算Wirtinger梯度：

$\nabla J(w^*) = Rw-\lambda a(\theta_0) = 0$

得：

$w = \lambda R^{-1}a(\theta_0)$

又：

$w^Ha(\theta_0)-1=0$

所以：

$\begin{aligned} &\left( \lambda R^{-1}a(\theta_0) \right)^H a(\theta_0) &= 1 \ &\lambda (a(\theta_0))^HR^{-1}a(\theta_0)&=1 \end{aligned}$

即：

$\lambda = \frac{1}{(a(\theta_0))^HR^{-1}a(\theta_0)}$

最后再带回$w$：

$\begin{aligned} w &= \frac{R^{-1}a(\theta_0)}{(a(\theta_0))^HR^{-1}a(\theta0)} \ P{\theta_0} &= \frac{1}{(a(\theta_0))^HR^{-1}a(\theta_0)} = \lambda \end{aligned}$

其中$(a(\theta_0))^HR^{-1}a(\theta_0)$为标量。

4. 代码

4.1 MVDR

4.1.1 计算导向矢量

$ a(\theta) = \begin{bmatrix} 1 \ e^{-j\frac{2\pi}{\lambda}d\sin\theta} \ e^{-j\frac{2\pi}{\lambda}2d\sin\theta} \ \vdots \ e^{-j\frac{2\pi}{\lambda}(M-1)d\sin\theta} \end{bmatrix} $

def steering_vector(M, d, wavelength, theta_deg):
    """
    :param M: ULA 阵元数
    :param d: 阵元间距: (标量)
    :param wavelength: 波长: (标量)
    :param theta_deg: 角度
    :return: steering vector: M x 1
    """
    theta = np.deg2rad(theta_deg)
    m = np.arange(M)  # [0, 1, ..., M-1]
    v = np.exp(-1j * 2*np.pi / wavelength * m * d * np.sin(theta))

    return v.reshape((M,1))

4.1.2 计算$\theta_0$的功率$P$