CN103106181B - 一种大点数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}{\Σ}}x\left(n\right)W_N^{nk},k,n=0,1,...,N-1---\left(1\right)$

其中，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=0,1,...,M-1;n_1=0,1,...,L-1\\ k_0=0,1,...,L-1;k_1=0,1,...,M-1\end{array}\right.---\left(2\right)$

将(2)式代入公式(1)中得：

$\begin{array}{c}X\left(k\right)=X(k_1,k_0)=\underset{n_0=0}{\overset{M-1}{\Σ}}\underset{n_1=0}{\overset{L-1}{\Σ}}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}{\Σ}}\left(\right\{\underset{n_1=0}{\overset{L-1}{\Σ}}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}\end{array}---\left(3\right)$

其中，与式(1)的形式类似，相当于对一点数据进行FFT，设 $X_s\left(n_0\right)=\underset{n_1=0}{\overset{L-1}{\Σ}}x(n_1,n_0)W_L^{n_1{k}_0},$ 则 $\underset{n_0=0}{\overset{M-1}{\Σ}}\left(X_s\right(n_0)\·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_P(b)＝e^-j2π/P(b)，W_Q(b)＝e^-j2π/Q(b)，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变换后得到的频域二维矩阵的列序号；

所述处理器为TS201处理器。

有益效果：

本发明提供了一种应改进型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}{\Σ}}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}{\Σ}}\left(X_s\right(n_0\left)W_N^{n_0{k}_0}{W}_M^{n_0{k}_1}\right)---\left(5\right)$

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

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

其中：

$\left\{\begin{array}{c}H\left(k_1\right)=\underset{h=0}{\overset{M/2-1}{\Σ}}{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}{\Σ}}{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₁＝M/2+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}{\Σ}}{X}_S\left(h\right)W_{N/2}^{{\mathrm{hk}}_0}{W}_{M/2}^{{\mathrm{hk}}_1}$

将变换序列按照奇偶分解

$\begin{array}{c}H\left(\mathrm{\μ}\right)=\underset{l=0}{\overset{M/4-1}{\Σ}}{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}{\Σ}}{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}\end{array}$

其中：

$\left\{\begin{array}{c}H\left(\mathrm{\η}\right)=\underset{l=0}{\overset{M/4-1}{\Σ}}{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}{\Σ}}{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在处理器上的实现方法，其特征在于，包括：

步骤一、设FFT变换前的序列为x(n)，将待处理的一维序列x(n)分L段存储为L×M的二维矩阵，n＝0，1，…，N-1，L为列的长度，M为行的长度，N为序列中元素总数；设定后续步骤二中每次读取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_P(b)＝e^-j2π/P(b)，W_Q(b)＝e^-j2π/Q(b)，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变换的行数据的行序号；

所述处理器为TS201处理器。