CN104142811B - 一种基于数字信号处理的高效并行处理优化方法 - Google Patents
一种基于数字信号处理的高效并行处理优化方法 Download PDFInfo
- Publication number
- CN104142811B CN104142811B CN201410341689.4A CN201410341689A CN104142811B CN 104142811 B CN104142811 B CN 104142811B CN 201410341689 A CN201410341689 A CN 201410341689A CN 104142811 B CN104142811 B CN 104142811B
- Authority
- CN
- China
- Prior art keywords
- sequence
- backward
- subsequence
- subsequences
- digital signal
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
Landscapes
- Compression Or Coding Systems Of Tv Signals (AREA)
- Complex Calculations (AREA)
Abstract
本发明涉及基于数字信号处理的高效并行处理优化方法，该方法包括下列顺序的步骤：（1）在原位逆序的基础上进行块原位部分逆序；（2）三阶/四阶合并；（3）中间的二阶合并循环；（4）最后两阶计算/最终逆序。本发明针对数字信号处理器的高数据并行处理能力，创新地提出基于块原位的部分逆序方法，在此基础上，构建高效的原位FFT算法，该FFT优化方法在保证时间效率的同时，大幅降低了对内存空间的需求，很好的满足数字信号处理领域对算法的空间受限的需求。
Description
技术领域
本发明涉及数字信号处理技术领域，尤其是一种基于数字信号处理的高效并行处理优化方法。
背景技术
数字信号处理器的内存一般不大，为了处理数据密集型的应用，必须在设计算法时考虑空间复杂度。FFT算法是数字信号处理领域的一个重要的算法，它分为时域抽取算法和频域抽取算法，但计算过程都涉及到逆序计算。传统的非FFT算法的空间复杂度为O(2*N)，单位为复数，其中N为FFT点数。而原位FFT的空间复杂度为O(N)，和非原位FFT相比，可以节约N个复数所占的内存空间，而数字信号处理器往往能够提供多个连续地址的同时访问能力。
为了能够利用这种强大的高数据并行的内存访问能力，降低FFT优化算法对数字信号处理器内存空间的需求，使得数字信号处理器支持更大点数FFT的处理，针对高数据处理能力的原位FFT算法就有了迫切的研究需要。
发明内容
本发明的目的在于提供一种既可以高效利用数字信号处理器提供的对多个连续地址的内存访问能力，又可以大幅降低算法对内存空间需求。一种基于数字信号处理的高效并行处理优化方法，该方法包括下列顺序的步骤：
(1)进行基于置换的原位逆序的除高/低part之外的部分逆序；
(2)三阶/四阶合并；
(3)中间的二阶合并循环；
(4)最后两阶计算/最终逆序；
所述基于置换的原位逆序分为以下迭代步骤：
步骤一：把N个复数的输入序列(A(0),A(1),A(2),…，A(N-1))看成一个步长因子为1组成的1个子序列，把前N/2序列的奇数序列和后N/2序列的偶数序列对应置换，即前半部分序列(A(0),A(1),A(2),…，A(N/2-1))的奇数部分A(1),A(3),A(5)，…，A(N/2-1)和后半部分序列A(N/2),A(N/2+1),A(N/2+2),…，A(N-1)的偶数部分A(N/2),A(N/2+2),…，A(N-2)依次对应置换；
步骤二：把步骤一置换后的序列看成是一个步长因子为2组成的4个序列,其中前两个子序列在前N/2位置上，两个子序列交迭分布,即这两个子序列为(A(0),A(2),…，A(N/2-2))和(A(1),A(3),…，A(N/2-1))；后两个子序列在后N/2位置上，两个序列交迭分布,即这两个子序列为(A(N/2),A(N/2+2),…，A(N-2))和(A(N/2+1),A(N/2+3),…，A(N-1))；把4个子序列各自的前半部分的奇数部分和后半部分的偶数部分依次对应置换；
步骤三：把步骤二置换后的序列看成是一个步长因子为4组成的16个子序列，其中0～3个子序列在位置0～(N/4-1)上，四个子序列交迭分布,即这四个子序列为(A(0),A(0+4),A(0+8)，…，A(N/4-4))，(A(1),A(1+4),A(1+8),…，A(N/4-3))，(A(2),A(2+4),A(2+8),…，A(N/4-2))，(A(3),A(3+4),A(3+8),…，A(N/4-1))；其中4～7个子序列在位置N/4～(N/2-1)上，四个子序列交迭分布,即这四个子序列为(A(N/4+0),A(N/4+0+4),A(N/4+0+8)，…，A(N/2-4))，(A(N/4+1),A(N/4+1+4),A(N/4+1+8),…，A(N/2-3))，(A(N/4+2),A(N/4+2+4),A(N/4+2+8),…，A(N/2-2))，(A(N/4+3),A(N/4+3+4),A(N/4+3+8),…，A(N/2-1))；其中8～11个子序列在位置N/2～(3*N/4-1)上，四个子序列交迭分布,即这4个子序列为(A(N/2+0),A(N/2+0+4),A(N/2+0+8)，…，A(3*N/4-4))，(A(N/2+1),A(N/2+1+4),A(N/2+1+8),…，A(3*N/4-3))，(A(N/2+2),A(N/2+2+4),A(N/2+2+8),…，A(3*N/4-2))，(A(N/2+3),A(N/2+3+4),A(N/2+3+8),…，A(3*N/4-1))；其中11～15个子序列在位置3*N/4～(N-1)上，四个子序列交迭分布,即这4个子序列为(A(3*N/4+0),A(3*N/4+0+4),A(3*N/4+0+8)，…，A(N-4))，(A(3*N/4+1),A(3*N/4+1+4),A(3*N/4+1+8),…，A(N-3))，(A(3*N/4+2),A(3*N/4+2+4),A(3*N/4+2+8),…，A(N-2))，(A(3*N/4+3),A(3*N/4+3+4),A(3*N/4+3+8),…，A(N-1))；把16个子序列各自的前半部分的奇数部分和后半部分的偶数部分依次对应置换；
步骤四：把步骤三置换后的序列按照步长因子是上一次迭代步长因子的2倍、子序列个数是上一次迭代子序列个数的4倍的方法划分出若干个子序列，依次对该若干个子序列进行各自的前半部分的奇数部分和后半部分偶数部分的进行置换；
步骤五：把整个序列的多个子序列每次所完成的置换作为一个迭代步骤，那么总迭代步骤为R/2次，该值向下取整，R为LogN。
所述原位部分逆序是指：首先把输入序列按照高part位连续划分成为几个大块；其次将低part位则作为一个整体小块参与块原位置换；最后，部分逆序转化为求解(R-2*part)个中间位的逆序操作，此时则可利用原位逆序方法完成输入序列的部分逆序过程；
在完成部分逆序后，对此时的序列首先进行三阶/四阶合并：若是奇数阶，则进行三阶合并，即从内存读取序列到数字信号处理器的寄存器进行前三阶运算之后，才写回内存；若是偶数阶，则进行四阶合并，即从内存读取序列到数字信号处理器的寄存器进行前四阶运算之后，才写回内存。
在完成三阶/四阶合并后，进行中间阶的两阶合并：从内存读取序列和旋转因子到数字信号处理器寄存器的寄存器进行两阶运算后，才写回内存。
在完成中间阶的两阶合并循环后，最后进行两阶计算/最终逆序：从内存读取序列和旋转因子到数字信号处理器寄存器的寄存器进行最后两阶运算后，写回内存时，完成最终的高低part位逆序。
由上述技术方案可知，本发明针对数字信号处理器要求的高数据并行处理能力，创新地提出基于块原位的部分逆序方法，在此基础上，构建高效的原位FFT算法，该FFT优化方法在保证时间效率的同时，大幅降低了对内存空间的需求，很好的满足数字信号处理领域对算法的空间受限的需求。
附图说明
图1为整块逆序示意图。
图2为数据块组合示意图。
图3为高/低2位不参与逆序的部分逆序示意图。
图4为部分逆序后的输入序列的位序状态示意图。
具体实施方式
基于数字信号处理的高效并行处理优化方法，包括：(1)进行基于置换的原位逆序的除高/低part之外的部分逆序；(2)三阶/四阶合并；(3)中间的二阶合并循环；(4)最后两阶计算/最终逆序。
所述基于置换的原位逆序分为以下迭代步骤：
步骤一：把N个复数的输入序列(A(0),A(1),A(2),…，A(N-1))看成一个步长因子为1组成的1个子序列，把前N/2序列的奇数序列和后N/2序列的偶数序列对应置换，即前半部分序列(A(0),A(1),A(2),…，A(N/2-1))的奇数部分A(1),A(3),A(5)，…，A(N/2-1)和后半部分序列A(N/2),A(N/2+1),A(N/2+2),…，A(N-1)的偶数部分A(N/2),A(N/2+2),…，A(N-2)依次对应置换；
步骤二：把步骤一置换后的序列看成是一个步长因子为2组成的4个序列,其中前两个子序列在前N/2位置上，两个子序列交迭分布,即这两个子序列为(A(0),A(2),…，A(N/2-2))和(A(1),A(3),…，A(N/2-1))；后两个子序列在后N/2位置上，两个序列交迭分布,即这两个子序列为(A(N/2),A(N/2+2),…，A(N-2))和(A(N/2+1),A(N/2+3),…，A(N-1))；把4个子序列各自的前半部分的奇数部分和后半部分的偶数部分依次对应置换；
步骤三：把步骤二置换后的序列看成是一个步长因子为4组成的16个子序列，其中0～3个子序列在位置0～(N/4-1)上，四个子序列交迭分布,即这四个子序列为(A(0),A(0+4),A(0+8)，…，A(N/4-4))，(A(1),A(1+4),A(1+8),…，A(N/4-3))，(A(2),A(2+4),A(2+8),…，A(N/4-2))，(A(3),A(3+4),A(3+8),…，A(N/4-1))；其中4～7个子序列在位置N/4～(N/2-1)上，四个子序列交迭分布,即这四个子序列为(A(N/4+0),A(N/4+0+4),A(N/4+0+8)，…，A(N/2-4))，(A(N/4+1),A(N/4+1+4),A(N/4+1+8),…，A(N/2-3))，(A(N/4+2),A(N/4+2+4),A(N/4+2+8),…，A(N/2-2))，(A(N/4+3),A(N/4+3+4),A(N/4+3+8),…，A(N/2-1))；其中8～11个子序列在位置N/2～(3*N/4-1)上，四个子序列交迭分布,即这4个子序列为(A(N/2+0),A(N/2+0+4),A(N/2+0+8)，…，A(3*N/4-4))，(A(N/2+1),A(N/2+1+4),A(N/2+1+8),…，A(3*N/4-3))，(A(N/2+2),A(N/2+2+4),A(N/2+2+8),…，A(3*N/4-2))，(A(N/2+3),A(N/2+3+4),A(N/2+3+8),…，A(3*N/4-1))；其中11～15个子序列在位置3*N/4～(N-1)上，四个子序列交迭分布,即这4个子序列为(A(3*N/4+0),A(3*N/4+0+4),A(3*N/4+0+8)，…，A(N-4))，(A(3*N/4+1),A(3*N/4+1+4),A(3*N/4+1+8),…，A(N-3))，(A(3*N/4+2),A(3*N/4+2+4),A(3*N/4+2+8),…，A(N-2))，(A(3*N/4+3),A(3*N/4+3+4),A(3*N/4+3+8),…，A(N-1))；把16个子序列各自的前半部分的奇数部分和后半部分的偶数部分依次对应置换；
步骤四：把步骤三置换后的序列按照步长因子是上一次迭代步长因子的2倍、子序列个数是上一次迭代子序列个数的4倍的方法划分出若干个子序列，依次对该若干个子序列进行各自的前半部分的奇数部分和后半部分偶数部分的进行置换；
步骤五：把整个序列的多个子序列每次所完成的置换作为一个迭代步骤，那么总迭代步骤为R/2次，该值向下取整，R为LogN。
所述原位部分逆序是指：首先把输入序列按照高part位连续划分成为几个大块；其次将低part位则作为一个整体小块参与块原位置换；最后，部分逆序转化为求解(R-2*part)个中间位的逆序操作，此时则可利用原位逆序方法完成输入序列的部分逆序过程。假定待逆序的序列A为(A(0),A(1),A(2),…，A(N-1))，其中A(i)始终表示为在i位置上的复数，在算法执行逆序和变换的过程中，所指并不是同一个复数，特化N等于1024，则阶数R＝10，即该序列为((A(0),A(1),A(2),…，A(1023)),该特化也并未使得该方法失去广泛性,其中R＝log N；
步骤一：位序高2位的划分；由于是部分逆序，并且高/低2位都不参与逆序，从高位看，(A(0),A(1),A(2),…，A(1023)把位置序号(0,1,2,…,1023)转换成10位二进制，则待逆序序列A可以看(A(00********),…，A(01********)，…，A(10********)，…，A(11********)，…)。即依次为位置序号高2位为00的序列(A(0),A(1),…，A(255))，位置序号高2位为01的序列(A(256),A(257),…，A(511))，位置序号高2位为10的序列(A(512),A(513),…，A(767))，位置序号高2位为11的序列(A(768),A(769),…，A(1023))。此时1024点待逆序的序列本质上被高两位等分成为了4个序列，分别记作B0,B1,B2,B3,序列大小都为256。由于部分逆序高/低二位不参与逆序，所以针对A序列的部分逆序就转换成为针对B0，B1，B2，B3的4个序列的分别部分逆序。由于此时不考虑高2位的话B0，B1，B2,B3四个序列的的低8位位序序号序列都为(0,1,2,…,255)。所以针对B0,B1，B2，B3的部分逆序方法是等价的。为了简化叙述，接下来的步骤只叙述B0的部分逆序方法；
步骤二：序列B0，高2位是00，针对低2位的划分。B0序列(A(0),A(1),…，A(255))，如果把每个连续4个点看成一个小块，则B0序列可以看成是64个小块组成的序列(B(0),B(1),B(2)，…，B(63)),其中B(0)表示由(A(0),A(1),A(2),A(3))组成的小块，B(1)表示由(A(4),A(5),A(6),A(7))组成的小块,B(2)表示由(A(8),A(9),A(10),A(11))组成的小块,…，B(63)表示由(A(252),A(253),A(254),A(255))组成的小块。这种划分方法本质上是考虑低2位的变换，每个小块B(i)的四个点低2位的变化规律都是(00，01，10，11)，以小块为单位，进行中间6位序列的逆序，则可知，高/低2位都不会改变，而针对B0序列(B(0),B(1),B(2)，…，B(63))的6位原位逆序置换，就完成了B0序列(A(0),A(1),…，A(255))的高/低2位不变的部分逆序置换。此时问题转化为针对B0序列(B(0),B(1),B(2)，…，B(63))的6位原位逆序置换，置换的基本单位是4个复数组成的小块,此时B(i)表示在i位置(i取值范围为0～63)的小块，所表示的实质内容随着置换的进行是变化的；
步骤三：把64个小块组成的序列(B(0),B(1),B(2),…，B(63))看成一个步长因子为1组成的1个子序列块，把该子序列块把前半个序列块(B(0),B(1),B(2),…，B(31))的奇数序列块B(1),B(3),B(5),…，B(31)和后半个序列块(B(32),B(33),B(34),…，B(63))的偶数序列块B(32),B(34),B(36),…，B(62)依次进行对应置换；
步骤四：把当前序列块看成是步长因子为2组成的4个子序列块,其中前2个序列在0～31序列块上，两个序列块交迭分布,即这两个序列为(B(0),B(2),…，B(30))和(B(1),B(3),…，B(31))；后2个序列块在32～63序列块上，两个序列交迭分布,即这两个序列为(B(32),B(34),…，B(62))和(B(33),A(35),…，B(63))；依次对4个子序列的各自的前半个序列块的奇数序列块和后半个序列块的偶数序列块进行置换。例如第1个序列块(B(0),B(2),…，B(30)),前半个序列块的奇数序列块是B(2),B(6),B(10),B(14),后半个序列块的偶数序列块是B(16),B(20),B(24),B(28),则依次进行B(2)和B(16),B(6)和B(20),B(10)和B(24),B(14)和B(28)的对应置换。其它三个序列块类似；
步骤五：把当前序列块看成是一个步长因子为4组成的16个子序列，其中0～3个子序列块分布在当前序列块第0～15个块上，即B(0),B(1),B(2),…，B(15)，4个子序列块交迭分布,即4个子序列块分别为(B(0),B(4),B(8)，B(12))，(B(1),B(5),B(9),B(13))，(B(2),B(6),B(10),B(14)),(B(3),B(7),B(11),B(15))；4～7个子序列块分布在当前序列块第16～31个块上，即B(16),B(17),B(18),…，B(31)，4个子序列块交迭分布，与0～4个子序列块的组成规律一致；8～11个子序列块分布在当前序列块第32～37个块上，即B(32),B(33),B(34),…，B(47)，4个子序列块交迭分布，与0～4个子序列块的组成规律一致；11～15个子序列块分布在当前序列块第48～63个块上，即B(48),B(49),B(50),…，B(63)，4个子序列块交迭分布，与0～4个子序列块的组成规律一致。依次对该16个序列块的各自的前半个序列块的奇数序列块和后半个序列块的偶数序列块依次进行对应置换。例如第1个序列块(B(0),B(4),B(8)，B(12))，B(4)和B(8)置换即可；
步骤六：此时B0序列块(B(0),B(1),B(2),…，B(63))完成了6位原位逆序。由于假定论述的对象是1024点的中间6位的部分逆序，进行了步骤三、步骤四、步骤五总共三个步骤，步长因子依次取值为1、2、4后，即结束了置换过程。这三个步骤本质上是类似，其核心步骤是：先根据步长因子的值，确定子序列块的个数和子序列的组成，然后针对每个子序列块各自的前半个序列块的奇数序列块和后半个序列块的偶数序列块进行对应置换，所以可以看成看成主要依赖于步长因子的迭代步骤。该核心步骤的迭代次数为(R-4)/2，该值向下取整，例R＝11时的核心步骤的迭代次数也为3。针对其它点数待逆序的序列，核心步骤根据(R-4)/2进行调整，步长因子都是从1开始，递增规律为上一次迭代步骤的步长因子的2倍，而对应的子序列块个数也是从1开始，递增规律为上一次迭代步骤的子序列块个数的4倍，对应的子序列组成根据步长因子和子序列个数从序列块中按照所述类似规律选取；
步骤七：B1，B2，B3也完成和B0一样的6位原位逆序后，则意味着原始序列(A(0),A(1),A(2),…，A(1023)完成了针对中间6位逆序高/低2位不变的部分逆序。
在完成部分逆序后，对此时的序列首先进行三阶/四阶合并：若是奇数阶，则进行三阶合并，即从内存读取序列到数字信号处理器的寄存器进行前三阶运算之后，才写回内存；若是偶数阶，则进行四阶合并，即从内存读取序列到数字信号处理器的寄存器进行前四阶运算之后，才写回内存。由于前三/四阶计算涉及到旋转因子是比较简单的，可以不用从内存中获得，而且计算量较小。
在完成三阶/四阶合并后，进行中间阶的两阶合并：从内存读取序列和旋转因子到数字信号处理器寄存器的寄存器进行两阶运算后，才写回内存。
在完成中间阶的两阶合并循环后，最后进行两阶计算/最终逆序：从内存读取序列和旋转因子到数字信号处理器寄存器的寄存器进行最后两阶运算后，写回内存时，完成最终的高低part位逆序。在进行最后两阶计算/最终逆序之前，序列是以部分逆序存储在内存中。A序列((A(0),A(1),A(2),…，A(1023))被等分成为高位分别为(00,01,10,11)四个子序列B0，B1，B2，B3。而B0，B1，B2，B3的低2位又都以(00,01,10,11)的规律变化；从中间6位看，B0，B1，B2，B3则全部完成了中间6位的完全逆序，它们位序中间6位的是依次相同的，依次为序列(0,1,2,3,,…,63)的逆序，以下以R(i)表示i的二进制逆数数，共6位。则此时的B0序列位序状态为
(00R(0)00,00R(0)01,00R(0)10,00R(0)11,
00R(1)00,00R(1)01,00R(1)10,00R(1)11,
00R(2)00,00R(2)01,00R(2)10,00R(2)11,
00R(3)00,00R(3)01,00R(3)10,00R(3)11,
…
00R(63)00,00R(63)01,00R(63)10,00R(63)11)
此时的B1序列位序状态为
(01R(0)00,01R(0)01,01R(0)10,01R(0)11,
01R(1)00,01R(1)01,01R(1)10,01R(1)11,
01R(2)00,01R(2)01,01R(2)10,01R(2)11,
01R(3)00,01R(3)01,01R(3)10,01R(3)11,
…
01R(63)00,01R(63)01,01R(63)10,01R(63)11)
此时的B2序列位序状态为
(10R(0)00,10R(0)01,10R(0)10,10R(0)11,
10R(1)00,10R(1)01,10R(1)10,10R(1)11,
10R(2)00,10R(2)01,10R(2)10,10R(2)11,
10R(3)00,10R(3)01,01R(3)10,01R(3)11,
…
10R(63)00,10R(63)01,10R(63)10,10R(63)11)
此时的B3的序列位序状态为
(11R(0)00,11R(0)01,11R(0)10,11R(0)11,
11R(1)00,11R(1)01,11R(1)10,11R(1)11,
11R(2)00,11R(2)01,11R(2)10,11R(2)11,
11R(3)00,11R(3)01,11R(3)10,11R(3)11,
…
11R(63)00,11R(63)01,11R(63)10,11R(63)11)
可见最终的高/低2位的逆序需要在位序中间6位的逆序数相同的R(i)上通过寄存器进行内部置换，即:
(00R(i)00,00R(i)01,00R(i)10,00R(i)11)，
(01R(i)00,01R(i)01,01R(i)10,01R(i)11)，
(10R(i)00,10R(i)01,10R(i)10,10R(i)11),
(11R(i)00,11R(i)01,11R(i)10,11R(i)11)
进行内部数据交换即可完成最终的逆序过程，即
(00R(i)00,10R(i)00,01R(i)00,11(i)00)，
(00R(i)10,10R(i)10,01R(i)10,11R(i)10)，
(00R(i)01,10R(i)01,01R(i)01,11R(i)01),
(00R(i)11,10R(i)11,01R(i)11,11R(i)11)
依次遍历R(0)，R(1),R(2),…，R(63),根据以上置换规律，即可完成最终的逆序过程。
在三阶/四阶合并、二阶合并的计算过程中涉及到取数规律的设置。如表1所示，以1042点FFT的高/低2位部分逆序为例，说明了基于部分逆序时各阶运算的寻址步长规律。
表1 1024点FFT部分逆序的寻址步长表
基于置换的原位逆序方法，是由以下两个推导步骤导出的。
推导步骤一：假定A点FFT逆序的方法是已知的，则获得4*A点FFT逆序的途径可以通过将4*A点转换成四个块序列的(0A0，0A1，1A0，1A1)的整体逆序，整体逆序完成之后；再对四个块的对应A点FFT用已知方法进行处理。整体逆序是通过图1所示的方法，进行整体置换。
推导步骤二：步骤一是可递归的。由A点已知的逆序方法既可以知道4*A点的逆序方法，也可以知道16*A(即4*(4*A))点的逆序方法。如图2所示，将16个块中首位和末尾相同的块依次组合成四个块(0B0，0B1，1B0，1B1)；用推导步骤一的方法对该四个块进行整体置换；然后对四个块中每个B块用推导步骤一的方法进行整体置换；最后再用已知的方法对A进行变换。
下面对基于置换的原位逆序算法的时间空间复杂度进行分析：
假定最内层循环的一次置换为基本单位，则算法复杂度为O((R/2)*(n/4))次，其中R/2为整体置换次数。当是1024点时，置换次数为5*256＝1280次；和非原位逆序运算相比，该原位逆序运算节约1022个空间，置换次数也多了((R/2)-1)＝4次。该算法空间复杂度则降低到最小，只需要2个单位的临时寄存器，即完成了整个逆序过程。
考虑到DSP处理器提供的连续地址空间的高数据访存并行性，提出面向高数据并行的块原位部分逆序方法。为了充分利用DSP处理器的连续个地址的数据访问能力，结合FFT旋转因子访问的特点，考虑部分逆序。部分逆序是指保持逆序数的高/低若干位不变，其它位数则对应逆序，如图3所示。
假定DSP处理器的连续多地址的数据访问能力为B个复数字，则部分逆序位数至少为logB，即保持高低LogB位不变，其它位数逆序。
该算法的思想是把首先把输入序列按照高part位连续划分成为2part个块；低part位则作为一个整体参与块置换；则部分逆序转化为求解(R-2*part)个中间位的逆序操作，此时则可利用原位逆序算法完成输入序列的部分逆序过程。该算法的时间复杂度为O(((R-2*part)/2)*2R-part-2)，单位为两个复数的置换操作。
部分逆序后输入序列在内存的位序如图4所示，该图是高/低位为2的部分逆序后的序列。输入序列被等分成为高位分别为(00,01,10,11)四个连续块，其中A，B，C，…，为2R -2*part个输入序列的逆序位序。可见完成最终的逆序过程只需依次遍历A，B，C，…，对中间(R-4)位位序相等的四组数据(00A00,00A01,00A10,00A11),(01A00,01A01,01A10,01A11)，(10A00,10A01,10A10,10A11),(11A00,11A01,11A10,11A11)进行内部数据交换即可完成最终的逆序过程。
假定一款数字信号处理器是具有支持8个连续地址(字寻址)的高数据并行处理能力、四个簇、支持多发射的VLIW结构的浮点运算信号处理器，每个核心有64个通用寄存器。其单核单周期可以读取两个浮点复数，存储一个浮点复数。由于单个形如(a+bj)±(c+dj)*(Wr+Wij)的蝶形运算只需要4次乘法和6次加减法即可完成，可以看出该数字信号处理器芯片单核可以在一个周期内实现一次浮点蝶形运算。下面看看数据访存的要求，为了达到单周期实现单个蝶形运算的目的，要求芯片能够支持单周期内读取三个复数(两个运算数、一个旋转因子)，单周期内存储两个复数(两个运算结果)，这已经超出了魂芯DSP芯片的能力。如果采取两阶合并策略的话，两个周期内要求芯片的I/O能力降低为：读取四个复数(两个运算数、两个旋转因子)，存储两个复数(两个运算结果)，则可以满足这一要求。故而，面向该款数字信号处理的高效FFT程序中，各处理循环至少为两阶合并以充分发挥芯片性能。
FFT算法的前两级为特殊旋转因子，单个蝶形运算只需要4次加减法。这样前后两级蝶形运算共8次加减法，一个周期即可完成。但再次遇到了数据访存能力不足的问题，继续采用多阶合并的解决办法。由于该款芯片核具有64个寄存器，合并阶数最高可以达到4阶。考虑到支持32点FFT点数，如果前几阶按4阶合并，则其最后只有独立1阶，无法实现最优。故而采用判断奇偶阶的处理办法。对于偶数阶(最小点数为64点)，采用4阶合并；对于奇数阶(最小点数32点)，采用3阶合并。
FFT的算法复杂度为N/2*log(N)次蝶形运算，由于设计的FFT算法经过阶合并后魂芯DSP的单核每周期可以进行一个蝶形运算，而总共有魂芯DSP有四个核心，所以该高效的FFT算法出去部分逆序的时间复杂度为O(N/8*log(N))。部分逆序的时间复杂度为O(((R-2*part)/2)*2R-part-2)＝O(((R-4)/2)*2R-4)。
所以基于置换的部分逆序的FFT算法的时间复杂度为O(N/8*log(N))+O(((R-4)/2)*2R-4)，空间复杂度为O(N)。如果采用非原址的部分逆序的FFT算法，则时间复杂度为N/8*log(N)，空间复杂度为O(2*N)。
和非原址的部分逆序的FFT算法相比，基于置换的原址部分逆序的FFT算法的性能降低的比率为：O(((R-4)/2)*2R-4)/O(N/8*log(N))＝(R-4)/4R(在该款数字信号处理平台上R的取值范围为[5，17])。
可见，随着FFT点数的增大，基于置换的原址部分逆序的FFT算法效率也随之降低，降低的百分比是从R＝5(32点FFT)的5％，逐步增大到R＝17(131072点FFT)的20％。而空间复杂度则从O(2*N)降低到O(N)。尤其是在大点数时R＝17(131072点FFT)时，虽然时间性能降低了20％，但空间复杂度则节约了131072*8bytes＝1M bytes。相比较于空间需求的大幅降低，算法时间效率降低的20％是可接受的；空间复杂度的降低对于信号处理器的较小的内存空间无疑是重要的。
综上所述，本发明针对数字信号处理器要求的高数据并行处理能力，创新地提出基于块原位的部分逆序方法，在此基础上，构建高效的原位FFT算法，该FFT优化方法在保证时间效率的同时，大幅降低了对内存空间的需求。
Claims (3)
1.一种基于数字信号处理的高效并行处理优化方法，其特征在于该方法包括下列顺序的步骤：
(1)进行基于置换的原位逆序的除高/低part之外的部分逆序；
(2)三阶/四阶合并；
(3)中间的二阶合并循环；
(4)最后两阶计算/最终逆序；
所述基于置换的原位逆序分为以下迭代步骤：
步骤一：把N个复数的输入序列(A(0),A(1),A(2),…，A(N-1))看成一个步长因子为1组成的1个子序列，把前N/2序列的奇数序列和后N/2序列的偶数序列对应置换，即前半部分序列(A(0),A(1),A(2),…，A(N/2-1))的奇数部分A(1),A(3),A(5)，…，A(N/2-1)和后半部分序列A(N/2),A(N/2+1),A(N/2+2),…，A(N-1)的偶数部分A(N/2),A(N/2+2),…，A(N-2)依次对应置换；
步骤二：把步骤一置换后的序列看成是一个步长因子为2组成的4个序列,其中前两个子序列在前N/2位置上，两个子序列交迭分布,即这两个子序列为(A(0),A(2),…，A(N/2-2))和(A(1),A(3),…，A(N/2-1))；后两个子序列在后N/2位置上，两个序列交迭分布,即这两个子序列为(A(N/2),A(N/2+2),…，A(N-2))和(A(N/2+1),A(N/2+3),…，A(N-1))；把4个子序列各自的前半部分的奇数部分和后半部分的偶数部分依次对应置换；
步骤三：把步骤二置换后的序列看成是一个步长因子为4组成的16个子序列，其中0～3个子序列在位置0～(N/4-1)上，四个子序列交迭分布,即这四个子序列为(A(0),A(0+4),A(0+8)，…，A(N/4-4))，(A(1),A(1+4),A(1+8),…，A(N/4-3))，(A(2),A(2+4),A(2+8),…，A(N/4-2))，(A(3),A(3+4),A(3+8),…，A(N/4-1))；其中4～7个子序列在位置N/4～(N/2-1)上，四个子序列交迭分布,即这四个子序列为(A(N/4+0),A(N/4+0+4),A(N/4+0+8)，…，A(N/2-4))，(A(N/4+1),A(N/4+1+4),A(N/4+1+8),…，A(N/2-3))，(A(N/4+2),A(N/4+2+4),A(N/4+2+8),…，A(N/2-2))，(A(N/4+3),A(N/4+3+4),A(N/4+3+8),…，A(N/2-1))；其中8～11个子序列在位置N/2～(3*N/4-1)上，四个子序列交迭分布,即这4个子序列为(A(N/2+0),A(N/2+0+4),A(N/2+0+8)，…，A(3*N/4-4))，(A(N/2+1),A(N/2+1+4),A(N/2+1+8),…，A(3*N/4-3))，(A(N/2+2),A(N/2+2+4),A(N/2+2+8),…，A(3*N/4-2))，(A(N/2+3),A(N/2+3+4),A(N/2+3+8),…，A(3*N/4-1))；其中11～15个子序列在位置3*N/4～(N-1)上，四个子序列交迭分布,即这4个子序列为(A(3*N/4+0),A(3*N/4+0+4),A(3*N/4+0+8)，…，A(N-4))，(A(3*N/4+1),A(3*N/4+1+4),A(3*N/4+1+8),…，A(N-3))，(A(3*N/4+2),A(3*N/4+2+4),A(3*N/4+2+8),…，A(N-2))，(A(3*N/4+3),A(3*N/4+3+4),A(3*N/4+3+8),…，A(N-1))；把16个子序列各自的前半部分的奇数部分和后半部分的偶数部分依次对应置换；
步骤四：把步骤三置换后的序列按照步长因子是上一次迭代步长因子的2倍、子序列个数是上一次迭代子序列个数的4倍的方法划分出若干个子序列，依次对该若干个子序列进行各自的前半部分的奇数部分和后半部分偶数部分的进行置换；
步骤五：把整个序列的多个子序列每次所完成的置换作为一个迭代步骤，那么总迭代步骤为R/2次，该值向下取整，R为LogN；
所述原位部分逆序是指：首先把输入序列按照高part位连续划分成为几个大块；其次将低part位则作为一个整体小块参与块原位置换；最后，部分逆序转化为求解(R-2*part)个中间位的逆序操作，此时则可利用原位逆序方法完成输入序列的部分逆序过程；
在完成部分逆序后，对此时的序列首先进行三阶/四阶合并：若是奇数阶，则进行三阶合并，即从内存读取序列到数字信号处理器的寄存器进行前三阶运算之后，才写回内存；若是偶数阶，则进行四阶合并，即从内存读取序列到数字信号处理器的寄存器进行前四阶运算之后，才写回内存。
2.根据权利要求1所述的一种基于数字信号处理的高效并行处理优化方法，其特征在于：在完成三阶/四阶合并后，进行中间阶的两阶合并：从内存读取序列和旋转因子到数字信号处理器寄存器的寄存器进行两阶运算后，才写回内存。
3.根据权利要求1所述的一种基于数字信号处理的高效并行处理优化方法，其特征在于：在完成中间阶的两阶合并循环后，最后进行两阶计算/最终逆序：从内存读取序列和旋转因子到数字信号处理器寄存器的寄存器进行最后两阶运算后，写回内存时，完成最终的高低part位逆序。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201410341689.4A CN104142811B (zh) | 2014-07-18 | 2014-07-18 | 一种基于数字信号处理的高效并行处理优化方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201410341689.4A CN104142811B (zh) | 2014-07-18 | 2014-07-18 | 一种基于数字信号处理的高效并行处理优化方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN104142811A CN104142811A (zh) | 2014-11-12 |
| CN104142811B true CN104142811B (zh) | 2017-02-01 |
Family
ID=51851995
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201410341689.4A Active CN104142811B (zh) | 2014-07-18 | 2014-07-18 | 一种基于数字信号处理的高效并行处理优化方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN104142811B (zh) |
Families Citing this family (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN104820581B (zh) * | 2015-04-14 | 2017-10-10 | 广东工业大学 | 一种fft和ifft逆序数表的并行处理方法 |
Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101211333A (zh) * | 2006-12-30 | 2008-07-02 | 北京邮电大学 | 一种信号处理方法、装置和系统 |
| CN101694648A (zh) * | 2009-08-28 | 2010-04-14 | 曙光信息产业(北京)有限公司 | 傅里叶变换处理方法和装置 |
| US7864900B2 (en) * | 2006-10-30 | 2011-01-04 | Al-Eidan Abdullah A | Communication system for sending and receiving digital data |
| CN102375805A (zh) * | 2011-10-31 | 2012-03-14 | 中国人民解放军国防科学技术大学 | 面向向量处理器的基于simd的fft并行计算方法 |
Family Cites Families (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20080052336A1 (en) * | 2006-08-22 | 2008-02-28 | Andre Kaufmann | Method and Apparatus for Providing FFT-Based Signal Processing with Reduced Latency |
-
2014
- 2014-07-18 CN CN201410341689.4A patent/CN104142811B/zh active Active
Patent Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US7864900B2 (en) * | 2006-10-30 | 2011-01-04 | Al-Eidan Abdullah A | Communication system for sending and receiving digital data |
| CN101211333A (zh) * | 2006-12-30 | 2008-07-02 | 北京邮电大学 | 一种信号处理方法、装置和系统 |
| CN101694648A (zh) * | 2009-08-28 | 2010-04-14 | 曙光信息产业(北京)有限公司 | 傅里叶变换处理方法和装置 |
| CN102375805A (zh) * | 2011-10-31 | 2012-03-14 | 中国人民解放军国防科学技术大学 | 面向向量处理器的基于simd的fft并行计算方法 |
Non-Patent Citations (3)
| Title |
|---|
| 利用逆序循环实现FFT运算中倒序算法的优化;方志红 等;《信号处理》;20041030;第20卷(第5期);533-535 * |
| 基于查找表的单基FFT原址倒序算法;王海兵 等;《清华大学学报》;20080115;第48卷(第1期);43-45 * |
| 逆序循环在原位FFT中的应用;方志红 等;《全国第21届计算机技术与应用学术会议（CACIS·2010）暨全国第2届安全关键技术与应用学术会议论文集》;20100820;43-46 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN104142811A (zh) | 2014-11-12 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Qiao et al. | AtomLayer: A universal ReRAM-based CNN accelerator with atomic layer computation | |
| Kestur et al. | Towards a universal FPGA matrix-vector multiplication architecture | |
| Inoue et al. | SIMD-and cache-friendly algorithm for sorting an array of structures | |
| CN102375805B (zh) | 面向向量处理器的基于simd的fft并行计算方法 | |
| US20230333815A1 (en) | Concurrent multi-bit adder | |
| CN103984527A (zh) | 优化稀疏矩阵向量乘提升不可压缩管流模拟效率的方法 | |
| CN104331497A (zh) | 一种利用向量指令并行处理文件索引的方法及装置 | |
| CN101986264A (zh) | 用于simd向量微处理器的多功能浮点乘加运算装置 | |
| CN101859241B (zh) | 基于全展开的全流水128位精度浮点累加器 | |
| TWI604379B (zh) | 用於ｋ最近相鄰者搜尋之系統、設備及方法 | |
| CN103699516A (zh) | 向量处理器中基于simd的并行fft/ifft蝶形运算方法及装置 | |
| CN102495721A (zh) | 一种支持fft加速的simd向量处理器 | |
| CN109240644B (zh) | 一种用于伊辛芯片的局部搜索方法及电路 | |
| CN110543291A (zh) | 有限域大整数乘法器及基于ssa算法的大整数乘法的实现方法 | |
| CN104142811B (zh) | 一种基于数字信号处理的高效并行处理优化方法 | |
| CN109614145A (zh) | 一种处理器核心结构及数据访存方法 | |
| CN110019184A (zh) | 一种压缩和解压缩有序整数数组的方法 | |
| Yang et al. | GQNA: Generic quantized DNN accelerator with weight-repetition-aware activation aggregating | |
| CN110990299A (zh) | 非规整组相联cache组地址映射方法 | |
| US20130262819A1 (en) | Single cycle compare and select operations | |
| CN209879493U (zh) | 乘法器 | |
| CN104317554A (zh) | 用于simd处理器的寄存器文件数据读写装置和方法 | |
| CN111368250B (zh) | 基于傅里叶变换/逆变换的数据处理系统、方法及设备 | |
| CN118377998B (zh) | 一种fft指令实现方法及系统 | |
| CN113031911A (zh) | 乘法器、数据处理方法、装置及芯片 |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| C06 | Publication | ||
| PB01 | Publication | ||
| C10 | Entry into substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| C14 | Grant of patent or utility model | ||
| GR01 | Patent grant | ||
| TR01 | Transfer of patent right | ||
| TR01 | Transfer of patent right |
Effective date of registration: 20191016 Address after: 5 / F, airborne center, 38 new area, No. 199, Xiangzhang Avenue, hi tech Zone, Hefei City, Anhui Province 230000 Patentee after: Anhui core Century Technology Co., Ltd. Address before: 230088, 199, camphor Road, hi tech Zone, Anhui, Hefei Patentee before: No.38 Inst., China Electronic Sci. & Tech. Group Co. |