用于ARM上的FFT与IFFT源代码C语言，不依赖特定平台

http://blog.csdn.net/syrchina/article/details/6670517

代码在2011年全国电子大赛结束后（2011年9月3日）发布，多个版本，注释详细。

[cpp]view plaincopy 
   
 /******************************************************************************* 
 ** 程序名称：快速傅里叶变换(FFT)  
** 程序描述：本程序实现快速傅里叶变换  
** 程序作者：宋元瑞  
** 最后修改：2011年4月5日  
*******************************************************************************/  
 #include <stdio.h>  
 #include <math.h>  
   
 #define PI 3.141592653589   //圆周率，12位小数   
 #define N 8                 //傅里叶变换的点数   
 #define M 3                 //蝶形运算的级数，N = 2^M   
 typedef double ElemType;    //原始数据序列的数据类型,可以在这里设置  
   
typedef struct              //定义复数结构体   
 {  
     ElemType real,imag;  
 }complex;  
 complex data[N];            //定义存储单元，原始数据与负数结果均使用之   
 ElemType result[N];         //存储FFT后复数结果的模  
//变址   
 void ChangeSeat(complex *DataInput)  
 {  
     int nextValue,nextM,i,k,j=0;  
     complex temp;  
       
     nextValue=N/2;                  //变址运算，即把自然顺序变成倒位序，采用雷德算法  
     nextM=N-1;  
     for (i=0;i<nextM;i++)  
     {  
         if (i<j)                 //如果i<j,即进行变址  
         {  
             temp=DataInput[j];  
             DataInput[j]=DataInput[i];  
             DataInput[i]=temp;  
         }  
         k=nextValue;                //求j的下一个倒位序  
         while (k<=j)             //如果k<=j,表示j的最高位为1  
         {  
             j=j-k;                  //把最高位变成0  
             k=k/2;                  //k/2，比较次高位，依次类推，逐个比较，直到某个位为0  
         j=j+k;                      //把0改为1  
     }                                         
 }  
/* 
//变址  
void ChangeSeat(complex *DataInput) 
{ 
    complex Temp[N]; 
    int i,n,New_seat; 
    for(i=0; i<N; i++) 
    { 
        Temp[i].real = DataInput[i].real; 
        Temp[i].imag = DataInput[i].imag; 
    } 
    for(i=0; i<N; i++) 
    { 
        New_seat = 0; 
        for(n=0;n<M;n++) 
        { 
            New_seat = New_seat | (((i>>n) & 0x01) << (M-n-1)); 
        } 
        DataInput[New_seat].real = Temp[i].real; 
        DataInput[New_seat].imag = Temp[i].imag; 
} 
*/  
//复数乘法   
 complex XX_complex(complex a, complex b)  
     complex temp;  
       
     temp.real = a.real * b.real-a.imag*b.imag;  
     temp.imag = b.imag*a.real + a.imag*b.real;  
return temp;  
 }  
//FFT  
void FFT(void)  
     int L=0,B=0,J=0,K=0;  
int step=0;  
     ElemType P=0,T=0;  
     complex W,Temp_XX;  
     //ElemType TempResult[N];  
     ChangeSeat(data);  
     for(L=1; L<=M; L++)  
     {  
         B = 1<<(L-1);//B=2^(L-1)  
for(J=0; J<=B-1; J++)  
             P = (1<<(M-L))*J;//P=2^(M-L) *J   
             step = 1<<L;//2^L  
             for(K=J; K<=N-1; K=K+step)  
             {  
                 W.real =  cos(2*PI*P/N);  
                 W.imag = -sin(2*PI*P/N);  
                   
                 Temp_XX = XX_complex(data[K+B],W);  
                 data[K+B].real = data[K].real - Temp_XX.real;  
                 data[K+B].imag = data[K].imag - Temp_XX.imag;  
                 data[K].real = data[K].real + Temp_XX.real;  
                 data[K].imag = data[K].imag + Temp_XX.imag;  
             }  
         }  
     }  
void IFFT(void)  
int step=0;  
     ElemType P=0,T=0;  
     complex W,Temp_XX;  
     //ElemType TempResult[N];  
     ChangeSeat(data);  
for(L=1; L<=M; L++)  
         B = 1<<(L-1);//B=2^(L-1)  
for(J=0; J<=B-1; J++)  
             P = (1<<(M-L))*J;//P=2^(M-L) *J   
             step = 1<<L;//2^L  
             for(K=J; K<=N-1; K=K+step)  
             {  
                 W.real =  cos(2*PI*P/N);  
                 W.imag =  sin(2*PI*P/N);//逆运算，这里跟FFT符号相反   
                   
                 Temp_XX = XX_complex(data[K+B],W);  
                 data[K+B].real = data[K].real - Temp_XX.real;  
                 data[K+B].imag = data[K].imag - Temp_XX.imag;  
                 data[K].real = data[K].real + Temp_XX.real;  
                 data[K].imag = data[K].imag + Temp_XX.imag;  
             }  
     }  
 int main(int argc, char *argv[])  
int i = 0;  
for(i=0; i<N; i++)//制造输入序列   
         data[i].real = sin(2*PI*i/N);  
         printf("%lf ",data[i]);  
     printf("\n\n");  
     FFT();//进行FFT计算   
for(i=0; i<N; i++)  
         {printf("%lf ",sqrt(data[i].real*data[i].real+data[i].imag*data[i].imag));}  
     IFFT();//进行FFT计算   
     printf("\n\n");  
for(i=0; i<N; i++)  
         {printf("%lf ",data[i].real/N);}  
     printf("\n");  
/*for(i=0; i<N; i++) 
        {printf("%lf ",data[i].imag/N);} 
    printf("\n");*/  
/*for(i=0; i<N; i++) 
  
return 0;  
 }  

copy 
    
   

性能

** 程序版本：V6.5
** 最后修改：2011年5月16日
#define PI 3.14159265358979323846264338327950288419716939937510 //圆周率，50位小数
#define PI2 6.28318530717958647692528676655900576839433879875021
#define N 1024 //傅里叶变换的点数
#define M 10 //蝶形运算的级数，N = 2^M
#define Npart2 512 //创建正弦函数表时取PI的1/2
#define Npart4 256 //创建正弦函数表时取PI的1/4
#define FFT_RESULT(x) (sqrt(data[x].real*data[x].real+data[x].imag*data[x].imag))
#define IFFT_RESULT(x) (data[x].real/N)
float ElemType; 在这里设置
//定义复数结构体
ElemType real,imag;
}complex;
complex data[N]; //定义存储单元，原始数据与负数结果均使用之
ElemType SIN_TABLE[Npart4+1];
//产生模拟原始数据输入
void InputData(int i;
data[i].real = sin(2*PI*i/N); //正弦波
data[i].imag = 0;
//创建正弦函数表
void CREATE_SIN_TABLE(int i=0;
for(i=0; i<=Npart4; i++)
SIN_TABLE[i] = sin(PI*i/Npart2);//SIN_TABLE[i] = sin(PI2*i/N);
ElemType Sin_find(ElemType x)
int i = (int)(N*x);
i = i>>1;
if(i>Npart4)//注意：i已经转化为0~N之间的整数了！
{//不会超过N/2
i = Npart2 - i;//i = i - 2*(i-Npart4);
}
return SIN_TABLE[i];
ElemType Cos_find(ElemType x)
if(i<Npart4) //i = Npart4 - i;
return SIN_TABLE[Npart4 - i];
}
else//i>Npart4 && i<N/2
//i = i - Npart4;
return -SIN_TABLE[i - Npart4];
//变址
void ChangeSeat(complex *DataInput)
nextValue=N/2; //变址运算，即把自然顺序变成倒位序，采用雷德算法
nextM=N-1;
for (i=0;i<nextM;i++)
temp=DataInput[j];
DataInput[j]=DataInput[i];
DataInput[i]=temp;
k=nextValue; //求j的下一个倒位序
j=j-k; //把最高位变成0
k=k/2; //k/2，比较次高位，依次类推，逐个比较，直到某个位为0
j=j+k; //把0改为1
}
}
//复数乘法
/*complex XX_complex(complex a, complex b)
{
complex temp;
temp.real = a.real * b.real-a.imag*b.imag;
temp.imag = b.imag*a.real + a.imag*b.real;
return temp;
}*/
//FFT运算函数
int step=0, KB=0;
//ElemType P=0;
ElemType angle;
ChangeSeat(data);//CREATE_SIN_TABLE();
step = 1<<L; B = step>>1;for(J=0; J<B; J++)
//P = (1<<(M-L))*J;//P=2^(M-L) *J
angle = (double)J/B; //这里还可以优化
W.imag = -Sin_find(angle); //用C++该函数课声明为inline
W.real = Cos_find(angle); //用C++该函数课声明为inline
//W.real = cos(angle*PI);
//W.imag = -sin(angle*PI);
for(K=J; K<N; K=K+step)
KB = K + B;
//Temp_XX = XX_complex(data[KB],W);
//用下面两行直接计算复数乘法，省去函数调用开销
Temp_XX.real = data[KB].real * W.real-data[KB].imag*W.imag;
Temp_XX.imag = W.imag*data[KB].real + data[KB].imag*W.real;
data[KB].real = data[K].real - Temp_XX.real;
data[KB].imag = data[K].imag - Temp_XX.imag;
//IFFT运算函数
W.imag = Sin_find(angle); W.real = Cos_find(angle); //W.real = cos(angle*PI);
//W.imag = -sin(angle*PI);
for(K=J; K<N; K=K+step)
KB = K + B;
//用下面两行直接计算复数乘法，省去函数调用开销
Temp_XX.real = data[KB].real * W.real-data[KB].imag*W.imag;
Temp_XX.imag = W.imag*data[KB].real + data[KB].imag*W.real;
data[KB].real = data[K].real - Temp_XX.real;
data[KB].imag = data[K].imag - Temp_XX.imag;
//主函数
CREATE_SIN_TABLE(); //创建正弦函数表 ,这句只需在程序开始时执行一次
InputData(); //输入原始数据，此处用公式模拟；实际应用时为AD采样数据
FFT(); //进行 FFT计算
{printf("%f ",FFT_RESULT(i));}/**/
//进行 IFFT计算
copy

用于ARM上的FFT与IFFT源代码C语言，不依赖特定平台

相关推荐