非周期序列设计（CAN算法）

本文主要内容来自于何浩的著作《有源感知系统波形设计算法》 建议阅读英文版，中文版有些符号有点问题。

He, H., Li, J., & Stoica, P. (2012). Waveform Design for Active Sensing Systems: A Computational Approach. Cambridge: Cambridge University Press. doi:10.1017/CBO9781139095174

一、设计准则

在雷达系统中，模糊函数是分析给定波形分辨率、副瓣性能以及多普勒和距离模型方面十分有用的工具。自相关函数可以看作零多普勒的模糊函数剖面。考虑发射信号$s(t)$由$N$ 个符号组成：
$$s(t)=\sum _{n=1}^{N}x(n)p_n(t)\text{\hspace{1em}(1)}$$
其中，$p_n(t)$为成形脉冲， { x ( n ) } n = 1 N \{x(n)\}_{n=1}^N {x(n)}n=1N​表示其中的$N$个符号。当考虑恒模约束，可以将符号函数表示为：
$$x(n)=e^{j\varphi (n)},n=1,\dots ,N\text{\hspace{1em}(2)}$$
进一步，$x(n)$的自相关函数可以表示为：
$$r(k)=\sum _{n=k+1}^{N}x(n)x^{\ast }(n-k)=r^{\ast }(-k),k=0,\dots ,N-1\text{. }\text{\hspace{1em}(3)}$$
上述自相关函数 { r ( k ) } k = − N + 1 N − 1 \{r(k)\}_{k=-N+1}^{N-1} {r(k)}k=−N+1N−1​更准确的叫法是非周期自相关函数（Aperiodic Auto Correlation Function, AACF），并且$r(0)$称为同相相关，总是等于信号能量。

如果仅考虑单个脉冲序列的性能，可以通过对序列的自相关函数进行设计，进而得到所需的脉冲序列。一般而言，非周期序列设计的目标是使 { r ( k ) } k ≠ 0 \{r(k)\}_{k\ne 0} {r(k)}k=0​尽可能小。传统准则主要有两个：积分旁瓣电平（Integrated Sidelobe Level, ISL）和峰值旁瓣电平 （Peak Sidelobe Level, PSL）。ISL 衡量的是波形 ACF 整体旁瓣电平，以 ISL 为准则优化得到的波形旁瓣整体性能较均衡；而 PSL 衡量的是波形 ACF 旁瓣峰值，码长越长，以 PSL 为准则优化得到的旁瓣峰值越低。

二、设计方法

Stoica P, He H, Li J. New algorithms for designing unimodular sequences with good correlation properties [J]. IEEE Transactions on Signal Processing, 2009, 57(4): 1415-1425.

CAN (Cyclic algorithm-new) 算法 采用积分旁瓣电平作为序列性能的衡量准则，其定义如下：
$$\mathrm{ISL}=\sum _{\begin{array}{c}k=-(N-1)\\ k\ne 0\end{array}}^{N-1}|r(k){|}^2=2\sum _{k=1}^{N-1}|r(k){|}^2\text{\hspace{1em}(4)}$$
另一种等价表述是品质因数（Merit Factor, MF）：
$$\mathrm{MF}=\frac{|r(0){|}^2}{{\sum }_{\begin{array}{c}k=-(N-1)\\ k\ne 0\end{array}}^{N-1}|r(k){|}^2}=\frac{N^2}{\mathrm{ISL}}\text{\hspace{1em}(5)}$$
在恒模约束的条件下，CAN是一种可以得到使ISL或与ISL相关的准则最小的恒模序列的高效优化方法。其算法流程如下所示。

2.1 CAN算法推导

CAN算法第一步是将ISL准则转换到频域表示。 对于任何一个频点$\omega \in \left[0,2\pi \right]$，
$${\left|\sum _{n=1}^{N}x(n)e^{-j\omega n}\right|}^2=\sum _{k=-(N-1)}^{N-1}r(k)e^{-j\omega k}\triangleq \Phi (\omega )\text{\hspace{1em}(6)}$$
此时，利用帕塞瓦尔定理，可以将ISL度量等价为：
$$\mathrm{ISL}=\frac{1}{2N}\sum _{p=1}^{2N}{\left[\Phi \left({\omega }_p\right)-N\right]}^2=\frac{1}{2N}\sum _{p=1}^{2N}{\left({\left|\sum _{n=1}^{N}x(n)e^{-j{\omega }_{p}n}\right|}^2-N\right)}^2\text{\hspace{1em}(7)}$$
其中， { ω p } \{\omega_p\} {ωp​}是以下傅里叶频率：
$${\omega }_p=\frac{2\pi }{2N}p,p=1,\dots ,2N\text{\hspace{1em}(8)}$$
这个结果有一个明显的直观解释：最小化ISL使序列表现得像白噪声，因此其周期图在频率上应该几乎是恒定的。

一个问题是，式(6)中的准则是 { x ( n ) } \left\{x(n)\right\} {x(n)}的四次函数。但是，Petre Stoica发现，对 { x ( n ) } \left\{x(n)\right\} {x(n)}的最小化几乎等同于以下 { x ( n ) } \left\{x(n)\right\} {x(n)}的二次函数：
min ⁡ { x ( n ) } n = 1 N ; { ψ p } p = 1 2 N ∑ p = 1 2 N ∣ ∑ n = 1 N x ( n ) e − j ω p n − N e j ψ p ∣ 2 . (9) \min _{\{x(n)\}_{n=1}^{N} ;\left\{\psi_{p}\right\}_{p=1}^{2 N}} \sum_{p=1}^{2 N}\left|\sum_{n=1}^{N} x(n) e^{-j \omega_{p} n}-\sqrt{N} e^{j \psi_{p}}\right|^{2} .\tag{9} {x(n)}n=1N​;{ψp​}p=12N​min​p=1∑2N​ ​n=1∑N​x(n)e−jωp​n−N ​ejψp​ ​2.(9)
（上式等价性的说明可以参考《有源感知系统波形设计算法》的附录2A）

当完成频域变换后，令
$${\mathbf{a}}_p^H=\left[e^{-j{\omega }_p}\cdots e^{-j2N{\omega }_p}\right],\text{\hspace{1em}(10)}$$,
设${\mathbf{A}}^H$是一个$2N\times 2N$的酉矩阵，由
$${\mathbf{A}}^H=\frac{1}{\sqrt{2N}}\left[\begin{array}{c}{\mathbf{a}}_1^H\\ ⋮\\ {\mathbf{a}}_{2N}^H\end{array}\right]\text{\hspace{1em}(11)}$$
设$\mathbf{z}$是序列 { x ( n ) } n = 1 N \left\{x(n)\right\}_{n=1}^N {x(n)}n=1N​进行补零后的序列，填充了$N$个$0$：
$$\mathbf{z}={\left[x(1),\cdots ,x(N),0,\cdots ,0\right]}_{2N\times 1}^T\text{\hspace{1em}(12)}$$
将（6）中的准则可以重写成更紧凑的形式
$$\underset{\mathbf{z}}{\min }\Vert {\mathbf{A}}^H\mathbf{z}-\mathbf{v}{\Vert }^2,\text{\hspace{1em}(13)}$$
其中，$\mathbf{v}$为相位辅助变量，
$$\mathbf{v}=\frac{1}{\sqrt{2}}{\left[e^{j{\psi }_1},\cdots ,e^{j{\psi }_{2N}}\right]}^T.\text{\hspace{1em}(14)}$$
对于给定的$\mathbf{z}$，（13）关于 { ψ p } \left\{\psi_p\right\} {ψp​}的最小值是直接计算的，即表示$\mathbf{z}$的傅里叶变换：
$$\mathbf{f}={\mathbf{A}}^H\mathbf{z}\text{\hspace{1em}(15)}$$
$${\psi }_p=\arg \left(f_p\right),p=1,\cdots ,2N\text{\hspace{1em}(16)}$$
同样地，对于一个给定的$\mathbf{v}$，利用$\mathbf{A}$的正交性，直接求解$\mathbf{v}$的逆傅里叶变换（IFFT）：
$$\mathbf{g}=\mathbf{Av}\text{\hspace{1em}(17)}$$

因为$\Vert {\mathbf{A}}^H\mathbf{z}-\mathbf{v}{\Vert }^2=\Vert \mathbf{z}-\mathbf{Av}{\Vert }^2$，所以最小化序列 { x ( n ) } \left\{x(n)\right\} {x(n)}是由
$$x(n)=e^{j\arg (g_n)},n=1,\cdots ,N.\text{\hspace{1em}(18)}$$

下表总结了循环局部最小化的CAN算法。由于只涉及简单的FFT和IFFT操作，即使在PC上也可以计算码长达$10^6$的序列。