粒子滤波算法的FPGA设计与实现

毛丽民，刘叔军，朱培逸，张水平

粒子滤波算法是一种基于蒙特卡罗方法和递推贝叶斯估计的滤波方法，主要用于非线性、非高斯系统，发展于20世纪90年代.粒子滤波虽然在理论上取得了很多研究成果，但在硬件实现方面，仅仅处于初级阶段，在实际信号处理领域还没有得到广泛的应用，主要是因为粒子滤波算法的结构复杂、运算量大，需要很大的存储空间［1］.虽然有些改进的粒子滤波算法在一定程度上提高了粒子滤波算法精度，但使得算法更加复杂，实现更加困难.

FPGA不仅存储容量大，而且能真正实现粒子滤波的并行运算，提高运算速度［2］.本文基于FPGA的粒子滤波算法硬件实现，以二维纯方位目标追踪为例验证了基于FPGA的粒子滤波算法，不仅在一定程度上简化粒子滤波算法，节省了硬件存储资源，而且提高了算法的运行速度.

1 粒子滤波算法简介

粒子滤波算法原理如图1所示［3］.图中粒子滤波算法主要由时间更新、观测更新和重采样3步组成.k-1时刻的先验概率由N个权值为1/N的粒子近似表示.在时间更新中，使用状态转移方程预测每个粒子在k时刻的状态.使用观测值更新粒子权值.在重采样过程中权值小的粒子将被抛弃，保留权值较大的粒子，权值被重新设为1/N.

图1 SIR算法示意图

2 ISE System Generator与Matlab联合开发FPGA

ISE System Generator提供了一种软/硬件相结合的设计环境，与传统的编写代码的开发模式不同，而且提供了丰富的模块库，只要能用Verilog语言实现功能都用系统提供的模块实现，并且模块中包含很多已经成熟的算法，可以减少工作量，提高开发效率［4］.

3 简化粒子滤波算法FPGA实现

基于FPGA的纯方位目标跟踪算法，如图2所示，主要由4部分组成:粒子存储及采样模块、粒子权值和阈值计算模块、粒子的重采样模块、粒子的状态估计及输出模块.每个变量都采用定点数表示，本文采用18位表示位置信息，用16位表示相位信息［5］.

图2 粒子滤波算法的实现框图

3.1 存储采样模块

存储模块主要用于存储采样前粒子的状态，其具体实现框图如图3所示.

图3 存储模块实现框图

图3 中xk RAM用于存储粒子X方向的位置，vxk RAM用于存储粒子X方向的速度.采样模块就主要将噪声生成模块产生的噪声和速度值相加;同时，噪声向右移一位后加上位置的更新值.其实现框图如图4所示［6］.

图4 采样更新单元结构

本文设计的噪声生成模块是根据查找表的噪声生成方式实现的，具体实现结构如图5所示.该图中，噪声生成模块主要由一个线性反馈移位寄存器(LFSR)与两个高斯白噪声生成器组成，其中LFSR用于生成L-bits的伪随机数，高斯白噪声生成器［1］用来产生零均值方差为1的高斯白噪声.图5表示噪声生成模块所生成的X方向的噪声随机数.

图5 噪声生成模块

高斯白噪声生成器的内部结构如图6所示:主要是由4个Box-Muller和3个加法器组成，用来生成零均值方差为1的高斯白噪声［7］.

图6 高斯白噪声生成器内部结构

3.2 权值及阀值计算模块

本文选取重要密度函数为先验概率密度［8］，权值更新公式简化为:

由于零均值的高斯白噪声为观测噪声，所以权值的计算公式为:

如果 N^eff值较小，k-1时刻进行粒子重采样，则=1/N;如果 N^eff值较大，k-1时刻不进行重采样，则=1/N，将权值归一化处理，

可以忽略，所以权值计算公式简化为

图7为简化粒子的权值计算实现框图.式中exp()与arctan()将用到Xilinx提供的高效IP核.位置信息通过IP核计算预测的观测值，然后与观测值比较，得到差值，依据公式3需要对差值进行平方.当差值大于255时，将该粒子的权值设为零，此时简化了硬件结构.delay是arctan IP核的运行时间，delay1是指数计算时间，exp指数运算的实现框图如图8所示［9］.

图7 权值计算模块框图

图8 指数运算实现框图

在计算中，指数计算由ex=shx+chx表示，其中sh与ch的范围为［-π/4，π/4］.因此，将指数运算的数值分为整数部分和小数部分.小数的指数运算由IP core cordic获取，整数使用查找表方法.最后将整数与小数相乘得到粒子的权值.

3.3 重采样模块

图9为简化重采样模块结构图.因为丢弃的粒子与选取的粒子数目是未知的，故将Single Port RAM2和Single Port RAM1的深度都设为256.

图9 简化重采样模块结构图

3.4 状态估计输出模块

状态估计值输出模块(图10)在使用FPGA实现时，与存储采样模块、权值及阀值计算模块并行执行，在存储采样模块与权值和阀值计算模块获得新的粒子以及权值后通过乘累加器来实现.

图10 状态估计值输出结构

要得到最终的状态估计值，还需对每个粒子及其权值的乘累加值进行比例调整.比例调整就是将各个状态值除以Wk和权值.二维目标跟踪中，有4个状态值，因此需要四次除法过程，但除法器占用Slice资源比较多.因此采用4个乘法器，先求出权值和的倒数l/Wk，最后将每个状态值输出相乘.

4 实验结果

System Generator把Xilinx Blockset的元件生成相应的IP核库组件，从系统设置的参数中推断出电路及控制信号，并把Simulink的模型对应生成ISE工程.为了确保多速率系统汇编的正确性，System Generator为执行工具生成了一些约束文件.很多信息都能在System Generator模块中配置相关参数.

将各个功能模块连接起来，在System Generator中配置好相关参数后，转换成硬件描述语言VHDL，下载到 FPGA，然后对整个模块进行调试.系统的资源使用情况如表1所示.

表1 资源使用情况对比

最后将FPGA实验结果导入Matlab，与Matlab仿真结果进行对比，对比结果如图11所示.

图11 粒子滤波算法的仿真结果图

图11 是系统运行50个时刻的某一次仿真结果对比，选取粒子数为256个.从图中可以看出，FPGA实验结果与真实轨迹比较吻合，已能较好地反映真实轨迹的走势.

5 结语

粒子滤波算法是一种以贝叶斯估计为依据的滤波算法，广泛用于非高斯、非线性、时变系统的滤波及参数估计.但是，粒子滤波算法结构复杂，运算量大，限制了在实时系统中的实际应用.本文在研究算法的基础上，重点研究了基于FPGA的硬件实现方法，给出了一种简化粒子滤波算法及在FPGA上的实现方案，并给出了整体结构及各个功能模块;最后进行下载调试，达到了预期效果.

［1］胡士强.粒子滤波原理及其应用［M］.北京：科学出版社，2010：1-7.

［2］刘刚，梁晓庚.遗传重采样粒子滤波的目标跟踪研究［J］. 计算机工程与应用，2010，46(19)：196-199.

［3］高展宏.多媒体处理FPGA实现［M］.北京：电子工业出版社，2010：67-70.

［4］洪少华，史治国.用于纯方位跟踪的简化粒子滤波算法及其硬件实现［J］.电子与信息学报，2009，31(1)：96-100.

［5］赵梅，张三同.辅助粒子滤波及其仿真应用［J］.北京交通大学学报，2006，11(2)：24-28.

［6］王来雄，黄士坦.一种新的粒子滤波算法［J］.武汉大学学报，2009，39(1)：118-120.

［7］邓文坛，张三同.一种改进的粒子滤波跟踪算法的研究［J］. 自动化技术与应用，2008，27(3)：84-87.

［8］唐涛，张涯辉.粒子滤波在被动定位跟踪中的应用［J］. 仪器仪表学报，2005，26(8)：211-214.

［9］刘刚，梁晓庚.遗传重采样粒子滤波的目标跟踪研究［J］. 计算机工程与应用，2010，46(19)：196-199.