利用相关函数进行元素识别的一个实验性观点

随着NXP智能车竞赛的发展，赛道元素种类越来越多(比如圆环)。在以往的比赛中我们一般采用特征识别的方法来检测这些赛道元素，这种方法缺点很明显，需要为不同的特征设计不同的识别方法。本文主要探讨利用相关函数的特性构建一个通用的识别方法。

流程简介

本方法的大概流程：

• 确定识别模板序列x
• 待检测的序列y
• 对x和y序列做傅里叶变换得到X和Y序列
• 对X序列就共轭序列得到X’序列
• X’序列和Y序列逐项相乘获得频域相关序列R
• 对R序列进行傅里叶反变换得到时域相关序列r
• 检测r序列的值，若r大于设定的阈值即可认为识别成功

从相关系数到互相关函数

相关系数的意义

好了我们从头的来理解这个方法，本方法的出发点是概率论中的协方差Cov和相关系数$\rho$. 首先回顾一下自相关的意义，我个人的表述是一个数列y能由另一个数列x线性表示(y=ax+b)的程度，如果如果相关系数等于1那么y当中的每一个值都能由对应的x值来线性表示。对于协方差和相关系数的理解可以参考知乎上的这个回答https://www.zhihu.com/question/20852004?sort=created, 注意一定要理解清楚相关系数的概念，这是本方法的核心！顺便再插一句变量相关和独立的区别，变量相关指的是y能由x线性表示，即y=ax+b。但独立指的是y完全和x没有关系，完全不能用x表示，即不能用$y=ax+b,y=ax^2$或者其他任何方式表示。

相关系数的定义与相关函数

好了，假设相关系数和协方差的意义我们已经了解清楚了。那么接下来请看协方差的公式：

$$ Cov=E{[x(n)-E(x)][y(n)-E(y)]} = 1/N \sum_1^N{[x(n)-E(x)][y(n)-E(y)]} $$

Cov的大小就可表示序列的相关性

我们将相关性公式进行进一步简化得到：

$$ r=\sum_1^N{[x(n)-E(x)][y(n)-E(y)]}$$

对于信号来说可以分为直流分量和交流分量，对于直流分量可由数学期望来表示，即E(x)。对于一般信号而言，我们只考虑信号的交流分量。所以表示相关性的公式就可以进一步简化：

$$ r=\sum_1^N{[x(n)][y(n)]}$$

好了，到目前为止我们已经可以获得一个r，用来表示两个序列(信号)的相关性了。对于这个公式我再给出另外一种理解方式：相关运算从线性空间的角度看其实是内积运算，而两个向量的内积在线性空间中表示一个向量向另一个向量的投影，表示两个向量的相似程度，所以相关运算就体现了这种相似程度。

于是我们检测赛道元素就只需要先确定一组特征，然后传感器采集回来的值一直按照上面的函数计算出一个r，只要r大于一个设定的阈值，就说明元素识别成功。但是很快就会发现问题，比如x={1,2,3,4,5,0,0},y={3,4,5,0,0,1,2}.那么按照上述的方法计算出来相关性就很差，但是实际上这两组序列其实是一样的，只是相位不同而已。要解决这个问题我们就需要进一步引入相关函数的定义：

$$r_{xy}(m)=\sum_1^N{[x(n)][y(n-m)]}$$

这个函数的定义就是不断平移其中的一个序列，计算出不同相位条件下的相关性。好了这个相位差的问题看似是解决了。但是仔细计算就能发现如果两个序列很长那么计算量就会变得非常的巨大！一个序列如果长度为n，平移n次那么就需要进行$4n^2$次乘法加法计算！对于80长度的图像来说就需要25,600‬次计算，这恐怖的计算量完全无法接受！于是我们就进一步引入了相关的快速算法。

相关，卷积与快速傅里叶变换

这部分内容需要学习信号与系统，数字信号处理的课程才能够完全理解，如果实在无法理解可暂时跳过，直接看快速算法的步骤。

1.1 首先我们来看看卷积的定义：
$$X(m)\ast Y(m)=\sum_1^N{[x(m)][y(n-m)]}$$
这个和相关函数长得好像是吧！

1.2 再看看相关的公式
$$r_{xy}(m)=\sum_1^N{[x(n)][y(n-m)]}=\sum_1^N{[x(n)][y(-(m-n))]}=X(m)\ast Y(-m)$$
好了我们现在知道了卷积和相关就差一个负号。

2.1 对于学过信号与系统的同学们应该都知道时域卷积对于频域乘积，频域卷积对于时域乘积。那么既然相关运算和卷积长得如此相像，那么我们能不能把信号弄到频域去，这样一个简单的相乘就能解决卷积的复杂计算的问题了吗。所以我们现在得到了思路：

• x(n)和y(n)序列变换到频域获得X(n)和Y(n)序列
• 把Y(n)序列取共轭，即反转一下得到Y(-n)，
• 频域两个序列逐项相乘，得到频域相关函数R(n)
• 对R(n)进行反傅里叶变换得到时域相关函数r(n)
• 检测时域r(n)，判断是否满足相关的阈值

3.1 于是现在的问题就变成了如何快速的求解序列的傅里叶变换，可能大家都听说过一个叫快速傅里叶变换FFT的算法。这个算法也不是很快就能讲清楚的，反正大家只要知道这个算法能够显著的减少傅里叶变换的计算量就行了，同时感谢cortex-M3系列内核内置的CMSIS-DSP算法库已经通过软硬件结合起来帮助我们完成了这个算法。总之，FFT和IFFT(逆FFT)都已实现。
https://arm-software.github.io/CMSIS_5/DSP/html/group__ConvolutionExample.html

相关算法的具体步骤

好了到目前为止这个算法就算是大概实现了，不过由于是取值都是离散过程，和连续函数还是有一些区别。整体流程如下：

1. 给定x(n)长度为N1，y(n)长度为N2
2. 将x(n)与y(n)两个序列分别在末尾补0，让序列都变成长度为L的序列，其中$L\geq N1+N2-1$
3. 对于长度为L点的x序列和y序列分别计算FFT得到X(k)和Y(k)序列
4. 对Y(k)序列共轭得到$Y^*(k)$序列
5. 计算$R_{xy}(k)=X(k)Y^*(k)$
6. 对$R_{xy}(k)$进行IFFT得到长度为L的$r_{xy}(m)$序列
7. 计算阈值，判断相关，结束

matlab的实现代码(简易实现)

x=[1 2 3 4 5];%待检测的输入序列
y=[1 2 3 4 5];%匹配特征序列
X=fft(x,10);%做10点的fft计算，注意点数要求满足大于等于N1+N2-1
Y=fft(y,10);%做10点的fft计算，注意点数要求满足大于等于N1+N2-1
Y1=conj(Y);%对Y取共轭,频域都是复数
R=X.*Y1;%频域乘积运算
r=ifft(R,10)%转换到时域输出相关函数

输出结果：

r =

   55.0000   40.0000   26.0000   14.0000    5.0000    0.0000    5.0000   14.0000   26.0000   40.0000

可以看到不同的相位情况下相关性不同，相差的值非常大，如果序列更长那么相关性会更好。