基于FPGA的快速并行FFT及其在空间太阳望远镜图像锁定系统中的应用

摘要：在空间太阳望远镜的在轨高速数据处理中，运算时间是影响系统性能的重要环节之一。利用FPGA丰富的逻辑单元实现快速傅里叶变换（FFT），解决 了在轨实时大数据量图像处理与航天级DSP运算速度不足之间的矛盾；利用溢出监测移位结构解决了定点运算的动态范围问题。经过实验验证，各项指标均达到了设计要求。 关键词：FFT FPGA 蝶形运算 空间太阳望远镜项目是我国太阳物理学家为了实现对太阳的高分辨率观测而提出的科学计划。它可以得到空间分辨率为0.1"的向量磁图和0.5"的X射线图像，实现这样高的观测精度的前提就是采用高精度的姿态控制系统和高精度的相关跟踪系统。从整个系统来看，相关运算所需的时间成为限制系统性能能否提高的一个重要环节。 目前，国际国内相关计算比较通用的实现方法有两种：用高速DSP或者专用（FFT）处理芯片。用DSP完成相关计算（关键是FFT）受到航天级DSP性能的限制，现有的航天级DSP（如ADSP21020）计算一个32%26;#215;32点8bit的二维FFT所用时间需要1.5ms以上，远远不能满足系统设计要求；而现有的FFT处理芯片在处理速度、系统兼容性、抗辐射能力等方面不能满足空间太阳望远镜所提出的要求。 为克服这一矛盾，本文利用FPGA资源丰富、易于实现并行流水的特点设计专用的FFT处理芯片来完成复杂的、大量的数据处理；并通过在运算中作溢出监测来保证定点运算的精度，从而大大缩短系统的响应时间，将极大地提高空间太阳望远镜的在轨实时图像处理能力；同时由于FPGA具有抗辐射能力，可以提高系统的可靠性，其在航天遥测遥感和星载高速数据处理等方面将有广泛的应用前景。 1 算法构成 1.1 FFT算法选择 提高FFT速度的两个主要途径是采用流水结构和并行运算[1]。采用高基数结构也可以提高速度，只是用FPGA实现时必须综合考虑系统要求、结构特点及片内资源。针对本系统自身特点，这里按时间抽选算法进行分析。由于32不满足N=4m，所以32点FFT算法不能采用基-4 FFT运算。当详细分析基-2蝶形图时，有些蝶形运算并不需要做乘法，例如 等[2]；对于32点DIT-FFT，一共80个蝶形运算，这种结构就有46个，极大地降低了运算复杂度。在一维FFT计算效率提高的基础上对二维FFT采用最常用的行列算法[3]，综合各项指标本系统采用基-2 DIT行列算法。 1.2 算术运算方案 本系统是针对32%26;#215;32点16bit的二维图像进行快速傅里叶变换（FFT），设计要求运算在0.5ms之内完成，所以采用定点运算更符合系统对时间的要求。对于定点运算，必须用定比例的方法防止溢出，即必须解决动态范围问题。下面对其进行理论分析： 若{x(n)}是-N点序列，其DFT为{X(K)}，由Parseval定理得[4]： 由式（1）可知变换结果的均方值是输入序列均方值的N倍。考虑基-2算法的第m级蝶形运算，用Xm(i)、Xm(j)表示原来的复数，则新的一对复数Xm+1(i)、Xm+1(j)为： Xm+1(i)=Xm(i)+Xm(j)%26;#215;W (2) Xm+1(j)=Xm(i)-Xm(j)%26;#215;W 其中，W为旋转因子。首先，考虑复数的均方程根值。由（2）式可得： 因此，从均方根意义上看，数据（实数或复数）复级都增加（2的平方根）倍。其次，再考虑复数的最大模。由（2）式可以证明[5]。 max{|Xm(i)|,|Xm(j)|}≤max{|Xm+1(i)|,|Xm+1(j)|}≤2max{|Xm(i)|,|Xm(j)|} 因此，复数数组的最大模是非减的。所以，对于DITFFT，其每一级的蝶形运算之后数值都会增加1+(2的平方根)≈2.414倍。在每一次运算完成之后，须将结果右移2bits以满足要求。 2 系统实现 系统原理如图1所示，整个FFT运算处理单元分为三部分：存储单元（两个输入/运算存储器、一个输出存储器及旋转因子存储器）、蝶形运算单元、地址产生器。 2.1 存储器 本系统实时接收前端CCD相机的图像。为保证CCD相机采集图像的准确率，图像的每一行、每一帧之间都必须有一定的时间间隔,故采用两个存储单元作为输入数据和中间数据的暂存单元（如图1所示），以节省时间实现实时处理。当系统工作时，将图像存入存储器、计算上一次采集的图像、将存储器中的结果输出，这三个工作同时进行，用简单的流水方式减少存储数据所需的时间。旋转因子则预先存储在器件的内置ROM中。根据级数不同选用不同的因子。 2.2 蝶形运算单元 一个基-2蝶形运算由一个复乘和两个复加（减）组成，采用完全并行运算，进一步分解为四个实数乘法，六个实数加（减）法，分三级并行完成，加上前后输入输出的数据锁存，共需要6个时钟周期。32点的FFT需要16%26;#215;5=80个基-2的蝶形运算，一幅图像一共是32行32列，不考虑不需要做乘法的蝶形运算，一路串行共需要6%26;#215;80%26;#215;32%26;#215;2=30720个时钟周期，采用频率为10MHz的时钟，即为3ms。对于蝶形运算的第一、第二级都可以由不带乘法器的蝶形结构来实现同步并行运算，每一个蝶形运算加上前后的数据锁存仅需4个时钟周期即可完成；对于第三、第四、第五级，由于带乘法器不带乘法器的两种蝶形运算结构同时存在，必须加入等待时间才可以实现严格同步。同时由于各级计算时间不同，所以不能实现深度流水。因此，采用多路并行及部分流水，在时间上即可满足系统要求。 上面讨论了当运算从一级转到另一级时，序列中数值的幅度一般会增大。因而，运算方法是在内循环中作溢出监测。如果没有溢出，则计算照常进行；若有溢出，则把产生溢出的数据右移，一直到没有溢出为止。记录下移位的次数（0、1或2），并把整个序列右移同样位数，移位总数进行累计，累计数的负值作为2的幂，由此得出最终序列的总的比例因子。比例因子s由下式定义[6][7]： 这里bi为比例参数。 k=0,1,2,…，N-1 （6） 根据公式(6)，FFT的最终结果要除以比例因子。式中x(n)为原始数据，X(k)为除以比例因子之前的结果，X"(k)为最终结果，1/s为比例因子的倒数。 如图2所示，对于一个基-2蝶形单元，当从存储器中读取的Bbit输入数据进入蝶形运算单元PE1后，经过乘法运算（MU1）乘以旋转因子，数据变为（B+Bω）bit，然后作加（减）法，得到蝶形运算结果（B+Bω+1）bit。为防止溢出，进行移位操作。M1、M2为比例选择器，根据不同的级数，选择不同的比例因子。最后，输出数据再放回到存储器中。 3 FPGA器件选择 本设计采用XILINX公司的VERIEX系列XCV300-4HQ240芯片。该芯片有丰富的可配置逻辑模块CLBs(Configurable Logic Blocks)、大量的触发器以及内置的不占系统资源的块RAM。系统最大工作频率可达200MHz，兼容多种接口标准，有相应的航天级产品，是目前市场上为数不多的能达到此项要求的高性能可编程逻辑器件。 VERTEX系列器件的一个显著特点是内置的延迟锁相环DLL（Delay-Locked Loop），它可以减少时钟传输的衰减，每一个DLL可以驱动两个全局时钟信号。DLL可以倍频，或者1.5、2、2.5、3、4、5、8以及16分频。VERTEX系列器件内部的4-输入查找表LUTs(Look-Up Tables)也具有多种功能：可以作为16 %26;#215;1bit的同步RAM，而且一个块（Slice）中的两个LUTs可以组合成一个16%26;#215;2bit或者一个32%26;#215;1bit的同步RAM或者一个16%26;#215;1bit的同步多口RAM。另外，LUTs还可作为一个16bit的移位寄存器使用，该寄存器用来获取高速或者突发数据非常理想，特别适用于数字图像处理中的数据存储[8]。 本设计充分利用了VERTEX器件的LUTs替代触发器和基本门电路搭建乘法器和加法器这两个显著的结构特点，节省大量触发器资源，避免了缺少触发器而LUTs大量剩余的尴尬；增加了器件利用率、布通率，降低布线延迟。由于本系统最终用于空间太阳望远镜，所以板上时钟频率不可超过20MHz。但基于地面测试的需要，特利用DLL对外部时钟信号进行了倍频，以提高芯片内部的运行速度。 本设计利用FPGA易于实现并行运算的特点实现专用的FFT处理芯片，解决了在轨实时大数据量图像处理与航天级DSP运算速度不足之间的矛盾，提高了系统实时处理能力。两维FFT不到400μs即可完成，高于航天级DSP(ADSP21020)1.5ms的处理速度。对太阳米粒组织图像进行处理（实验数据如表1所示），结果显示数据误差都在1%左右。这样的误差满足空间太阳望远镜中的相关摆镜的系统要求。实验证明用高性能FPGA实现空间化的FFT处理芯片是完全可行的。 表1 实验数据 序 号原始数据浮点数据定点数据误差% 序 号原始数据浮点数据定点数据误差% 1 -32768 420104 419328 0.18 17 19072 5048 5120 1.43 2 -32768 154298 155008 0.46 18 18704 5943 6016 1.23 3 7216 97679 98304 0.64 19 18796 9445 9600 1.64 4 14800 106224 106688 0.44 20 19152 26885 27200 1.17 5 16080 93313 93760 0.48 21 16304 31028 31232 0.66 6 19840 76389 76800 0.54 22 14752 33592 33600 0.02 7 23104 65472 65672 0.68 23 14576 38986 39168 0.47 8 20000 604403 60672 0.45 24 15200 58786 59136 0.60 9 19104 69973 70400 0.61 25 8160 39973 70400 0.61 10 23376 58786 59136 0.60 26 7048 60403 60736 0.55 11 20128 38986 39104 0.30 27 14384 65472 65920 0.68 12 18688 33592 33600 0.02 28 11584 76389 76800 0.54 13 19680 31028 31232 0.66 29 11472 93313 93760 0.48 14 19536 26885 27200 1.17 30 11456 106224 106752 0.50 15 19680 9445 9600 1.64 31 11760 97679 98304 0.64 16 18720 5943 6016 1.23 32 11488 154298 154944 0.42