最近在做飞腾上的FFT优化，记录一下以后用。
目前实现了基2FFT，使用arm提供的neon接口做了并行计算。算法原理网上很多，这里就不讲了，记录复数正向优化方法。
优化思路：
第一层蝶形计算：
第一层的蝶形因子都是1，蝶形下半部分就无需乘蝶形因子，直接与蝶形上半部分做加减运算即可。
由于第一层蝶形上下部分跨度为1（即连续），直接加载4个数据到向量寄存器无法实现并行计算，改用vld2q_f32()对输入信号进行加载，存入到float32x4x2类型的变量中。具体形式如下图：
第一行4个数据为蝶形的上半部分，第二行四个数据为蝶形的下半部分，可直接进行结果的计算。
复数部分与实数部分一样。
第二层蝶形计算
第二层蝶形有两种旋转因子分别为1和-i，蝶形上下部分间隔为2，可以使用vld4q_f32()将数据加载到float32x4x4_t的变量中，具体形式如下图：
第一行与第三行为同组蝶形的上下两部分（间隔为2），第二行与第四行为第二种蝶形的上下部分。对于旋转因子为1的部分，计算方式与第一层相同。对于第二种蝶形，旋转因子为-i，可以直接得到下半部分与蝶形因子的计算结果而无需乘法。
第三层及之后的
蝶形计算
第三层旋转因子有4个，蝶形上下部分间隔为4，可以选择使用vld1q_f32()加载数据到float32x4变量中，再加载之后的4个到第二个float32x4变量中。两个变量分别存储蝶形上下两部分的数据。然后进行常规的复数乘法运算与加减法。
后来发现这种方式的访存次数太多，考虑到armV8有32个向量寄存器，cache line大小为64B即16个float，利用vld1q_f32_x4()接口加载一个缓存行的数据，用4个向量寄存器存着，计算完再利用vst1q_f32_x4()存回去，将8次访存（4次取4次存）减少为2次访存（1次取1次存）。
实现了这个想法后，在飞腾CPU上利用g++编译发现不认识vld1q_f32_x4()与vst1q_f32_x4()接口!!
在我本机（macbook pro 2017)编码是可以看到<arm_neon.h>库中有这些接口的。
最后使用4个vld1q_f32()对16个数据进行了加载，4个vst1q_f32()进行存储。
在飞腾CPU上进行测试，编译器打开O3优化，进行了1000次8192点复数FFT变换，性能对比FFTW慢50%左右。仍然有很大的提升空间。
附上代码：（旋转因子是另一个函数专门算的，作为入参数传到了FFT函数中）
#include <iostream>
#include "FFTNeon.h"
#include <arm_neon.h>
#include <math.h>using namespace std;#define PI 3.1415926535inline float32x4x4_t vld1q_f32_x4(float *A){float32x4x4_t ret;ret.val[0] = vld1q_f32(A);ret.val[1] = vld1q_f32(A+4);ret.val[2] = vld1q_f32(A+8);ret.val[3] = vld1q_f32(A+12);return ret;
}inline float32x4x2_t vld1q_f32_x2(float *A){float32x4x2_t ret;ret.val[0] = vld1q_f32(A);ret.val[1] = vld1q_f32(A+4);return ret;
}inline void vst1q_f32_x4(float *A,float32x4x4_t VA){vst1q_f32(A,VA.val[0]);vst1q_f32(A+4,VA.val[1]);vst1q_f32(A+8,VA.val[2]);vst1q_f32(A+12,VA.val[3]);
}
void FFTC2Neon(float *dataR,float *dataI,float *TwiddleR,float *TwiddleI,int N,int M)
{int i,j,k,r;int p,L,Stride; //L表示第L级蝶形计算，stride步长unsigned int J,K,m,n;float Tr,Ti,temp;float *PTwiddleR;float *PTwiddleI;float TwiddleL3[8] = {1,0.707107,0,-0.707107,0,-0.707107,-1,-0.707107};float32x4x2_t VL1,VTempL1,VL4RUp,VL4RDown,VL4IUp,VL4IDown;float32x4x4_t VL2R,VTempL2,VL2I,VTwiddleR,VTwiddleI,VTempI,VRDown,VIDown;float32x4_t VL3RUp,VL3RDown,VL3IUp,VL3IDown,VTempL3;//存储4组蝶形计算的输入float32x4_t VL3TwiddleR,VL3TwiddleI;//存储旋转因子//倒序for(i=0;i< N;i++) { //获取I的反序j J=0;for(k=0;k<(M/2+0.5);k++) //k表示操作{//反序m=1;//m位最低位为1的二进制数//n=(unsigned int)pow(2,M-1);//n是第m位为1的二进制数n = 1<<(M-1);m <<= k; n >>= k; // F0=i & n;// F1=i & m;if(i & n) J=J | m; if(i & m) J=J | n;}if(i<J){temp = dataR[i];dataR[i] = dataR[J];dataR[J] = temp;//虚部 temp = dataI[i];dataI[i] = dataI[J];dataI[J] = temp;} } if(N>31){//level1:相邻2个为一组计算，W=1;实数level1都是实数计算,一次计算8个数据for( i = 0;i < N/8; ++i){VL1 = vld2q_f32(dataR + 8*i);VL4IDown = vld2q_f32(dataI + 8*i);VTempL1.val[0] = vaddq_f32(VL1.val[0],VL1.val[1]);VTempL1.val[1] = vsubq_f32(VL1.val[0],VL1.val[1]);vst2q_f32(dataR + 8*i,VTempL1);VL4RDown.val[0] = vaddq_f32(VL4IDown.val[0],VL4IDown.val[1]);VL4RDown.val[1] = vsubq_f32(VL4IDown.val[0],VL4IDown.val[1]);vst2q_f32(dataI + 8*i,VL4RDown);}//level2:计算输入间隔为2，同一蝶形间隔为4，W有2种分别为1和，同一种蝶形N/4组。一次计算16个数据for(i = 0;i < N/16; ++i){VL2R = vld4q_f32(dataR + 16 * i);//第一个蝶形旋转因子为1VTempL2.val[0] = vaddq_f32(VL2R.val[0],VL2R.val[2]);VL2I = vld4q_f32(dataI + 16 * i);VTempL2.val[2] = vsubq_f32(VL2R.val[0],VL2R.val[2]);VTempL3 = VL2I.val[0]; //VL2I.val[0] = vaddq_f32(VL2I.val[0],VL2I.val[2]);VL2I.val[2] = vsubq_f32(VTempL3,VL2I.val[2]);//下一个蝶形旋转因子为0-jVTempL2.val[1] = vaddq_f32(VL2R.val[1],VL2I.val[3]);VTempL2.val[3] = vsubq_f32(VL2R.val[1],VL2I.val[3]);vst4q_f32(dataR + 16 * i,VTempL2);VTempL3 = VL2I.val[1];VL2I.val[1] = vsubq_f32(VTempL3,VL2R.val[3]);VL2I.val[3] = vaddq_f32(VTempL3,VL2R.val[3]);vst4q_f32(dataI + i * 16,VL2I);}//Level3:4个w，跨度为4，全虚数运算for(i = 0; i < N/16;++i){//取出数据// VL3RUp = vld1q_f32(dataR + i * 8);// VL3RDown = vld1q_f32(dataR + i * 8 + 4);VL2R = vld1q_f32_x4(dataR + 16 * i);//位置：[0]UP,[1]Down,[2]Up,[3]Down// VL3IUp = vld1q_f32(dataI + i * 8);// VL3IDown = vld1q_f32(dataI + i * 8 +4);VL2I = vld1q_f32_x4(dataI + 16 * i);VL3TwiddleR = vld1q_f32(TwiddleL3);VL3TwiddleI = vld1q_f32(TwiddleL3+4);;VTempL3 = VL2R.val[1];VL2R.val[1] = vsubq_f32(vmulq_f32(VL2R.val[1],VL3TwiddleR),vmulq_f32(VL2I.val[1],VL3TwiddleI));VL2I.val[1] = vaddq_f32(vmulq_f32(VTempL3,VL3TwiddleI),vmulq_f32(VL2I.val[1],VL3TwiddleR));VL3TwiddleR = vaddq_f32(VL2R.val[1],VL2R.val[0]);//蝶形上部分实数VL3TwiddleI = vaddq_f32(VL2I.val[1],VL2I.val[0]);//蝶形上部分虚数VL2R.val[1] = vsubq_f32(VL2R.val[0],VL2R.val[1]);VL2I.val[1] = vsubq_f32(VL2I.val[0],VL2I.val[1]);VL2R.val[0] = VL3TwiddleR;VL2I.val[0] = VL3TwiddleI;// vst1q_f32(dataR + i * 8,VL3TwiddleR);// vst1q_f32(dataR + i * 8 + 4,VL3RDown);// vst1q_f32(dataI + i * 8,VL3TwiddleI);// vst1q_f32(dataI + i * 8 + 4,VL3IDown);VTempL3 = VL2R.val[3];VL2R.val[3] = vsubq_f32(vmulq_f32(VL2R.val[3],VL3TwiddleR),vmulq_f32(VL2I.val[3],VL3TwiddleI));VL2I.val[3] = vaddq_f32(vmulq_f32(VTempL3,VL3TwiddleI),vmulq_f32(VL2I.val[3],VL3TwiddleR));VL3TwiddleR = vaddq_f32(VL2R.val[3],VL2R.val[2]);//蝶形上部分实数VL3TwiddleI = vaddq_f32(VL2I.val[3],VL2I.val[2]);//蝶形上部分虚数VL2R.val[3] = vsubq_f32(VL2R.val[2],VL2R.val[1]);VL2I.val[3] = vsubq_f32(VL2I.val[2],VL2I.val[1]);VL2R.val[2] = VL3TwiddleR;VL2I.val[2] = VL3TwiddleI;vst1q_f32_x4(dataR + 16 * i,VL2R);vst1q_f32_x4(dataI + 16 * i,VL2I);}PTwiddleR = TwiddleR;PTwiddleI = TwiddleI;//Level4,8个旋转因子，上下跨度为8,一次加载16个为8组蝶形for(i=0;i<N/16;++i){VL2R = vld1q_f32_x4(dataR + 16 * i);//0,1Up,2,3DownVL3TwiddleR = VL2R.val[2];//暂存下半部分实数VL3TwiddleI = VL2R.val[3];//暂存下半部分实数VL2I = vld1q_f32_x4(dataI + 16 * i);VL4RDown = vld1q_f32_x2(PTwiddleR);//旋转因子实数VL4IDown = vld1q_f32_x2(PTwiddleI);//旋转因子虚数//计算下部分乘旋转因子VL2R.val[2] = vsubq_f32(vmulq_f32(VL3TwiddleR,VL4RDown.val[0]),vmulq_f32(VL4IDown.val[0],VL2I.val[2]));VL2I.val[2] = vaddq_f32(vmulq_f32(VL3TwiddleR,VL4IDown.val[0]),vmulq_f32(VL2I.val[2],VL4RDown.val[0]));VL2R.val[3] = vsubq_f32(vmulq_f32(VL3TwiddleI,VL4RDown.val[1]),vmulq_f32(VL4IDown.val[1],VL2I.val[3]));VL2I.val[3] = vaddq_f32(vmulq_f32(VL3TwiddleI,VL4IDown.val[1]),vmulq_f32(VL2I.val[3],VL4RDown.val[1]));//计算实数部分VTempL2.val[0] = vaddq_f32(VL2R.val[0],VL2R.val[2]);VTempL2.val[2] = vsubq_f32(VL2R.val[0],VL2R.val[2]);VTempL2.val[1] = vaddq_f32(VL2R.val[1],VL2R.val[3]);VTempL2.val[3] = vsubq_f32(VL2R.val[1],VL2R.val[3]);vst1q_f32_x4(dataR + i * 16,VTempL2);//计算虚数部分VL2R.val[0] = vaddq_f32(VL2I.val[0],VL2I.val[2]);VL2R.val[2] = vsubq_f32(VL2I.val[0],VL2I.val[2]);VL2R.val[1] = vaddq_f32(VL2I.val[1],VL2I.val[3]);VL2R.val[3] = vsubq_f32(VL2I.val[1],VL2I.val[3]);vst1q_f32_x4(dataI + i * 16,VL2R);}//Level 5-M:L表示第L层，每层2^(L-1)个蝶形PTwiddleR = TwiddleR + 8;PTwiddleI = TwiddleI + 8;for( L = 5;L<=M;++L){//旋转因子cos(2PI*j*2(M-L)/N)+sin()=cos(j*PI/2(L-1))+sin(),j = 0,1,2,...,2^(L-1) - 1J = 1<<(L-5);//每次循环计算16个蝶形Stride = 1<<(L-1);//跨度for(i = 0;i<(1<<(M-L));++i)//i控制多少种重复的蝶形计算for(k = 0;k<J; ++k){//J次循环算完所有不同蝶形//16个下半部分VRDown = vld1q_f32_x4(dataR + 16*k + i*(1<<L) + Stride);VIDown = vld1q_f32_x4(dataI + 16*k + i*(1<<L) + Stride);//16个蝶形因子VTwiddleR = vld1q_f32_x4(PTwiddleR);VTwiddleI = vld1q_f32_x4(PTwiddleI);//计算下半部分乘蝶形因子VTempL2.val[0] = vsubq_f32(vmulq_f32(VRDown.val[0],VTwiddleR.val[0]),vmulq_f32(VIDown.val[0],VTwiddleI.val[0]));VTempL2.val[1] = vsubq_f32(vmulq_f32(VRDown.val[1],VTwiddleR.val[1]),vmulq_f32(VIDown.val[1],VTwiddleI.val[1]));VTempL2.val[2] = vsubq_f32(vmulq_f32(VRDown.val[2],VTwiddleR.val[2]),vmulq_f32(VIDown.val[2],VTwiddleI.val[2]));VTempL2.val[3] = vsubq_f32(vmulq_f32(VRDown.val[3],VTwiddleR.val[3]),vmulq_f32(VIDown.val[3],VTwiddleI.val[3]));//16个上半部分VL2R = vld1q_f32_x4(dataR + 16*k + i*(1<<L));VL2I = vld1q_f32_x4(dataR + 16*k + i*(1<<L));//虚部VL4RDown.val[0] = vaddq_f32(vmulq_f32(VRDown.val[0],VTwiddleI.val[0]),vmulq_f32(VIDown.val[0],VTwiddleR.val[0]));VL4RDown.val[1] = vaddq_f32(vmulq_f32(VRDown.val[1],VTwiddleI.val[1]),vmulq_f32(VIDown.val[1],VTwiddleR.val[1]));VL4IDown.val[0] = vaddq_f32(vmulq_f32(VRDown.val[2],VTwiddleI.val[2]),vmulq_f32(VIDown.val[2],VTwiddleR.val[2]));VL4IDown.val[1] = vaddq_f32(vmulq_f32(VRDown.val[3],VTwiddleI.val[3]),vmulq_f32(VIDown.val[3],VTwiddleR.val[3]));//上半部分结果VRDown.val[0] = vaddq_f32(VL2R.val[0],VTempL2.val[0]);VIDown.val[0] = vaddq_f32(VL2I.val[0],VL4RDown.val[0]);VRDown.val[1] = vaddq_f32(VL2R.val[1],VTempL2.val[1]);VIDown.val[1] = vaddq_f32(VL2I.val[1],VL4RDown.val[1]);VRDown.val[2] = vaddq_f32(VL2R.val[2],VTempL2.val[2]);VIDown.val[2] = vaddq_f32(VL2I.val[2],VL4IDown.val[0]);VRDown.val[3] = vaddq_f32(VL2R.val[3],VTempL2.val[3]);VIDown.val[3] = vaddq_f32(VL2I.val[3],VL4IDown.val[1]);vst1q_f32_x4(dataR + 16*k + i*(1<<L),VRDown);vst1q_f32_x4(dataI + 16*k + i*(1<<L),VIDown);//下半部分结果VTwiddleR.val[0] = vsubq_f32(VL2R.val[0],VTempL2.val[0]);VTwiddleI.val[0] = vsubq_f32(VL2I.val[0],VL4RDown.val[0]);VTwiddleR.val[1] = vsubq_f32(VL2R.val[1],VTempL2.val[1]);VTwiddleI.val[1] = vsubq_f32(VL2I.val[1],VL4RDown.val[1]);VTwiddleR.val[2] = vsubq_f32(VL2R.val[2],VTempL2.val[2]);VTwiddleI.val[2] = vsubq_f32(VL2I.val[2],VL4IDown.val[0]);VTwiddleR.val[3] = vsubq_f32(VL2R.val[3],VTempL2.val[3]);VTwiddleI.val[3] = vsubq_f32(VL2I.val[3],VL4RDown.val[1]);vst1q_f32_x4(dataR + 16*k + i*(1<<L) + Stride,VTwiddleR);vst1q_f32_x4(dataI + 16*k + i*(1<<L) + Stride,VTwiddleI);}//旋转因子指针指向next层第一个旋转因子PTwiddleR = TwiddleR + Stride;PTwiddleI = TwiddleI + Stride;}return ;}
}
