摘 要：二维傅立叶变换" title="傅立叶变换">傅立叶变换" title="快速傅立叶变换" title="快速傅立叶变换">快速傅立叶变换">快速傅立叶变换(FFT)在一个传统概念的处理机上实现时，需要芯片具有更多的逻辑资源。本文给出了基于FPGA的自定义处理机（CCM）的二维FFT算法和实现。在CCM的Splash-2平台上实现了二维FFT，计算速度达到180Mflops，最快速度超过Sparc-10工作站的23倍。同时，对于一个N×N图像，这种实现方法可以满足二维FFT所需要的O(N²log₂N)次的浮点算术运算。
关键词：自定义处理机（CCM）；二维；快速傅立叶变换（FFT）； 图像；Splash-2

　　当被处理的图像像素相对较小时，可以直接采用空间域滤波的方法。但随着像素的增多，计算量也随之增大。针对这种情况，本文对像素较大的图像，采用了一种将图像先从空域转换到频域，执行点到点的乘法滤波，然后通过反傅立叶变换，将频域内滤波器图像再转换到空间域的方法。在此基础上，对二维快速傅立叶变换采用了基于FPGA的CCM的Splash-2平台实现方式，目的是为了对视频信号进行实时滤波。实现过程中需要对浮点数" title="浮点数">浮点数进行计算，采用了CCM可以自定义数据格式和算法数据流以满足应用的需要。
1 用傅立叶变换进行图像滤波
　　图1是用快速傅立叶变换进行图像滤波处理的实现过程^[1]，其中，预处理阶段包括确定图像大小、获得填充参数和生成一个滤波函数等过程；后处理阶段包括计算结果的实部" title="实部">实部，修剪图像以及将图像进行转换等操作，其他部分将在以下环节分别进行介绍。

1.1 离散傅立叶变换（DFT）
　　一个N×N图像的二维的DFT，f(x，y)定义为：

　　二维的DFT变换可以分解成多个一维的傅立叶变换，则(1)式等同于：

　　公式(1)、(2)表明一个N×N的二维DFT可以通过先进行N个一维DFT（行），然后再进行另外N个一维DFT（列）来计算，即所谓的行列分解法^[2]。
1.2 快速傅立叶变换（FFT）
　执行一个N点DFT所需的复杂乘法和加法运算次数与N2成正比。由于一维DFT可以被分解，所以进行乘法与加法运算的次数正比于N×log₂N。而FFT算法通过“分而治之”的方法来实现其较高的计算效率。通过时间与频率采样的分组，具有N个值的DFT运算可以表示成两点DFT的组合方式。这两点DFT称为蝶形运算，它需要一次复乘和两次复加来实现。蝶形结构的符号如图2所示（左侧：基2的蝶形结构）。

　　用FFT“分而治之”的方法，一个8点的FFT的计算如图2所示。N点FFT的每一级运算由N/2基2的蝶形组成，总共有log₂N次运算。因此，每个FFT共有(N/2)×log₂N个蝶形结构。一个二维FFT可以拆分为两组一维的DFT运算，且每个DFT可以用一维的FFT来计算。此外，数据以倒位序输入而以正常的线性顺序输出。
2 浮点算法
2.1 浮点格式的表示
　　为了实现本文给出的FFT运算，在此设计了一个较小的18bit的浮点格式。用这种格式来满足两种特定要求：(1)动态范围需要非常大，以便能够准确地表示一个实数的大小、正负；(2)进入Splash-2平台的XILINX 4010处理器的数据通道宽度为36bit，且每一个时钟周期" title="时钟周期">时钟周期内都输入复数的实部与虚部。基于这些要求，使用的浮点格式如图3所示。

2.2 浮点加法减法与乘法
　　为了在每个时钟周期产生一个结果，开发浮点加减法程序来实现管道化单元，所实现的浮点加法与减法算法与最传统的处理器相似。
　　浮点乘法与整数乘法类似。因为浮点数是以“符号－数字”的形式存储，乘法器仅需要处理无符号整数，与浮点加法器的结构类似^[3]，浮点乘法器单元也是每个时钟周期产生一次结果，这种设计的瓶颈是整数乘法器。
　　本文将VHDL语言放在XILINX芯片相应开发工具ISE中进行编译，同时给出时延、速度等相关参数。
　　浮点运算单元还与FIR滤波器组合一起使用。FFT运算以10MHz工作，转换的结果存储在Splash-2开发板的存储器中。算术单元的最大时钟速度至少可以工作在10MHz。
3 FFT的实现
　　本文从一个帧缓冲区得到连续的视频图像提供给FFT和IFFT参与并行计算从而构建出滤波算法。以下讨论用于实现FFT、FFT中的蝶形操作以及滤波处理的再循环算法。
3.1 FFT再循环算法
　　为了计算二维的FFT，需要设计一个通过蝶形运算的数据循环算法。图4给出了相应的结构图。

　　两列存储器用来存储FFT的每一级运算的输入输出数据。每一列存储器包括三个处理单元，用来把两个18bit浮点数的实部与虚部存储到它们的本地内存。由于本地内存只有16bit宽，先将两个18bit浮点数值进行分割，然后按顺序放在三个存储单元中。为了计算输入图像的FFT，来自于帧缓冲区中的数据帧从8bit的整数被转换为18bit的浮点数存储在第一列中。在计算二维FFT时，首先计算图像每一行的一维FFT，然后计算行变换后每一列的FFT。一维FFT的计算与图2中所示方法相同。FFT的第一级通过从第一列存储区以倒位序读取数据采样点的每一行，然后将它送给蝶形操作进行计算。蝶形计算的结果存储在第二列存储区。当第一级中的蝶形计算完成时，开始第二级计算，从第二列存储区中顺序读取数据然后送入蝶形处理器，这一操作的结果存储在第一列存储区。这一循环方法通过从一个存储区中读数而在另一外存储区存储计算结果来进行。当图中每一行一维FFT计算完时循环终止。第二组一维FFT与第一组以相同的方法进行计算，直到第二组的FFT计算完第一组的所有结果。当最后一次FFT的最后一级计算完时，数据从X11传到X15进行滤波。一次完整的二维FFT过程包括2×N²×log₂N次蝶形操作。
3.2 蝶形实现
　　蝶形操作是FFT算法的核心。蝶形运算框图如图2所示，它包括了一次复乘与两次浮点加法、减法。复乘包括4次乘法和两次加减。每个时钟周期内以10MHz进行8次浮点操作。蝶形操作的计算量为80Mflops。图5表示在Splash-2平台上如何将蝶形运算在5个处理单元之间进行分割计算。
复乘的实部与虚部分别用下式来表示：

　　图2中蝶形运算输入A与输入B如图5中的f（x）所示。A值先送入计算，在下个时钟周期再送入B的值。输入A不与旋转因子相乘而是通过多路选择开关对它乘1加0无改变地通过循环运算。当B的实部与虚部送入时，它们依次通过4个处理单元来得到复乘结果。处理单元1（PE1）读出相应旋转因子的实部并将它与B的实部相乘，计算结果与旋转因子送到处理单元3（PE3），B的实部与虚部送到处理单元2（PE2）。PE2将B的虚部与相应旋转因子相乘结果送入PE3。PE3从PE1中读取B.re.re，把它与从PE2中读取的B.im.im相加，得到复乘实部的最终结果.re；复乘虚部.im可用同样的方法在PE3与PE4中得到。在第一个时钟周期内将A与相加得到X，在第二个时钟周期内将A减去得到Y，这样就完成了蝶形运算。

　　由于存储器数据总线的宽度仅为16bit，在本地PE1、PE2的存储器中无法存储18bit格式的旋转因子。因此从18bit浮点指数域中减去2bit来产生一个较小的16bit格式，而用欧拉定理可以将旋转因子表示成正余弦函数的形式，因此浮点数的值将不会有大于0的指数。因此，指数域范围从0～－31变化，这样指数域的长度从7bit减到5bit。当旋转因子被读入处理单元时，在相应运算单元中进行一次从16bit格式到18bit格式的转换即可。
3.3 滤波
　　当输入图像转换到频域时，矩阵滤波系数H(u，v)和转换后的图像可以通过滤波器得到滤波后的图像。矩阵H(u，v)元素值的范围为0～1.0，它与蝶形运算旋转因子以相同的的存储方式存储在本地的存储器中。滤波系数存储在X15和X16芯片的本地存储器中。滤波器芯片包括一个浮点乘法单元和滤波系数地址逻辑单元，X15与X16用来对转换后图像的实部和虚部分别进行滤波。
　　不同类型的滤波器如理想滤波器以及其他类型滤波器可以先计算，然后从宿主机下载到开发平台的存储器中。
　　由于CCM的灵活性，浮点格式可以通过自定义方式以最小的比特数来达到最大精度。利用Splash-2的并行结构，地址计算、蝶形运算以及滤波可以并行地进行计算。通过蝶形运算操作，每个时钟周期内可以以10MHz的速度得到一个实数或复数结果。这种应用的性能类似于一个典型的DSP处理器所达到的性能。

参考文献
[1] RAFAE1 C G，RICHARD E W，STEVEN L E.数字图像处理.北京：电子工业出版社，2006.
[2] DAN E D，RUSSELL M M.多维数字信号处理.北京：科学出版社，1991.
[3] SHIRAZI N，WALTERS A P A.Quantitative analysis of floating point arithmetic on FPGA based custom computing machines.IEEE Symposium on FPGAs for Custom Computing Machines，1995，(4).