CN103106181A - 一种大点数fft在处理器上的实现方法 - Google Patents

Abstract

本发明公开了一种大点数FFT在处理器上的实现方法，能够解决传统FFT算法在处理器上实现大点数快FFT时没有充分考虑Cache丢失对执行效率影响的问题，改进了传统Winograd算法处理速度有限的问题。该方法包括：将一维序列存储为二维矩阵；处理器先列FFT：每次从二维矩阵中读取i列数据，读取的i列数据分

次处理，则处理器共读取并处理次；其中，在保证列长度的基础上，使行长度M小于或等于处理器所用Cache的容量；处理器再进行行FFT，一次一行，且采用新的旋转因子，并将结果按照列方向输出。

Description

一种大点数FFT在处理器上的实现方法

技术领域

本发明涉及一种改进型快速傅立叶变换（FFT）算法在处理器上实现大点数FFT的方法，属于信号处理领域。

背景技术

快速傅里叶变换(FFT)广泛应用于雷达、通信和图像处理等科学技术领域，这使得FFT的工程实现具有十分重要的意义。特别是，在雷达系统中高分辨大测绘带宽的合成孔径雷达的飞速发展，对信号处理系统大数据的高速实时处理提出了更高的要求，这就需要信号处理中大点数FFT的快速实现。在实际应用中，一般采用专用数字信号处理器(DSP)来实现。

TS201是美国模拟器件公司的一款高性能、高并行的静态超标量处理器。在TS201处理器中，其内部嵌入了24Mbit的嵌入式DRMA，整个DRAM划分为6个存储块，每个存储块都通过交叉连接器分别连接4套128位宽的内部总线，因此处理器能够在同一个时钟内实现对4个存储块的访问。这些交叉连接器包含预取数缓冲、读缓冲、回存缓冲和高速缓冲，其连接图1所示。TS201通过地址总线和数据总线对DRAM进行读写数据操作时，首先会将数据缓存到缓冲区（Cache）内，内核在读取数据时会先直接从Cache中读取数据，如果不能在Cache中命中数据，再从DRAM中读取数据。因此，通过Cache的预缓存作用可以提高内核对DRAM的读写效率。但是Cache的大小有限，当进行大点数FFT处理时，Cache不能容纳整个序列的数据，那么一部分在Cache中另一部分在DRAM中，将会带来读取速度较慢，以及访问错误等问题。

现有多种FFT算法在TS201上实现，例如Winograd算法。现将该算法介绍如下：设FFT变换前的序列为x(n)，FFT变换后的序列为X(k)：

$X\left(k\right)=\underset{n=0}{\overset{N-1}{\mathrm{\Σ}}}x\left(n\right)W_N^{\mathrm{nk}},$ k,n＝0,1,...,N-1  (1)

其中，W_n为铰链因子，W_n=e^-j2π/n，

N为序列中元素总数，FFT变换前后序列元素数目不变。

传统的Winograd算法实现FFT的思想是将大点数的FFT尽量拆分成小点数来计算，将一维序列x(n)和X(k)分别在时域和频域映射成二维矩阵形式。以时域x(n)为例，将其拆分为L×M的二维矩阵，L为行数，M为列数，那么FFT变换后频域X(k)表示为M×L的二维矩阵。

设n₁和n₀分别为时域二维矩阵的行、列序号，k₀和k₁分别为频域二维矩阵的行、列序号，则有如下关系：

$\left\{\begin{array}{c}n={\mathrm{Mn}}_1+n_0;k={\mathrm{Lk}}_1+k_0\\ n_0=\mathrm{0,1},...,M-1;n_1=\mathrm{0,1},...,L-1\\ k_0=\mathrm{0,1},...,L-1;k_1=\mathrm{0,1}...,M-1\end{array}\right.---\left(2\right)$

将（2）式代入公式（1）中得：

$X\left(k\right)=X(k_1,k_0)=\underset{n_0=0}{\overset{M-1}{\mathrm{\Σ}}}\underset{n_1=0}{\overset{L-1}{\mathrm{\Σ}}}x({\mathrm{Mn}}_1+n_0)W_N^{({\mathrm{Mn}}_1+n_0)({\mathrm{Lk}}_1+k_0)}$

$=\underset{n_0=0}{\overset{M-1}{\mathrm{\Σ}}}\left(\right\{\underset{n_1=0}{\overset{L-1}{\mathrm{\Σ}}}x(n_1,n_0)W_L^{n_1{k}_0}\left\}{W}_N^{n_0{k}_0}\right)W_M^{n_0{k}_1}---\left(3\right)$

其中，

与式（1）的形式类似，相当于对一点数据进行FFT，设 $X_s\left(n_0\right)=\underset{n_1=0}{\overset{L-1}{\mathrm{\Σ}}}x(n_1,n_0)W_L^{n_1{k}_0},$ 则 $\underset{n_0=0}{\overset{M-1}{\mathrm{\Σ}}}(X_s\left(n_0\right)\·W_N^{n_0{k}_0})W_M^{n_0{k}_1}$ 也与式（1）的形式类似，相当于对X_s(n₀)乘以铰链因子

后再进行一次FFT。

可见，由公式（3）可以得到Winograd算法实现FFT的步骤：

1)将一维序列x(n)拆分成L×M的二维矩阵序列，并转置为M×L；

2)对二维矩阵，按照行方向计算L点FFT，共处理M条线；

3)对步骤2)的变换结果乘以铰链因子

4)对步骤3)的处理结果进行转置，变换为L×M的矩阵；

5)对步骤4)的变换结果，按照行方向计算M点FFT，共处理L条线；

6)将处理结果进行转置，得到计算结果。

该Winograd算法的问题是，需要单独执行一个乘以铰链因子的步骤，并进行3次显性转置，从而引入了额外运算，因此会降低处理速度。

发明内容

有鉴于此，本发明提供了一种大点数FFT在处理器上的实现方法，能够解决传统FFT算法在处理器上实现大点数快FFT时没有充分考虑Cache丢失对执行效率影响的问题，并且通过优化矩阵行列读取方法，并进行蝶形运算重构，改进了传统Winograd算法处理速度有限的问题。

该大点数FFT在处理器上的实现方法包括如下步骤：

步骤一、将待处理的一维序列x(n)分L段存储为L×M的二维矩阵，L为列的长度，M为行的长度；设定后续步骤二中每次读取i列数据，i为正整数，则在保证

的基础上，使行长度M小于或等于CacheLength；CacheLength为处理器所用Cache的容量；

步骤二、处理器进行列FFT；

处理器每次从L×M二维矩阵中读取i列数据，通过Cache的缓存放入内存的指定空间，然后从指定空间读取数据并进行列FFT，并将结果原位存回L×M二维矩阵；假设根据处理器数据宽度的限制，处理器每次处理数据量为w列，则i的取值为w的整数倍；读取的i列数据分

次处理；处理器共读取次数据并进行列FFT，从而实现M次L点的列FFT；

步骤三、处理器进行行FFT；

处理器每次从步骤二处理后的L×M二维矩阵中读取一行数据，通过Cache的缓存放入内存的指定空间，然后从指定空间读取缓存数据并进行行FFT，并将结果按照列方向输出；处理器共读取L次数据并进行行FFT，从而实现L次M点的行FFT；

本步骤的行FFT运算时第b级蝶形运算所用的旋转因子W(b,u)由下式确定：

$W(b,u)=W_{P\left(b\right)}^{k_0}\·W_{Q\left(b\right)}^u$

其中，W_n=e^-j2π/n，P(b)＝N/2^c-b，Q(b)＝M/2^c-b；

b表示FFT算法中蝶形运算的当前级数；

c表示FFT算法中所包含蝶形运算的总级数，c＝log₂(M)；

u表示b级蝶形运算输出序列的序号，取值范围是u＝0,1,...,Q(b)-1；

k₀为当前进行FFT变换的行数据的行序号，即FFT变换后得到的频域二维矩阵的列序号。

有益效果：

本发明提供了一种应改进型Winograd算法处理器上实现大点数FFT的方法，与现有的FFT实现方法相比，具有以下优点：通过优化矩阵行列读写效率，来避免数据读写时Cache的频繁丢失对FFT执行效率的影响，同时采用改进型Winograd算法进行蝶形运算重构将乘铰链因子隐藏在蝶形运算中，减少了数据相乘带来的时间开销，因此，本发明方法可以显著提高了FFT的运行效率。

附图说明

图1为TS201处理器内部框图。

图2为矩阵列方向访问时的示意图。

图3为矩阵行方向访问时的示意图。

图4为本发明流程图。

图5为本发明方法与现有几种算法实现不同浮点数FFT时，执行时间的曲线对比图。

具体实施方式

本发明通过以下两个方面对现有技术进行改进：

（一）相对于传统Winograd算法中引入单独一个乘以铰链因子的步骤的问题，本发明采取重构蝶形运算将铰链因子

运算隐藏在第二次FFT处理中。

（二）相对于传统Winograd算法有三次显性转置的问题，本发明改变对二维矩阵的读取顺序来实现转置，但考虑到按列读取时，行切换会带来额外时间开销，则本发明通过优化矩阵行列读写效率来避免，且行列读写时尽量充分利用Cache，从而提高Cache命中率，以提高FFT执行效率。

下面结合附图对上述两个改进点进行详细说明，并且采用TS201处理器为例。在实际中，本发明不限于采用哪种处理器。

（一）将铰链因子

运算隐藏在第二次FFT处理中的方案设计。

由于FFT算法中包括多级蝶形运算，每次蝶形运算都需要乘以旋转因子，因此，可以通过重构FFT算法中的蝶形运算将第（3）步和第（5）步进行合并，只是蝶形运算所乘的旋转因子不同。新旋转因子的推导过程如下：

设：

$X_s\left(n_0\right)=\underset{n_1=0}{\overset{L-1}{\mathrm{\Σ}}}x(n_1,n_0)W_L^{n_1{k}_0}---\left(4\right)$

将式（4）代入公式（3）有：

$X(k_1,k_0)=\underset{n_0=0}{\overset{M-1}{\mathrm{\Σ}}}\left(X_s\left(n_0\right)W_N^{n_0{k}_0}{W}_M^{n_0{k}_1}\right)---\left(5\right)$

对公式（5）的变换序列按照奇偶分解，可表示为：

$X(k_1,k_0)=\underset{h=0}{\overset{M/2-1}{\mathrm{\Σ}}}{X}_S\left(2h\right)W_N^{{2\mathrm{hk}}_0}{W}_M^{{2\mathrm{hk}}_1}$

$+\underset{h=0}{\overset{M/2-1}{\mathrm{\Σ}}}{X}_S(2h+1)W_N^{(2h+1)k_0}{W}_M^{(2h+1)k_1}---\left(6\right)$

$=H\left(k_1\right)+H^{\′}\left(k_1\right)W_N^{k_0}{W}_M^{k_1}$

其中：

$\left\{\begin{array}{c}H\left(k_1\right)=\underset{h=0}{\overset{M/2-1}{\mathrm{\Σ}}}{X}_S\left(2h\right)W_N^{{2\mathrm{hk}}_0}{W}_{M/2}^{{\mathrm{hk}}_1}\\ H^{\′}\left(k_1\right)=\underset{h=0}{\overset{M/2-1}{\mathrm{\Σ}}}{X}_S(2h+1)W_N^{{2\mathrm{hk}}_0}{W}_{M/2}^{{\mathrm{hk}}_1}\end{array}\right.---\left(7\right)$

由以上推导可以看出：

当k₁＝0,1...M/2-1时：

$X(k_1,k_0)=H\left(k_1\right)+H^{\′}\left(k_1\right)W_N^{k_0}{W}_M^{k_1}---\left(8\right)$

当k₁＝M2+0,...,M/2+u,…M-1时：

$X(\frac{M}{2}+u,k_0)=H\left(u\right)-H^{\′}\left(u\right)W_N^{k_0}{W}_M^u---\left(9\right)$

其中，u＝0,1,2...M/2-1。

通过公式（7）可以求出序列H(k₁)和H'(k₁)在区间(0~M/2-1)内的各个值，进而利用公式（8）和（9）得到所有的X(k₁,k₀)。对于序列H(k₁)和H'(k₁)的计算可根据前面推导过程继续分解，直到分解为两个数据参与运算为止，该推导过程与基二FFT变换相同。

由公式（8）和（9）可以看出，第b级旋转因子为

此时，u的取值为u＝0,1，2...M/2-1；

继续推导，由公式（7）可知，令k₁=u，代入可知

$H\left(\mathrm{\μ}\right)=\underset{h=0}{\overset{M/2-1}{\mathrm{\Σ}}}{X}_S\left(h\right)W_{N/2}^{{\mathrm{hk}}_0}{W}_{M/2}^{{\mathrm{hk}}_1}$

将变换序列按照奇偶分解

$H\left(\mathrm{\μ}\right)=\underset{l=0}{\overset{M/4-1}{\mathrm{\Σ}}}{X}_S\left(2l\right)W_{N/2}^{{2\mathrm{lk}}_0}{W}_{M/2}^{{2\mathrm{lk}}_1}$

$+\underset{l=0}{\overset{M/4-1}{\mathrm{\Σ}}}{X}_S(2l+1)W_{N/2}^{(2l+1)k_0}{W}_{M/2}^{(2l+1)k_1}$

$=H\left(\mathrm{\η}\right)+H^{\′}\left(\mathrm{\η}\right)W_{N/2}^{k_0}{W}_{M/2}^{k_1}$

其中：

$\left\{\begin{array}{c}H\left(\mathrm{\η}\right)=\underset{l=0}{\overset{M/4-1}{\mathrm{\Σ}}}{X}_S\left(2l\right)W_{N/2}^{{2\mathrm{lk}}_0}{W}_{M/2}^{{2\mathrm{hk}}_1}\\ H^{\′}\left(\mathrm{\η}\right)=\underset{l=0}{\overset{M/4-1}{\mathrm{\Σ}}}{X}_S(2l+1)W_{N/2}^{{2\mathrm{lk}}_0}{W}_{M/2}^{{2\mathrm{lk}}_1}\end{array}\right.$

由此可以看出，第（b-1）级旋转因子为

此时u的取值为（0～M/4-1）；

以此类推，直到两个数据参与蝶形运算，设第b级蝶形计算所乘的旋转因子，那么该旋转因子可表示为：

$W(b,u)=W_{P\left(b\right)}^{k_0}{W}_{Q\left(b\right)}^u---\left(10\right)$

其中，u＝0,1,...,Q(b)-1，Q(b)＝M/2^c-b，P(b)＝N/2^c-b，c＝log₂(M)。

W(b,u)表示第b级蝶形计算所乘的旋转因子；

b表示蝶形运算的级数；

c表示FFT算法中所包含蝶形运算的总级数；

u表示b级蝶形运算时输出序列的序号；第b级蝶形运算时u会取遍u＝0,1,...,Q(b)-1。

k₀为当前进行行FFT变换的行数据的行序号，即FFT变换后得到的频域二维矩阵的列序号。

可见，通过重构第二次FFT中蝶形运算所乘的旋转因子W(b,u)，从而将铰链因子

隐藏在新的蝶形运算中，从而解决了引入额外运算的问题。

（二）解决三次显性转置的方案设计。

为了解决显性转置对处理效率的影响，可以采用变换行列号方式实现对二维矩阵的访问。具体来说，第一步，将一维序列x(n)拆分为L×M的二维矩阵，在不进行转置的情况下，处理器按列从L×M二维矩阵中读取列数据，通过Cache的缓存放入内存的指定空间，再指定空间读取列数据并进行列FFT，并将结果原位存回L×M二维矩阵。第二步，处理器按照行方向进行L个M点的行FFT，行FFT运算中各级蝶形运算所用的旋转因子由公式（10）确定，Cache的使用方法与第一步相同。

从以上处理可以看出，并没有进行显性转置操作。但是其带来的问题是，读取一列数据时，需要每读一个数据均进行行切换，而行切换需要切换时间。读取一行数据时，行数据不能太长，如果超出Cache的容量，则会出现因无法从Cache命中数据带来的各种问题。基于此，对二维矩阵的访问优化为如下方案：

对于列处理：

为了尽量少的行切换，可以在第一步中一次读取几列数据，读取列数的多少以及每列含数据量大小均需要根据Cache的容纳量计算，以实现对Cache的充分利用。为了实现高效并行处理和减少Cache丢失，必须满足两个条件：

（1）一次读取i列数据，如图2所示。考虑到处理器数据宽度的限制，每次可以处理w列数据，则读取的数据列数i需要是w的整数倍，本实施例中，设w=4，则i为4的整数倍。处理器一次的处理量w和该处理器总线带宽有关系，根据处理器每次可以读取数据的长度来决定w的值。

（2）总和考虑i和列长度L，使得i×L小于CacheLength，但尽可能贴近CacheLength，那么CacheLength为Cache的容量。这样，可以保证每打开一个页面（page）尽可能多读数据到缓存，以减少开/关页面引起的时延。

图2右图示出了从二维矩阵中读取i列数据的示意图，图2左图示出了i列数据在Cache中的分布。

对于行处理：

在第二步中每次读取一行进行FFT时，为了充分利用Cache，且行方向的点数在划分时也需要考虑Cache的容纳量。只有当读一行的数据尽可能填充Cache空间才能充分利用Cache的优势，协助数据访问。图3右图示出了从二维矩阵中读取一行数据的示意图，图3左图示出了该行数据在Cache中的分布。

可见，由于经过拆分后列方向数据访问跳变步长较小，处理数据缓存在Cache中，因此通过行列号变换方式访问矩阵，仅在每列处理的第一级和最后一级存在Cache丢失，后续处理中因数据位于缓存中而不会有额外开销。同时采取变换行列号的方法属于原位操作，可省去额外的转置空间，故采用该方法访问二维矩阵，要优于直接使用三次显性转置。然而，直接通过行列号进行访问矩阵数据时仍会因行列拆分不同引起Cache频繁丢失，这时就需要限制行列的拆分规则来发挥行列读取对Cache的利用率。

根据以上分析，本发明改进型FFT算法在TS201上实现的具体流程包括以下步骤，参见图4：

步骤一、将待处理的一维序列x(n)拆分为L段，存储L×M的二维矩阵，L为列的长度，M为行的长度；设定后续步骤中每次读取i列数据，i为正整数，且i为w的倍数，w为处理器每次处理数据的列数值；在保证

的基础上，使行长度M小于或等于TS201Cache大小，以M＝CacheLength为最佳。

步骤二、处理器按照列方向进行列FFT处理。

在处理过程中，处理器每次从L×M二维矩阵中读取i列的数据，通过Cache的缓存放入内存的指定空间，该指定空间是为计算中间量分配的存储空间；然后，处理器从指定空间读取数据并进行列FFT，并将结果原位存回L×M二维矩阵。由于（i×L）小于CacheLength，因此i×L点的数据均存储在Cache中，处理器每次从指定空间读取所需数据时，都能从Cache中命中。

由于处理器数据总线宽度有限，因此处理器在每次处理i列数据时，分多次处理，假设每次处理4列，则需处理

次。处理器共读取

次数据并进行FFT，从而实现M次L点的列FFT；

步骤三、对步骤二的结果，处理器按照行方向进行行FFT处理。

在处理过程中，处理器每次从L×M二维矩阵中读取一行数据，通过Cache的缓存放入内存的指定空间；然后处理器从指定空间读取缓存数据并进行行FFT，并将结果按照列方向顺序输出到内存中用于存储FFT变换结果的位置。处理器共读取L次数据并进行行FFT，从而实现L次M点的行FFT。同理，由于M小于或等于CacheLength，因此M点数据均存储在Cache中，处理器每次从指定空间读取所需数据时，都能从Cache中命中。

并且，本步骤的行FFT所用的旋转因子由公式（10）决定。

至此，本流程结束。

对不同浮点数FFT利用现有方法和本实施例中的方法得到的执行时间对比曲线，如图5所示；相比于传统的Winograd算法的执行效率提高了至少30%，比SingLeton算法的执行效率也提高了近15%。

综上所述，以上仅为本发明的较佳实施例而已，并非用于限定本发明的保护范围。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (1)

1.一种大点数FFT在处理器上的实现方法，其特征在于，包括：

步骤一、将待处理的一维序列x(n)分L段存储为L×M的二维矩阵，L为列的长度，M为行的长度；设定后续步骤二中每次读取i列数据，i为正整数，则在保证

的基础上，使行长度M小于或等于CacheLength；CacheLength为处理器所用Cache的容量；

步骤二、处理器进行列FFT；

处理器每次从L×M二维矩阵中读取i列数据，通过Cache的缓存放入内存的指定空间，然后从指定空间读取数据并进行列FFT，并将结果原位存回L×M二维矩阵；假设根据处理器数据宽度的限制，处理器每次处理数据量为w列，则i的取值为w的整数倍；读取的i列数据分次处理；处理器共读取

次数据并进行列FFT，从而实现M次L点的列FFT；

步骤三、处理器进行行FFT；

处理器每次从步骤二处理后的L×M二维矩阵中读取一行数据，通过Cache的缓存放入内存的指定空间，然后从指定空间读取缓存数据并进行行FFT，并将结果按照列方向输出；处理器共读取L次数据并进行行FFT，从而实现L次M点的行FFT；

本步骤的行FFT运算时第b级蝶形运算所用的旋转因子W(b,u)由下式确定：

$W(b,u)=W_{P\left(b\right)}^{k_0}\·W_{Q\left(b\right)}^u$

其中，W_n=e^-j2π/n，P(b)＝N/2^c-b，Q(b)＝M/2^c-b；

b表示FFT算法中蝶形运算的当前级数；

c表示FFT算法中所包含蝶形运算的总级数，c＝log₂(M)；

u表示b级蝶形运算输出序列的序号，取值范围是u＝0,1,...,Q(b)-1；

k₀为当前进行FFT变换的行数据的行序号。

Images