CN104615582A - 面向gpdsp的大点数一维fft向量化计算的方法 - Google Patents

Abstract

本发明公开了一种面向GPDSP的大点数一维FFT向量化计算的方法，在计算D＝2^d点一维FFT时，将d级FFT蝶形单元计算分两阶段完成：阶段I：前(d-m)级FFT蝶形单元计算的每一级FFT蝶形单元由DSP核的所有向量处理阵列计算单元按一维FFT蝶形单元计算方式按照向量化计算完成；直到2^m点序列数据能够全部存放在GPDSP的片内共享存储阵列中；阶段II：DSP核的所有向量处理阵列计算单元依次计算2^d-m次2^m点FFT计算；由DSP核的向量处理阵列采用一维转二维的计算方法，分割为更小点数的FFT计算，并由DSP核的向量处理阵列分别采用并行化和向量化计算方法计算完成。本发明能够显著提高FFT计算效率、减少数据传输时间开销。

Description

面向GPDSP的大点数一维FFT向量化计算的方法

技术领域

本发明主要涉及通用计算数字信号处理器(General-Purpose Digital Signal Processor，简称GPDSP)，特指一种适用于GPDSP的大点数一维FFT向量化计算的方法。

背景技术

离散傅里叶变换(Discrete Fourier Transform，DFT)在现代信号处理系统领域里应用广泛，如雷达信号处理、SAR图像处理、声纳计算、视频图像算法、频谱分析、语音识别等。傅里叶变换计算是典型的计算密集和访存密集型应用，例如N点的DFT变换的计算复杂度为o(N²)。1965年Cooley和Turkey提出一种快速傅立叶变换(Fast Fourier Transform，FFT)计算方法，可显著地减少运算量，计算复杂度由原来的o(N²)降到o(Nlog₂N)。信号处理应用通常对计算的实时性要求很高，FFT计算效率越高，信号处理的实时性就越好。

为提高FFT的计算性能，许多文献提出了不同的加速FFT计算的方法。专利申请号：201210218588.9的文献提供一种基于多核DSP平台的FFT并行方法，将需要FFT变换的原始数据均匀分配到所有处理器上并行处理。专利申请号：201010607219.X的文献提供一种通用DSP处理器中FFT计算实现装置和方法。专利申请号：200910179924.1的文献提供一种实现FFT及IFFT运算的装置和方法。专利申请号：201110163600.6的文献提供一种基于并行处理的FFT装置及其方法。然而，这些文献都没有针对大点数FFT的计算提供高效的计算方法。专利申请号：201110337733.0的文献提供了一种在向量处理器上基于SIMD实现FFT并行计算的方法，但是数据存放在片内向量阵列存储器中，不适合大点数(处理数据超过片内向量阵列存储器容量)FFT计算。专利申请号：201210448784.5的文献提供了一种大点数FFT的实现方法，是一种硬件实现方法，硬件开销大，不灵活。专利申请号：201310034812.3的文献提供了一种大点数FFT在处理器上的实现方法，其主要考虑Cache对执行效率的影响，行、列的划分也是与Cache行的长度相关。这两种大点数FFT计算方法不适合GPDSP的非Cache的向量阵列存储访存模式和向量处理阵列并发向量处理的体系结构特征。

在专利申请号为201310725118.6的文献(处于实审阶段)中提供了一种通用计算数字信号处理器(General-Purpose Digital Signal Processor，简称GPDSP)，它包含CPU核单元和DSP核单元，CPU核单元主要用于负责包括存储管理、文件控制、进程调度、中断管理任务在内的通用事务管理以及提供对通用操作系统的完整支持；DSP核单元包含若干强大计算能力的64位向量处理阵列，用于支持高密集运算任务的解算。

对于面向GPDSP的大点数FFT计算，由于DSP核的片内向量阵列存储器的容量不够大，需要进行FFT变换处理的原始序列数据不能完全存储在DSP核的片内向量阵列存储器上，通常存储在容量更大，但是访存速度更慢的片外DDR存储器中。因此，完成FFT计算的总执行时间主要包括：(1)数据在片外DDR存储器与DSP核的片内向量阵列存储器之间的数据传输时间；(2)DSP核的计算时间。在数据传输带宽有限的高性能处理系统上处理大点数FFT计算时，上述(1)所占的时间比重甚至远远超过上述(2)所占时间比重。因此，对于给定的GPDSP和点数的FFT计算，不同的FFT计算方法，涉及不同的数据搬移策略和计算流程，导致计算性能差异很大，对于面向GPDSP的大点数FFT计算而言，减少数据传输时间开销尤为重要。

发明内容

本发明要解决的技术问题就在于：针对现有技术存在的技术问题，本发明提供一种能够显著提高FFT计算效率、减少数据传输时间开销的高效大点数FFT向量化计算方法。

为解决上述技术问题，本发明采用以下技术方案：

一种面向GPDSP的大点数一维FFT向量化计算的方法，在计算D＝2^d点一维FFT时，将d级FFT蝶形单元计算分两阶段完成：

阶段I：前(d-m)级FFT蝶形单元计算的每一级FFT蝶形单元由DSP核的所有向量处理阵列计算单元按一维FFT蝶形单元计算方式按照向量化计算完成；直到2^m点序列数据能够全部存放在GPDSP的片内共享存储阵列中；

阶段II：DSP核的所有向量处理阵列计算单元依次计算2^d-m次2^m点FFT计算；由DSP核的向量处理阵列采用一维转二维的计算方法，分割为更小点数的FFT计算，并由DSP核的向量处理阵列分别采用并行化和向量化计算方法计算完成。

作为本发明的进一步改进：在所述阶段II中，将序列x(n)分组为N₁个长度为N₂的子序列，n＝0...N-1，将原来的N点一维FFT计算分三个子阶段完成：

(1)由DSP核的各个向量处理阵列计算单元并行的按列计算N₂个N₁点FFT计算；

(2)将计算结果传输到片外DDR存储器之前，将计算的结果与一个系数矩阵相乘；

(3)由DSP核的所有向量处理阵列计算单元协作按行计算N₁个N₂点FFT计算。

作为本发明的进一步改进：在所述阶段I中，先根据DSP核的向量数据加载能力和FFT处理数据的类型，确定DSP核的向量处理阵列每次处理的蝶形单元数为u个，片内向量阵列存储器能够存储的最大蝶形单元数量为y，y取值为u的整数倍，其中蝶形单元数据含序列数据和旋转因子；点数为D＝2^d的一维FFT的每一级蝶形单元数为D/2个，DSP核启动DMA从片外DDR存储器传输D/(2y)次数据到片内向量阵列存储器，每次传输y个蝶形单元数据量，传输的y个蝶形单元数据由DSP核分y/u次向量阵列处理；每次计算完的结果由DMA原位存回片外DDR存储器，最终完成D＝2^d的一维FFT的前(d-m)级FFT蝶形单元计算。

作为本发明的进一步改进：所述参数u的确定方法是：设DSP核的向量LOAD/STORE指令加载向量数据最大为p*w字节，FFT处理的序列数据的一个数据大小为z字节，则DSP核的向量处理阵列每次处理的蝶形单元数为p*w/z个。

作为本发明的进一步改进：所述DSP核的向量处理阵列在处理D＝2^d的一维FFT的每一级FFT蝶形单元计算时，根据DSP核的体系结构特征采用如下的双缓冲的乒乓方式进行：

片内向量阵列存储器能够存储的最大蝶形单元数量为y，在片内向量阵列存储器设立两个缓冲区，缓冲区的大小为y/2所需存储量，y/2取值为u的整数倍；采用双缓冲的方式进行蝶形单元的计算，即在一个缓冲区中依次计算y/2个蝶形单元计算的同时，通过DMA将上一次缓冲区的计算结果传输至片外DDR存储器中，以及将下一次缓冲区计算所需要的序列数据和旋转因子数据传输至片内向量阵列存储器中；直到完成该级全部蝶形单元的计算。

作为本发明的进一步改进：在所述阶段II中，由GPDSP的DSP核进行N＝2^m点的一维FFT计算的具体流程为：

(a)：根据DSP核的向量数据加载能力和FFT处理数据的类型，确定同一计算单元可同时计算t个N₁点FFT；根据DSP核的向量处理阵列计算单元数量和片内向量阵列存储器容量特征，将序列x(n)(n＝0,...,N-1)分组为N₁个长度为N₂的子序列；

(b)：按计算要求构建一个行数为N₁，列数为pt的N₁*(pt)规模的系数表和一个长度为N₁的列向量表；由GPDSP的DSP核按列进行N₂个N₁点FFT计算，计算结果与一个系数矩阵相乘；其中，N₁点FFT的计算采用原位计算，计算所需的旋转因子为N₁点的，每一个N₁点FFT的计算在同一个计算单元上完成，同一计算单元同时计算t个N₁点FFT，DSP核的向量处理阵列同时并行计算pt个N₁点FFT。pt个N₁点FFT的计算结果与系数矩阵的对应部分相乘；

(c)：由DSP核的向量处理阵列计算单元并行的按列计算N₂个N₁点FFT；DSP核依次并行计算pt个N₁点FFT，同一计算单元同时计算t个N₁点FFT；pt个N₁点FFT的计算结果与上述系数表的对应元素进行相乘操作，更新系数表；直到完成全部N₁点FFT的计算。

(d)：由DSP核的向量处理阵列计算单元协作按行计算N₁个N₂点FFT，每一个N₂点FFT的计算由DSP核的所有向量处理阵列计算单元协作完成，直到完成全部N₁点FFT的计算。

作为本发明的进一步改进：所述步骤(a)中N＝N₁N₂，N₁＝2ⁿ¹和N₂＝2ⁿ²，并且N₁和N₂的参数根据如下三个条件优选：(1)、1份或以上的N₂点FFT计算的序列数据和1份N₂点的旋转因子数据能够存放在DSP核的片内向量阵列存储器中；(2)、在满足条件(1)时，N₂点FFT的计算能够充分地发挥出DSP核的计算性能；(3)、N₁和N₂相等或尽量接近。

作为本发明的进一步改进：所述参数t的确定方法是：根据DSP核的向量数据加载能力和FFT处理数据的类型，确定同一计算单元可同时完成多少个N₁点FFT的计算；设DSP核的向量LOAD/STORE指令加载向量数据最大为p*w字节，对应到每一个计算单元最大为w字节，FFT处理的序列数据的一个数据大小为z字节，则同一计算单元可同时完成w/z个N₁点FFT的计算。

作为本发明的进一步改进：所述步骤(b)中，构建N₁*(pt)规模的系数表和长度为N₁的列向量表的方法是：设对任意的下标号i,k(0≤i<N₁,0≤k<pt)，其对应的系数表矩阵元素为b_ik,下标号i对应的列向量表元素为c_i；令s为i的位反序值，则系数表矩阵元素为 列向量表元素为 计算完pt个N₁点FFT以后，更新系数表，更新的方法是：对任意的下标号i,k(0≤i<N₁,0≤k<pt)，其对应的系数表矩阵元素b_ik更新为b_ik＝b_ik*c_i。

作为本发明的进一步改进：所述步骤(b)中按列计算N₂个N₁点FFT和步骤(c)中按行计算N₁个N₂点FFT时，根据DSP核的体系结构特征采用双缓冲的乒乓方式进行。

与现有技术相比，本发明的优点在于：

1、本发明的面向GPDSP的大点数一维FFT向量化计算的方法，原理简单，操作方便、能充分利用GPDSP中DSP核向量处理阵列的强大并行计算和高带宽向量数据加载能力，并显著降低片内与片外数据传输开销的大点数一维FFT向量化计算的方法，尤其适用于基于GPDSP结构实现大点数一维FFT的计算。

2、本发明的面向GPDSP的大点数一维FFT向量化计算的方法，为一种高效的向量化计算方法，尤其适合在高性能的GPDSP上计算大点数FFT，将数倍的提高计算性能。相比通常的一维FFT计算方法，本发明的技术方案增加了与系数矩阵相乘的计算量，但是却大幅度减少了片内向量阵列存储器和片外DDR存储器之间的数据传输量，因此将显著地降低总的FFT计算执行时间，尤其是在高性能的GPDSP上计算FFT，数据传输时间所占比重时较大，将数倍的提高计算性能。

附图说明

图1是本发明面向的GPDSP计算系统的简化存储模型示意图。

图2是本发明的流程示意图。

图3是本发明在具体应用实例中将一维序列分组为二维子序列的具体实施例流程示意图。

图4是本发明在具体应用实例中生成系数表和列向量表以及更新系数表的具体流程示意图。

图5是本发明在具体应用实例中按列进行子序列FFT计算的流程示意图。

图6是本发明在具体应用实例中按行进行子序列FFT计算的流程示意图。

图7是本发明在具体应用实例中利用双缓冲区进行子序列FFT计算的流程示意图。

具体实施方式

以下将结合说明书附图和具体实施例对本发明做进一步详细说明。

矩阵傅里叶算法(Matrix Fourier Algorithm，MFA)将大点数的一维FFT计算转化为多个小点数的FFT计算，其基本原理如下：

序列x(n)(n＝0，...，N-1)的离散傅立叶变换X(k)(k＝0，...，N-1)定义为：

$X\left(k\right)=\underset{n=0}{\overset{N-1}{\mathrm{\&Sigma;}}}x\left(n\right)W_N^{\mathrm{kn}}(k=0,...,N-1)$

其中 $W_N^{\mathrm{kn}}=e^{-j(2\mathrm{\&pi;}/N)\mathrm{kn}}(j=\sqrt{-1})$ 是旋转因子。

令N＝N₁xN₂，将序列x(n)分组为N₁个长度为N₂的子序列，即将一维序列x(n)转换为如下形式的二维数组序列：

$\left[\begin{array}{cccc}x\left(0\right)& x\left(1\right)& \&CenterDot;\&CenterDot;\&CenterDot;& x(N_2-1)\\ x\left(N_2\right)& x(N_2+1)& \&CenterDot;\&CenterDot;\&CenterDot;& x(2N_2-1)\\ \&CenterDot;& \&CenterDot;& & \&CenterDot;\\ \&CenterDot;& \&CenterDot;& \&CenterDot;\&CenterDot;\&CenterDot;& \&CenterDot;\\ \&CenterDot;& \&CenterDot;& & \&CenterDot;\\ x\left((N_1-1)N_2\right)& x((N_1-1)N_1+1)& \&CenterDot;\&CenterDot;\&CenterDot;& x(N_1{N}_2-1)\end{array}\right]$

令n和k的序号映射如下：

$\left\{\begin{array}{cc}{n=N}_2{n}_1+n_2,& \left\{\begin{array}{c}0\&le;n_1\&le;N_1-1\\ 0\&le;n_2\&le;N_2-1\end{array}\right.\\ k=k_1+N_1{k}_2,& \left\{\begin{array}{c}0\&le;k_1\&le;N_1-1\\ 0\&le;k_2\&le;N_2-1\end{array}\right.\end{array}\right.$

则X(k)可以进行如下变换：

$\begin{array}{c}X\left(k\right)=X(k_1+N_1{k}_2)=\underset{n_2=0}{\overset{N_2-1}{\mathrm{\&Sigma;}}}\underset{n_1=0}{\overset{N_1-1}{\mathrm{\&Sigma;}}}x(N_2{n}_1+n_2)W_N^{(k_1+N_1{k}_2)(N_2{n}_1+n_2)}\\ =\underset{n_2=0}{\overset{N_2-1}{\mathrm{\&Sigma;}}}\left\{\right[\underset{n_1=0}{\overset{N_1-1}{\mathrm{\&Sigma;}}}x(N_2{n}_1+n_2)W_{N_1}^{k_1{n}_1}\left]W_N^{k_1{n}_2}\right\}{W}_{N_2}^{k_2{n}_2}\end{array}---\left(1\right)$

从上式可以看出，计算N点一维DFT可以转化类似二维DFT的计算，即首先按列计算N₂个N₁点DFT，然后将计算的结果与一个系数矩阵相乘，接着按行计算N₁个N₂点DFT。

在通常的一维FFT计算方案中，对于每一级FFT计算，需要将N点的数据从片外DDR存储器传输到片内向量阵列存储器一次(本发明主要针对大点数FFT计算，假定原始数据只能存放在片外DDR存储器)，计算完以后还需要将计算结果从片内向量阵列存储器传输到片外DDR存储器；在下一级FFT计算时，重复上述过程，直到所有级数的FFT计算完成，如1M点FFT，包括20级FFT蝶形单元计算，则需要来回传输20次，数据传输时间开销很大。同时还要考虑N点的旋转因子数据的传输，因此数据传输开销很大，数据传输时间远大于计算时间，导致总的FFT计算效率很低。

本发明的面向GPDSP的大点数一维FFT向量化计算的方法，在计算D＝2^d点一维FFT时，将d级FFT蝶形单元计算分两阶段完成：

阶段I：前(d-m)级FFT蝶形单元计算的每一级FFT蝶形单元由DSP核的所有向量处理阵列计算单元按一维FFT蝶形单元计算方式按照向量化计算完成；直到2^m点序列数据能够全部存放在GPDSP的片内共享存储阵列中；

阶段II：DSP核的所有向量处理阵列计算单元依次计算2^d-m次2^m点FFT计算。

其中，上述阶段I当中由DSP核的向量处理阵列采用一维FFT蝶形单元计算方式按向量化计算方法完成；上述阶段II当中由DSP核的向量处理阵列采用一维转二维的计算方法，分割为更小点数的FFT计算，并由DSP核的向量处理阵列分别采用并行化和向量化计算方法计算完成。

更进一步，阶段II的计算当中采用如下优化计算方法：

将序列x(n)(n＝0...N-1)分组为N₁个长度为N₂的子序列，将原来的N点一维FFT计算分三阶段完成：

(1)由DSP核的各个向量处理阵列计算单元并行的按列计算N₂个N₁点FFT计算，这时，每次的N₁点FFT计算所需的数据和旋转因子总数据量比较小，只有原来数据存储量的N₂分之一，并且在该N₁点FFT计算中，每一级的FFT计算都是在片内向量阵列存储器中完成，数据不需要在片内向量阵列存储器和片外DDR存储器之间反复传输，直到该N₁点FFT计算完成，才将计算结果传输到片外DDR存储器中。

(2)将计算结果传输到片外DDR存储器之前，将计算的结果与一个系数矩阵相乘。本发明提供构建一个系数表和一个列向量表的巧妙方法，将原来与系数矩阵相乘的计算转化为与该规模更小的系数表的相乘计算，能够显著降低系数矩阵的数据存储量。后续计算所需要的系数表可以通过前面的系数表与该列向量表相乘计算得到。该方法所需要的存储量是原来的(pt+1)/N₂，对于大点数的FFT计算而言，通常(pt+1)远小于N₂，因此本方法能够节省大量存储空间。

(3)由DSP核的所有向量处理阵列计算单元协作按行计算N₁个N₂点FFT计算，这时，每次的N₂点FFT计算所需的数据和旋转因子总数据量比较小，只有原来数据存储量的N₁分之一，并且在该N₂点FFT计算中，每一级的FFT计算都是在片内向量阵列存储器中完成，数据不需要在片内向量阵列存储器和片外DDR存储器之间反复传输，直到该N₂点FFT计算完成，才将计算结果传输到片外DDR存储器中。

在一个具体应用实例中，如图1所示，为本发明在具体应用实例中所面向的GPDSP计算系统的简化存储模型示意图。GPDSP计算系统包括CPU核和DSP核，DSP核包含若干64位向量处理阵列计算单元，存储系统包括DSP核专用的片内向量阵列存储器，CPU核和DSP核共享的片内共享存储阵列、大容量的片外DDR存储器。

如图2所示，设GPDSP中DSP核的向量处理阵列计算单元数量为p个，片内向量阵列存储器容量为q字节，GPDSP的片内共享存储阵列容量为r字节，片内共享存储阵列能够存储的最大序列数据点数为N＝2^m，计算的一维FFT点数为D＝2^d。由于本发明面向大点数FFT计算，其中d，m均为大于10的整数。

本发明方法的详细流程为：

S1：判断d，m的大小，若d>m，则转步骤S2，否则，令m＝d，转步骤S4。

S2：由GPDSP的DSP核进行D＝2^d点的一维FFT的前(d-m)级FFT蝶形单元计算。

根据DSP核的向量数据加载能力和FFT处理数据的类型，确定DSP核的向量处理阵列每次处理的蝶形单元数为u个，片内向量阵列存储器能够存储的最大蝶形单元数量为y，y取值为u的整数倍，其中蝶形单元数据含序列数据和旋转因子。

点数为D＝2^d的一维FFT的每一级蝶形单元数为D/2个，DSP核需要启动DMA从片外DDR存储器传输D/(2y)次数据到片内向量阵列存储器，每次传输y个蝶形单元数据量(含序列数据和旋转因子)，传输的y个蝶形单元数据由DSP核分y/u次向量阵列处理；每次计算完的结果由DMA原位存回片外DDR存储器。最终完成D＝2^d的一维FFT的前(d-m)级FFT蝶形单元计算。

S3：由GPDSP的DSP核进行2^d-m个N＝2^m点的一维FFT计算。

在步骤S2当中完成D＝2^d点的一维FFT的前(d-m)级FFT蝶形单元计算后，由GPDSP的DSP核进行2^d-m个2^m点的一维FFT计算。N＝2^m点的一维FFT计算按步骤S 4执行。

S4：由GPDSP的DSP核进行N＝2^m点的一维FFT计算。

S4.1：根据DSP核的向量数据加载能力和FFT处理数据的类型，确定同一计算单元可同时计算t个N₁点FFT。

根据DSP核的向量处理阵列计算单元数量和片内向量阵列存储器容量特征，将序列x(n)(n＝0，...，N-1)分组为N₁个长度为N₂的子序列；

其中，N＝N₁N₂，N₁＝2ⁿ¹和N₂＝2ⁿ²，并且N₁和N₂的参数根据如下三个条件优选：(a)、1份或以上的N₂点FFT计算的序列数据和1份N₂点的旋转因子数据能够存放在DSP核的片内向量阵列存储器中；(b)、在满足条件(a)时，N₂点FFT的计算能够充分地发挥出DSP核的计算性能；(c)、N₁和N₂相等或尽量接近。

例如，DSP核的向量处理阵列计算单元数量为16个，片内向量阵列存储器容量为1M字节。对于1M点双精度浮点FFT计算而言，序列数据存储量(16MB)超过了片内向量阵列存储器容量。令N＝1024*1024，根据实际计算，1份1024点的双精度浮点FFT的序列数据和1份1024点的旋转因子数据能够存放在DSP核的片内向量阵列存储器中，并且在序列数据和旋转因子数据已经存放在片内向量阵列存储器中时，1024点的双精度浮点FFT能够充分地发挥出DSP核的计算性能。所以，可以选择N₁＝1024，N₂＝1024，N₁＝N₂。对于512K点双精度浮点FFT计算而言，N＝512*1024，可选择N₁＝512，N₂＝1024，N₁和N₂尽量接近。当然，这仅仅是本发明的优选参数方式之一，本领域的技术人员可以根据上述说明选择合适的N₁和N₂参数。

S4.2：按计算要求构建一个行数为N₁，列数为pt的N₁*(pt)规模的系数表和一个长度为N₁的列向量表；

由GPDSP的DSP核按列进行N₂个N₁点FFT计算，计算结果与一个系数矩阵相乘。其中，N₁点FFT的计算采用原位计算，计算所需的旋转因子为N₁点的，每一个N₁点FFT的计算在同一个计算单元上完成，同一计算单元同时计算t个N₁点FFT，DSP核的向量处理阵列同时并行计算pt个N₁点FFT。pt个N₁点FFT的计算结果与系数矩阵的对应部分相乘。

上述参数t的确定方法是：根据DSP核的向量数据加载能力和FFT处理数据的类型，确定同一计算单元可同时完成多少个N₁点FFT的计算。设DSP核的向量LOAD/STORE指令加载向量数据最大为p*w字节，对应到每一个计算单元最大为w字节，FFT处理的序列数据的一个数据大小为z字节，则同一计算单元可同时完成w/z个N₁点FFT的计算。

举例来说，DSP核的向量处理阵列计算单元数量为16个，向量LOAD/STORE指令加载的向量数据最大为16*16字节，对16位定点FFT计算，一个数据大小为4字节(实部和虚部为16位，各占2字节)，则同一计算单元可同时完成16/4＝4个N₁点FFT的计算；对单精度浮点FFT计算，一个数据大小为8字节(实部和虚部为单精度，各占4字节)，则同一计算单元可同时完成16/8＝2个N₁点FFT的计算；对双精度浮点FFT计算，一个数据大小为16字节(实部和虚部为双精度，各占8字节)，则同一计算单元可同时完成16/16＝1个N₁点FFT的计算。

根据前面的计算公式，计算结果与系数矩阵相乘是一个N₁*N₂的系数矩阵，矩阵下标号为k，n的元素为W_N ^kn(0≤k<N₁,0≤n<N₂)，若按照通常的系数矩阵计算，该矩阵所需要的存储容量为N＝N₁*N₂个数据量，计算不同列数的N₁点FFT需要传输该系数矩阵的不同列数。

S4.3：由DSP核的向量处理阵列计算单元并行的按列计算N₂个N₁点FFT。DSP核依次并行计算pt个N₁点FFT，同一计算单元同时计算t个N₁点FFT。pt个N₁点FFT的计算结果与上述系数表的对应元素进行相乘操作，更新系数表。直到完成全部N₁点FFT的计算。

S4.4：由DSP核的向量处理阵列计算单元协作按行计算N₁个N₂点FFT。每一个N₂点FFT的计算由DSP核的所有向量处理阵列计算单元协作完成。直到完成全部N₁点FFT的计算。

S5：判断DSP核是否处理完2^d-m个N＝2^m点的一维FFT计算，若没有，按照步骤S 4依次处理剩下的N＝2^m点的一维FFT计算，直到DSP核处理完全部2^d-m个N＝2^m点的一维FFT计算。

S6：完成D＝2^d点的一维FFT计算。

作为较佳的应用实例中，本实例对于上述步骤S2中参数u的确定方法是：根据DSP核的向量数据加载能力和FFT处理数据的类型，确定DSP核的向量处理阵列每次处理的蝶形单元数为多少。设DSP核的向量LOAD/STORE指令加载向量数据最大为p*w字节，FFT处理的序列数据的一个数据大小为z字节，则DSP核的向量处理阵列每次处理的蝶形单元数为p*w/z个。

举例来说，DSP核的向量处理阵列计算单元数量为16个，向量LOAD/STORE指令加载的向量数据最大为16*16字节，对16位定点FFT计算，一个数据大小为4字节(实部和虚部为16位，各占2字节)，则DSP核的向量处理阵列每次处理的蝶形单元数为16*16/4＝64个；对单精度浮点FFT计算，一个数据大小为8字节(实部和虚部为单精度，各占4字节)，则DSP核的向量处理阵列每次处理的蝶形单元数为16*16/8＝32个；对双精度浮点FFT计算，一个数据大小为16字节(实部和虚部为双精度，各占8字节)，则DSP核的向量处理阵列每次处理的蝶形单元数为16*16/16＝16个。

作为较佳的应用实例中，本实例在上述步骤S2中DSP核的向量处理阵列在处理D＝2^d的一维FFT的每一级FFT蝶形单元计算时，可以根据DSP核的体系结构特征采用如下的双缓冲的乒乓方式进行：

由步骤S2知道，片内向量阵列存储器能够存储的最大蝶形单元数量为y，则在片内向量阵列存储器设立两个缓冲区，缓冲区的大小为y/2所需存储量，y/2取值为u的整数倍。采用双缓冲的方式进行蝶形单元的计算，即在一个缓冲区中依次计算y/2个蝶形单元计算的同时，通过DMA将上一次缓冲区的计算结果传输至片外DDR存储器中，以及将下一次缓冲区计算所需要的序列数据和旋转因子数据传输至片内向量阵列存储器中。直到完成该级全部蝶形单元的计算。

作为较佳的应用实例中，本实例中基于系数矩阵所具有显著的规律，进一步提供一种能够显著降低系数矩阵数据存储量的改进方法。所述方法是：构建一个行数为N₁，列数为pt的N₁*(pt)规模的系数表和一个长度为N₁的列向量表，将原来与系数矩阵相乘的计算转化为与该系数表的相乘计算。后续计算所需要的系数表可以通过前面的系数表与该列向量表相乘计算得到。这种方法所需要的系数表和列向量表存储量只有原来的(pt+1)/N₂，对于大点数的FFT计算而言，通常(pt+1)远小于N₂，因此本方法能够节省大量存储空间，并且共享该数据空间。

构建N₁*(pt)规模的系数表和长度为N₁的列向量表的方法是：设对任意的下标号i，k(0≤i<N₁,0≤k<pt)，其对应的系数表矩阵元素为b_ik，下标号i对应的列向量表元素为c_i。令s为i的位反序值(即s的二进制码由i的二进制码反转得到)，则系数表矩阵元素为 列向量表元素为 $c_i=W_N^{\mathrm{pts}}(W_N^{\mathrm{pts}}=e^{-j(2\mathrm{\&pi;}/N)\mathrm{pts}}).$

计算完pt个N₁点FFT以后，需要更新系数表，更新的方法是：对任意的下标号i，k(0≤i<N₁,0≤k<pt)，其对应的系数表矩阵元素b_ik更新为b_ik＝b_ik*c_i。

作为较佳的应用实例中，本实例在步骤S4.2中按列计算N₂个N₁点FFT时，可以根据DSP核的体系结构特征采用如下的双缓冲的乒乓方式进行：

由步骤S4.2知道，每次DSP核同时并行计算pt个N₁点FFT，设pt个N₁点FFT的序列数据存储需求空间为s₁字节，N₁点旋转因子数据存储需求空间为s₂字节，系数表和列向量表存储需求空间为s₃字节。若片内向量阵列存储器容量q满足：q大于等于(2*v*s₁+s₂+s₃)，v为大于0的整数，则在片内向量阵列存储器设立两个缓冲区，缓冲区的大小为v*s₁。采用双缓冲的方式进行N₂个N₁点FFT的计算，即在一个缓冲区中依次计算vpt个N₁点FFT计算的同时，通过DMA将上一次缓冲区的计算结果传输至片外DDR存储器中，以及将下一次缓冲区计算所需要的序列数据传输至片内向量阵列存储器中，其中旋转因子、系数表和列向量表是数据共用的，只需要传输一次。若N₂不是vpt的整数倍，则最后若干个N₁点FFT的计算只用到部分计算单元。若片内向量阵列存储器容量q不满足：q大于等于(2*s₁+s₂+s₃)，则只设定一个数据缓冲区，依次在该缓冲区上传输数据和计算，直到完成全部N₁点FFT的计算。

作为较佳的应用实例中，本实例在步骤S4.3按行计算N₁个N₂点FFT时，可以根据DSP核的体系结构特征采用如下的双缓冲的乒乓方式进行：

由步骤S4.3知道，每一个N₂点FFT的计算由DSP核的所有向量处理阵列计算单元协作完成，设N₂点FFT的序列数据存储需求空间为u₁字节，N₂点旋转因子数据存储需求空间为u₂字节。若片内向量阵列存储器容量q满足：q大于等于(2*y*u₁+u₂)，y为大于0的整数，则在片内向量阵列存储器设立两个缓冲区，缓冲区的大小为y*u₁。采用双缓冲的方式进行N₁个N₂点FFT的计算，即在一个缓冲区中依次计算y个N₂点FFT计算的同时，通过DMA将上一次缓冲区的计算结果传输至片外DDR存储器中，以及将下一次缓冲区计算所需要的序列数据传输至片内向量阵列存储器中，其中旋转因子数据共用的，只需要传输一次。若片内向量阵列存储器容量q满足：q大于等于(2*u₁+u₂)，则只设定一个数据缓冲区，依次在该缓冲区上传输数据和计算，直到完成全部N₂点FFT的计算。

如图3所示，设计算1M点的FFT，N＝1024*1024＝1048576，原始一维序列数据x(n)(n＝0,1,2,,1024*1024-1)：x(0),x(1),x(2),,x(1024*1024-1)，令N₁＝N₂＝1024，N＝N₁N₂，则序列x(n)(n＝0…N-1)分组为1024个长度为1024的子序列：

$\left[\begin{array}{cccc}x\left(0\right)& x\left(1\right)& \&CenterDot;\&CenterDot;\&CenterDot;& x\left(1023\right)\\ x\left(1024\right)& x\left(1025\right)& \&CenterDot;\&CenterDot;\&CenterDot;& x\left(1027\right)\\ \&CenterDot;& \&CenterDot;& \&CenterDot;& \&CenterDot;\\ \&CenterDot;& \&CenterDot;& \&CenterDot;& \&CenterDot;\\ \&CenterDot;& \&CenterDot;& \&CenterDot;& \&CenterDot;\\ x(1023*1024)& x(1023*1024+1)& \&CenterDot;\&CenterDot;\&CenterDot;& x(1024*1024-1)\end{array}\right]$

如图4所示，设DSP核的向量处理阵列计算单元数量为4个，计算64点的双精度浮点FFT(选择64点仅仅用于说明生成系数表和列向量表的具体实施例流程)，N＝64，N₁＝N₂＝8。DSP核的向量LOAD/STORE指令加载向量数据最大为4*16字节，对双精度浮点FFT计算，一个数据大小为16字节(实部和虚部为双精度，各占8字节)，则同一计算单元可同时完成16/16＝1个8点FFT的计算。系数表的行数是8，列数是4*1＝4，列向量表的元素个数为8。经计算得知0，1，2，3，4，5，6，7的位反序值分别为0，4，2，6，1，5，3，7。因此，生成的系数表是：

$\left[\begin{array}{cccc}1& 1& 1& 1\\ 1& W_{64}^4& W_{64}^8& W_{64}^{12}\\ 1& W_{64}^2& W_{64}^4& W_{64}^6\\ 1& W_{64}^6& W_{64}^{12}& W_{64}^{18}\\ 1& W_{64}^1& W_{64}^2& W_{64}^3\\ 1& W_{64}^5& W_{64}^{10}& W_{64}^{15}\\ 1& W_{64}^3& W_{64}^6& W_{61}^9\\ 1& W_{64}^7& W_{64}^{14}& E_{64}^{21}\end{array}\right]$

生成的列向量表是：

$\left[\begin{array}{c}1\\ W_{64}^{16}\\ W_{64}^8\\ W_{64}^{24}\\ W_{64}^4\\ W_{64}^{20}\\ W_{64}^{12}\\ W_{64}^{28}\end{array}\right]$

该系数表的每一列向量与上述列向量表的对应元素相乘操作，得到更新的系数表是：

$\left[\begin{array}{cccc}1& 1& 1& 1\\ W_{64}^{16}& W_{64}^{20}& W_{64}^{24}& W_{64}^{28}\\ W_{64}^8& W_{64}^{10}& W_{64}^{12}& W_{64}^{14}\\ W_{64}^{24}& W_{64}^{30}& W_{64}^{36}& W_{64}^{42}\\ W_{64}^4& W_{64}^5& W_{64}^6& W_{64}^7\\ W_{64}^{20}& W_{64}^{25}& W_{64}^{30}& W_{64}^{35}\\ W_{64}^{12}& W_{64}^{15}& W_{64}^{18}& W_{64}^{21}\\ W_{64}^{28}& W_{64}^{34}& W_{64}^{42}& W_{64}^{49}\end{array}\right]$

如图5所示，是本发明的按列进行子序列FFT计算的流程示意图。设DSP核的向量处理阵列计算单元数量为4个，计算1M点的双精度浮点FFT，N＝1024*0124，N₁＝N₂＝1024。DSP核的向量LOAD/STORE指令加载数据最大为4*16字节，对双精度浮点FFT计算，一个数据大小为16字节(实部和虚部为双精度，各占8字节)，则同一计算单元可同时完成16/16＝1个1024点FFT的计算。DSP核同时进行4个1024点FFT的计算。

第一次计算时，DSP核将原始序列数据中的首4列序列数据和按要求计算的1024点旋转因子、系数表和列向量表数据传输到片内向量阵列存储器上，计算时，通过向量LOAD指令加载到向量寄存器，由DSP的4个计算单元并行的按列计算1024点FFT计算。其中第0个计算单元按原位计算第0列数据的FFT：x(0),x(1024),,x(1023*1024)，第1个计算单元按原位计算第1列数据的FFT：x(1),x(1025),,x(1023*1024+1)，第2个计算单元按原位计算第2列数据的FFT：x(2),x(1026),,x(1023*1024+2)，第3个计算单元按原位计算第3列数据的FFT：x(3),x(1027),,x(1023*1024+3)。

$\left[\begin{array}{cccc}x\left(0\right)& x\left(1\right)& x\left(2\right)& x\left(3\right)\\ x\left(1024\right)& x\left(1025\right)& x\left(1026\right)& x\left(1027\right)\\ \&CenterDot;& \&CenterDot;& \&CenterDot;& \&CenterDot;\\ \&CenterDot;& \&CenterDot;& \&CenterDot;& \&CenterDot;\\ \&CenterDot;& \&CenterDot;& \&CenterDot;& \&CenterDot;\\ x(1023*1024)& x(1023*1024+1)& x(1023*1024+2)\&CenterDot;& x(1024*1024+3)\end{array}\right]$

此时的系数表规模是1024*4的矩阵，列FFT的计算结果与该系数表的对应元素进行相乘操作，该操作完毕后更新系数表，将计算结果传出到片外DDR存储器相应位置。

依次循环下去，直到1024列的1024点的FFT全部计算完成。

如图6所示，是本本实施例中按行进行子序列FFT计算的流程示意图。

设DSP核的向量处理阵列计算单元数量为4个，计算1M点的双精度浮点FFT，N＝1024*0124，N₁＝N₂＝1024。

第一次计算时，DSP核将原始序列数据中的首行序列数据x(0),x(1),x(2),,x(1023)和按要求计算的1024点旋转因子传输到片内向量阵列存储器上，计算时，由DSP核的所有计算单元按原位计算协作完成该行的FFT计算。计算完毕后将计算结果传出到片外DDR存储器相应位置。

依次循环下去，直到1024行的1024点的FFT全部计算完成。

如图7所示，是本实施例中利用双缓冲区进行子序列FFT计算的流程示意图。这里以列子序列FFT计算为例进行说明。设每次DSP核同时并行计算pt个N₁点FFT，设pt个N₁点FFT的序列数据存储需求空间为s₁字节，旋转因子数据存储需求空间为s₂字节，小系数矩阵和更新列向量存储需求空间为s₃字节。并且片内向量阵列存储器容量q满足：q大于等于(2*v*s₁+s₂+s₃)，v为大于0的整数，则在片内向量阵列存储器设立两个缓冲区：缓冲区A和缓冲区B，缓冲区的大小为v*s₁。采用双缓冲的方式进行N₂个N₁点FFT的计算，即在一个缓冲区中依次计算vpt个N₁点FFT计算的同时，通过DMA将上一个缓冲区的计算结果传输至片外DDR存储器中，以及将下一个缓冲区计算所需要的序列数据传输至片内向量阵列存储器中，其中旋转因子、系数表和列向量表是数据共用的，只需要传输一次。依次在缓冲区上传输数据和计算，直到完成全部N₁点FFT的计算。

以上仅是本发明的优选实施方式，本发明的保护范围并不仅局限于上述实施例，凡属于本发明思路下的技术方案均属于本发明的保护范围。应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理前提下的若干改进和润饰，应视为本发明的保护范围。

Claims (10)

1.一种面向GPDSP的大点数一维FFT向量化计算的方法，其特征在于，在计算D＝2^d点一维FFT时，将d级FFT蝶形单元计算分两阶段完成：

阶段I：前(d-m)级FFT蝶形单元计算的每一级FFT蝶形单元由DSP核的所有向量处理阵列计算单元按一维FFT蝶形单元计算方式按照向量化计算完成；直到2^m点序列数据能够全部存放在GPDSP的片内共享存储阵列中；

阶段II：DSP核的所有向量处理阵列计算单元依次计算2^d-m次2^m点FFT计算；由DSP核的向量处理阵列采用一维转二维的计算方法，分割为更小点数的FFT计算，并由DSP核的向量处理阵列分别采用并行化和向量化计算方法计算完成。

2.根据权利要求1所述的面向GPDSP的大点数一维FFT向量化计算的方法，其特征在于，在所述阶段II中，将序列x(n)分组为N₁个长度为N₂的子序列，n＝0...N-1，将原来的N点一维FFT计算分三个子阶段完成：

(1)由DSP核的各个向量处理阵列计算单元并行的按列计算N₂个N₁点FFT计算；

(2)将计算结果传输到片外DDR存储器之前，将计算的结果与一个系数矩阵相乘；

(3)由DSP核的所有向量处理阵列计算单元协作按行计算N₁个N₂点FFT计算。

3.根据权利要求1或2所述的面向GPDSP的大点数一维FFT向量化计算的方法，其特征在于，在所述阶段I中，先根据DSP核的向量数据加载能力和FFT处理数据的类型，确定DSP核的向量处理阵列每次处理的蝶形单元数为u个，片内向量阵列存储器能够存储的最大蝶形单元数量为y，y取值为u的整数倍，其中蝶形单元数据含序列数据和旋转因子；点数为D＝2^d的一维FFT的每一级蝶形单元数为D/2个，DSP核启动DMA从片外DDR存储器传输D/(2y)次数据到片内向量阵列存储器，每次传输y个蝶形单元数据量，传输的y个蝶形单元数据由DSP核分y/u次向量阵列处理；每次计算完的结果由DMA原位存回片外DDR存储器，最终完成D＝2^d的一维FFT的前(d-m)级FFT蝶形单元计算。

4.根据权利要求3所述的面向GPDSP的大点数一维FFT向量化计算的方法，其特征在于，所述参数u的确定方法是：设DSP核的向量LOAD/STORE指令加载向量数据最大为p*w字节，FFT处理的序列数据的一个数据大小为z字节，则DSP核的向量处理阵列每次处理的蝶形单元数为p*w/z个。

5.根据权利要求3所述的面向GPDSP的大点数一维FFT向量化计算的方法，其特征在于，所述DSP核的向量处理阵列在处理D＝2^d的一维FFT的每一级FFT蝶形单元计算时，根据DSP核的体系结构特征采用如下的双缓冲的乒乓方式进行：

片内向量阵列存储器能够存储的最大蝶形单元数量为y，在片内向量阵列存储器设立两个缓冲区，缓冲区的大小为y/2所需存储量，y/2取值为u的整数倍；采用双缓冲的方式进行蝶形单元的计算，即在一个缓冲区中依次计算y/2个蝶形单元计算的同时，通过DMA将上一次缓冲区的计算结果传输至片外DDR存储器中，以及将下一次缓冲区计算所需要的序列数据和旋转因子数据传输至片内向量阵列存储器中；直到完成该级全部蝶形单元的计算。

6.根据权利要求1或2所述的面向GPDSP的大点数一维FFT向量化计算的方法，其特征在于，在所述阶段II中，由GPDSP的DSP核进行N＝2^m点的一维FFT计算的具体流程为：

(a)：根据DSP核的向量数据加载能力和FFT处理数据的类型，确定同一计算单元可同时计算t个N₁点FFT；根据DSP核的向量处理阵列计算单元数量和片内向量阵列存储器容量特征，将序列x(n)(n＝0,...,N-1)分组为N₁个长度为N₂的子序列；

(b)：按计算要求构建一个行数为N₁，列数为pt的N₁*(pt)规模的系数表和一个长度为N₁的列向量表；由GPDSP的DSP核按列进行N₂个N₁点FFT计算，计算结果与一个系数矩阵相乘；其中，N₁点FFT的计算采用原位计算，计算所需的旋转因子为N₁点的，每一个N₁点FFT的计算在同一个计算单元上完成，同一计算单元同时计算t个N₁点FFT，DSP核的向量处理阵列同时并行计算pt个N₁点FFT；pt个N₁点FFT的计算结果与系数矩阵的对应部分相乘；

(c)：由DSP核的向量处理阵列计算单元并行的按列计算N₂个N₁点FFT；DSP核依次并行计算pt个N₁点FFT，同一计算单元同时计算t个N₁点FFT；pt个N₁点FFT的计算结果与上述系数表的对应元素进行相乘操作，更新系数表；直到完成全部N₁点FFT的计算；

(d)：由DSP核的向量处理阵列计算单元协作按行计算N₁个N₂点FFT，每一个N₂点FFT的计算由DSP核的所有向量处理阵列计算单元协作完成，直到完成全部N₁点FFT的计算。

7.根据权利要求6所述的面向GPDSP的大点数一维FFT向量化计算的方法，其特征在于，所述步骤(a)中N＝N₁N₂，N₁＝2ⁿ¹和N₂＝2ⁿ²，并且N₁和N₂的参数根据如下三个条件优选：(1)、1份或以上的N₂点FFT计算的序列数据和1份N₂点的旋转因子数据能够存放在DSP核的片内向量阵列存储器中；(2)、在满足条件(1)时，N₂点FFT的计算能够充分地发挥出DSP核的计算性能；(3)、N₁和N₂相等或尽量接近。

8.根据权利要求6所述的面向GPDSP的大点数一维FFT向量化计算的方法，其特征在于，所述参数t的确定方法是：根据DSP核的向量数据加载能力和FFT处理数据的类型，确定同一计算单元可同时完成多少个N₁点FFT的计算；设DSP核的向量LOAD/STORE指令加载向量数据最大为p*w字节，对应到每一个计算单元最大为w字节，FFT处理的序列数据的一个数据大小为z字节，则同一计算单元可同时完成w/z个N₁点FFT的计算。

9.根据权利要求6所述的面向GPDSP的大点数一维FFT向量化计算的方法，其特征在于，所述步骤(b)中，构建N₁*(pt)规模的系数表和长度为N₁的列向量表的方法是：设对任意的下标号i,k(0≤i<N₁,0≤k<pt)，其对应的系数表矩阵元素为b_ik,下标号i对应的列向量表元素为c_i；令s为i的位反序值，则系数表矩阵元素为列向量表元素为计算完pt个N₁点FFT以后，更新系数表，更新的方法是：对任意的下标号i,k(0≤i<N₁,0≤k<pt)，其对应的系数表矩阵元素b_ik更新为b_ik＝b_ik*c_i。

10.根据权利要求6所述的面向GPDSP的大点数一维FFT向量化计算的方法，其特征在于，所述步骤(b)中按列计算N₂个N₁点FFT和步骤(c)中按行计算N₁个N₂点FFT时，根据DSP核的体系结构特征采用双缓冲的乒乓方式进行。