问题引入:自相关(auto-correlation)

问题引入:自相关(auto-correlation)

在感知系统中(比如雷达声纳),常会遇到这样的问题:当我们发射一串信号后,我们需要从反射的信号中分析出我们感兴趣的内容。比如可以通过检测返回声波信号的多普勒频移来判断目标的速度。在感知系统中,有两个因素对系统的准确性起关键作用,分别是接受滤波器传输波形。常用的接受滤波器有匹配滤波器,匹配滤波器能迅速从杂乱的信号中提取我们感兴趣的波形。常用的传输波形有Golomb波形,它的自相关性比较低,这样的波形不容易被混淆。本文主要研究第二个因素,设计一个低自相关性的传输波形。

  • 自相关

    有两种自相关定义,分别是周期自相关非周期自相关,周期自相关相当于对信号做了周期延拓,而非周期自相关相当于对信号做了补零(padding)。非周期自相关定义如下:
    image-20200613204309780
    其中,x^*(n)代表x(n)的共轭,为了区别周期自相关,用\tilde r(k)表示周期自相关,其定义如下:
    image-20200613204324742
    可以发现,两者都是长度为N的序列。若我们将周期自相关的定义式拆开,可以得到非周期自相关与周期自相关的关系。
    image-20200613204342354

  • Intergrated Sidelobe Level(ISL)

    为了更加直接地表示自相关性,我们引入ISL,其定义如下:
    {\rm ISL} = \sum_{k=-(N-1)\\k\neq0}^{N-1}|r(k)|^2=2\sum_{k=1}^{N-1}|r(k)|^2
    还有另一种与ISL相似的参数,merit factor(MF),其定义如下:
    {\rm MF} = \frac{|r(0)|^2}{2\sum_{k=1}^{N-1}|r(k)|^2}=\frac{N^2}{\rm ISL}

    我们的目标是设计一个高效的算法,使得一个幺模序列(即|x(n)|=1)的ISL最小,或是MF最大。

常见的低自相关性波形

  • Chirp信号(线性调频信号)

    Chirp是一种常用的雷达信号,其作用距离远,分辨率高。它的连续表达式为:
    s(t)=\frac{1}{\sqrt{T}}e^{j\pi \frac{B}{T}t}

    T=100s,B=1{\rm Hz}s(t)的实部如下:

image-20200611205729021

​ 它的自相关函数r (\tau)如下,其旁瓣峰值为-13.4dB。

image-20200611205741393

  • Golomb序列

    我们对s(t)进行采样,取采样间隔为t_s=\frac{1}{B},则总采样数N=BT,据此,我们可以得出以下序列:
    x(n)\triangleq s(nt_s)=e^{j\pi \frac{B}{T}(\frac{n}{B})^2}=e^{j\pi \frac{n^2}{N}},n=1,…,N
    上述定义只有在n为偶数时成立。当n为奇数时,其定义如下:
    x(n)\triangleq s(nt_s)=e^{j\pi \frac{n(n-1)}{N}},n=1,…,N
    N=100时,它的自相关函数r (\tau)如下,其旁瓣峰值为-16.3dB。

    image-20200611211630488

  • 其他

    还有一些与Golomb序列相似的序列,比如Frank序列、Zadoff-Chu序列、P4序列等,这里不再一一介绍。这些序列的周期自相关性都有很好的表现,但非周期自相关的序列设计鲜有人设计,以下我们将主要对非周期自相关进行讨论。

Cyclic Algorithm-New(CAN)

  • 目标函数

    在频域中进行相关操作能减少时间复杂度,因此首先将序列x(n)变换到频域,任取w \in [0,2\pi],定义功率谱\Phi(x)为:
    \Phi(w)\triangleq |\sum_{n=1}^Nx(n)e^{-jwn}|
    由牛顿二项式定理可将\Phi(x)变换到如下形式:
    \Phi(w)\triangleq \sum_{k=-(N-1)}^{N-1}r(k)e^{-jwk}
    对比ISL的定义式,我们可以用包含\Phi(x)的式子ISL等价式:
    {\rm ISL}=\frac{1}{2N}\sum_{p=1}^{2N}[\Phi(w_p)-N]^2\tag{1}
    其中,w_p=\frac{2\pi}{2N}p,p=1,…,2N

    若要最小化(1)所表示的ISL,就是要让目标序列的功率谱为常数,就像白噪声一样。

    优化(1)式很困难,CAN提出了改进的目标函数:
    \min_{x(n)^N_{n=1},\psi(p)_{p=1}^{2N}}\sum^{2N}_{p=1}|\sum^N_{n=1}x(n)e^{-jw_pn}-\sqrt{N}e^{j\psi(p)}|^2 \tag{2}
    (2)式越小,其优化值与(1)式越接近。

  • 优化方法

    先对x(n)进行填零,定义序列z为:
    {\bold z}=[x(1)…x(N)\space 0 … 0]^T_{2N \times1}
    z的傅里叶变化标记为\tilde z
    {\tilde z}=[\tilde x(1)…\tilde x(N)\space 0 … 0]^T_{2N \times1}
    定义\psi(p)=arg(\tilde z),arg表示角度。定v为如下形式:
    \bold v=\frac{1}{\sqrt 2}[e^{j\psi(1)}…e^{j\psi(2N)}]^T
    到此,(2)可以写成如下形式:
    ||\frac{1}{\sqrt{2N}}\tilde z-\bold v||^2
    那么,可以求出最小序列x(n)
    x(n)=e^{j \space \arg(\tilde v_n)},n=1,…,N

  • Matlab实现

    按上述分析设计算法,求出x,直到||x^{(i)}-x^{(i+1)}||<\epsilon,其中\epsilon是一个很小的值。为减小计算量,迭代初始值x_0初始化为Golomb序列。

    function x = cansiso(N, x0)
      if nargin == 2
          x = x0;
      else
          x = exp(1i * 2*pi * rand(N,1));
      end
    
      xPre = zeros(N, 1);
      iterDiff = norm(x - xPre);
    
      while (iterDiff>1e-3)
          xPre = x;
    
          z = [x; zeros(N, 1)]; 
          f = 1/sqrt(2*N) * fft(z); 
          v = sqrt(1/2) * exp(1i * angle(f)); 
    
          g = sqrt(2*N) * ifft(v);   
          x = exp(1i * angle(g(1:N))); 
    
          iterDiff = norm(x - xPre);
      end
    
  • 实验验证

    下图为长度为100的CAN序列的相关性图。其旁瓣峰值比Golomb序列更低。

    image-20200613203118036

    下图是使用MF(N^2/ISL)参数的对比图:(越大越好)

    image-20200613203335085