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

时间：2023-09-01 20:14:08

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

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

/*******************************************************************************

** 程序名称：快速傅里叶变换(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=;

    complex temp;

    nextValue=N/;                  //变址运算，即把自然顺序变成倒位序，采用雷德算法

    nextM=N-;

    for (i=;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/;                    //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=,B=,J=,K=;

    int step=;

    ElemType P=,T=;

    complex W,Temp_XX;

    //ElemType TempResult[N];

    ChangeSeat(data);

    for(L=; L<=M; L++)

    {

        B = <<(L-);//B=2^(L-1)

        for(J=; J<=B-; J++)

        {

            P = (<<(M-L))*J;//P=2^(M-L) *J

            step = <<L;//2^L

            for(K=J; K<=N-; K=K+step)

            {

                W.real =  cos(*PI*P/N);

                W.imag = -sin(*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 L=,B=,J=,K=;

    int step=;

    ElemType P=,T=;

    complex W,Temp_XX;

    //ElemType TempResult[N];

    ChangeSeat(data);

    for(L=; L<=M; L++)

    {

        B = <<(L-);//B=2^(L-1)

        for(J=; J<=B-; J++)

        {

            P = (<<(M-L))*J;//P=2^(M-L) *J

            step = <<L;//2^L

            for(K=J; K<=N-; K=K+step)

            {

                W.real =  cos(*PI*P/N);

                W.imag =  sin(*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 = ;

    for(i=; i<N; i++)//制造输入序列

    {

        data[i].real = sin(*PI*i/N);

        printf("%lf ",data[i]);

    }

    printf("\n\n");

    FFT();//进行FFT计算

    printf("\n\n");

    for(i=; 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=; 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++)

        {printf("%lf ",sqrt(data[i].real*data[i].real+data[i].imag*data[i].imag)/N);}*/

    return ;

}
/*******************************************************************************

** 程序名称：快速傅里叶变换(FFT)

** 程序描述：本程序实现快速傅里叶变换

** 性能提升：修正了IFFT的bug,使用宏定义改变N大小

** 程序版本：V6.5

** 程序作者：宋元瑞

** 最后修改：2011年5月16日

*******************************************************************************/

#include <stdio.h>

#include <math.h>

#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)

typedef float ElemType;        //原始数据序列的数据类型,可以在这里设置

typedef struct                //定义复数结构体

{

    ElemType real,imag;

}complex;

complex data[N];            //定义存储单元，原始数据与负数结果均使用之

ElemType SIN_TABLE[Npart4+];

//产生模拟原始数据输入

void InputData(void)

{

    int i;

    for(i=; i<N; i++)//制造输入序列

    {

        data[i].real = sin(*PI*i/N);    //正弦波

        data[i].imag = ;

    }

}

//创建正弦函数表

void CREATE_SIN_TABLE(void)

{

    int i=;

    for(i=; 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>>;

    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)

{

    int i = (int)(N*x);

    i = i>>;

    if(i<Npart4)//注意：i已经转化为0~N之间的整数了！

    {//不会超过N/2

        //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)

{

    int nextValue,nextM,i,k,j=;

    complex temp;

    nextValue=N/;                  //变址运算，即把自然顺序变成倒位序，采用雷德算法

    nextM=N-;

    for (i=;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/;                    //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运算函数

void FFT(void)

{

    int L=,B=,J=,K=;

    int step=, KB=;

    //ElemType P=0;

    ElemType angle;

    complex W,Temp_XX;

    ChangeSeat(data);//变址

    //CREATE_SIN_TABLE();

    for(L=; L<=M; L++)

    {

        step = <<L;//2^L

        B = step>>;//B=2^(L-1)

        for(J=; 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;

                data[K].real = data[K].real + Temp_XX.real;

                data[K].imag = data[K].imag + Temp_XX.imag;

            }

        }

    }

}

//IFFT运算函数

void IFFT(void)

{

    int L=,B=,J=,K=;

    int step=, KB=;

    //ElemType P=0;

    ElemType angle;

    complex W,Temp_XX;

    ChangeSeat(data);//变址

    //CREATE_SIN_TABLE();

    for(L=; L<=M; L++)

    {

        step = <<L;//2^L

        B = step>>;//B=2^(L-1)

        for(J=; 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;

                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 = ;

    CREATE_SIN_TABLE();    //创建正弦函数表 ,这句只需在程序开始时执行一次 

    InputData();        //输入原始数据 ，此处用公式模拟；实际应用时为AD采样数据

    FFT();                //进行 FFT计算 

    for(i=; i<N; i++)

        {printf("%f ",FFT_RESULT(i));}/**/

    IFFT();//进行 IFFT计算

    for(i=; i<N; i++)

        {printf("%f ",IFFT_RESULT(i));}/**/

    return ;

}
/*******************************************************************************

** 程序名称：快速傅里叶变换(FFT)

** 程序描述：本程序实现快速傅里叶变换及其逆变换

** 性能提升：修正了FFT的bug,变量重新命名,并将 N_FFT改为动态值

** 程序版本：V6.6

** 程序作者：宋元瑞

** 最后修改：2011年5月16日

*******************************************************************************/

#include <stdio.h>

#include <stdlib.h>

#include <math.h>

//#define OUTPRINT printf

//#define DEL /##/

#define RESULT(x) sqrt(data_of_N_FFT[x].real*data_of_N_FFT[x].real+data_of_N_FFT[x].imag*data_of_N_FFT[x].imag)

#define PI  3.14159265358979323846264338327950288419716939937510    //圆周率，50位小数

#define PI2 6.28318530717958647692528676655900576839433879875021

int N_FFT=;                //傅里叶变换的点数

int M_of_N_FFT=;            //蝶形运算的级数，N = 2^M

int Npart2_of_N_FFT=;        //创建正弦函数表时取PI的1/2

int Npart4_of_N_FFT=;        //创建正弦函数表时取PI的1/4

typedef float ElemType;        //原始数据序列的数据类型,可以在这里设置

typedef struct                //定义复数结构体

{

    ElemType real,imag;

}complex_of_N_FFT,*ptr_complex_of_N_FFT;

ptr_complex_of_N_FFT data_of_N_FFT=NULL;//开辟存储单元，原始数据与负数结果均使用之

ElemType* SIN_TABLE_of_N_FFT=NULL;

//产生模拟原始数据输入

void InputData(void)

{

    int i;

    for (i=; i<N_FFT; i++)//制造输入序列

    {

        data_of_N_FFT[i].real = sin(*PI*i/N_FFT);    //正弦波

        data_of_N_FFT[i].imag = ;

        printf("%f ",data_of_N_FFT[i].real);

    }

}

//创建正弦函数表

void CREATE_SIN_TABLE(void)

{

    int i=;

    for (i=; i<=Npart4_of_N_FFT; i++)

    {

        SIN_TABLE_of_N_FFT[i] = sin(PI*i/Npart2_of_N_FFT);//SIN_TABLE[i] = sin(PI2*i/N);

    }

}

ElemType Sin_find(ElemType x)

{

    int i = (int)(N_FFT*x);

    i = i>>;

    if (i>Npart4_of_N_FFT)//注意：i已经转化为0~N之间的整数了！

    {

        //不会超过N/2

        i = Npart2_of_N_FFT - i;//i = i - 2*(i-Npart4);

    }

    return SIN_TABLE_of_N_FFT[i];

}

ElemType Cos_find(ElemType x)

{

    int i = (int)(N_FFT*x);

    i = i>>;

    if (i<Npart4_of_N_FFT)//注意：i已经转化为0~N之间的整数了！

    {

        //不会超过N/2

        //i = Npart4 - i;

        return SIN_TABLE_of_N_FFT[Npart4_of_N_FFT - i];

    }

    else//i>Npart4 && i<N/2

    {

        //i = i - Npart4;

        return -SIN_TABLE_of_N_FFT[i - Npart4_of_N_FFT];

    }

}

//变址

void ChangeSeat(complex_of_N_FFT *DataInput)

{

    int nextValue,nextM,i,k,j=;

    complex_of_N_FFT temp;

    nextValue=N_FFT/;                  //变址运算，即把自然顺序变成倒位序，采用雷德算法

    nextM=N_FFT-;

    for (i=;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/;                    //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运算函数

void FFT(void)

{

    int L=,B=,J=,K=;

    int step=, KB=;

    //ElemType P=0;

    ElemType angle;

    complex_of_N_FFT W,Temp_XX;

    ChangeSeat(data_of_N_FFT);//变址

    //CREATE_SIN_TABLE();

    for (L=; L<=M_of_N_FFT; L++)

    {

        step = <<L;//2^L

        B = step>>;//B=2^(L-1)

        for (J=; 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_FFT; K=K+step)

            {

                KB = K + B;

                //Temp_XX = XX_complex(data[KB],W);

                //用下面两行直接计算复数乘法，省去函数调用开销

                Temp_XX.real = data_of_N_FFT[KB].real * W.real-data_of_N_FFT[KB].imag*W.imag;

                Temp_XX.imag = W.imag*data_of_N_FFT[KB].real + data_of_N_FFT[KB].imag*W.real;

                data_of_N_FFT[KB].real = data_of_N_FFT[K].real - Temp_XX.real;

                data_of_N_FFT[KB].imag = data_of_N_FFT[K].imag - Temp_XX.imag;

                data_of_N_FFT[K].real = data_of_N_FFT[K].real + Temp_XX.real;

                data_of_N_FFT[K].imag = data_of_N_FFT[K].imag + Temp_XX.imag;

            }

        }

    }

}

//IFFT运算函数

void IFFT(void)

{

    int L=,B=,J=,K=;

    int step=, KB=;

    //ElemType P=0;

    ElemType angle;

    complex_of_N_FFT W,Temp_XX;

    ChangeSeat(data_of_N_FFT);//变址

    //CREATE_SIN_TABLE();

    for (L=; L<=M_of_N_FFT; L++)

    {

        step = <<L;//2^L

        B = step>>;//B=2^(L-1)

        for (J=; 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_FFT; K=K+step)

            {

                KB = K + B;

                //Temp_XX = XX_complex(data[KB],W);

                //用下面两行直接计算复数乘法，省去函数调用开销

                Temp_XX.real = data_of_N_FFT[KB].real * W.real-data_of_N_FFT[KB].imag*W.imag;

                Temp_XX.imag = W.imag*data_of_N_FFT[KB].real + data_of_N_FFT[KB].imag*W.real;

                data_of_N_FFT[KB].real = data_of_N_FFT[K].real - Temp_XX.real;

                data_of_N_FFT[KB].imag = data_of_N_FFT[K].imag - Temp_XX.imag;

                data_of_N_FFT[K].real = data_of_N_FFT[K].real + Temp_XX.real;

                data_of_N_FFT[K].imag = data_of_N_FFT[K].imag + Temp_XX.imag;

            }

        }

    }

}

//初始化FFT程序

//N_FFT是 FFT的点数，必须是2的次方

void Init_FFT(int N_of_FFT)

{

    int i=;

    int temp_N_FFT=;

    N_FFT = N_of_FFT;                    //傅里叶变换的点数 ，必须是 2的次方

    M_of_N_FFT = ;                    //蝶形运算的级数，N = 2^M

    for (i=; temp_N_FFT<N_FFT; i++)

    {

        temp_N_FFT = *temp_N_FFT;

        M_of_N_FFT++;

    }

    printf("\n%d\n",M_of_N_FFT);

    Npart2_of_N_FFT = N_FFT/;        //创建正弦函数表时取PI的1/2

    Npart4_of_N_FFT = N_FFT/;        //创建正弦函数表时取PI的1/4

    //ptr_complex_of_N_FFT data_of_N_FFT=NULL;//开辟存储单元，原始数据与负数结果均使用之

    data_of_N_FFT = (ptr_complex_of_N_FFT)malloc(N_FFT * sizeof(complex_of_N_FFT));

    //ptr_complex_of_N_FFT SIN_TABLE_of_N_FFT=NULL;

    SIN_TABLE_of_N_FFT = (ElemType *)malloc((Npart4_of_N_FFT+) * sizeof(ElemType));

    CREATE_SIN_TABLE();                //创建正弦函数表

}

//结束 FFT运算，释放相关内存

void Close_FFT(void)

{

    if (data_of_N_FFT != NULL)

    {

        free(data_of_N_FFT);          //释放内存

        data_of_N_FFT = NULL;

    }

    if (SIN_TABLE_of_N_FFT != NULL)

    {

        free(SIN_TABLE_of_N_FFT);     //释放内存

        SIN_TABLE_of_N_FFT = NULL;

    }

}

//主函数

int main(int argc, char *argv[])

{

    int i = ;

    Init_FFT();        //初始化各项参数，此函数只需执行一次

    InputData();        //输入原始数据

    FFT();                //进行 FFT计算

    printf("\n\n");

    for (i=; i<N_FFT; i++)

    {

        printf("%f ",RESULT(i));

    }

    IFFT();//进行 IFFT计算

    printf("\n\n");

    for (i=; i<N_FFT; i++)

    {

        printf("%f ",data_of_N_FFT[i].real/N_FFT);

    }

    Close_FFT();        //结束 FFT运算，释放相关内存

    return ;

}
/*******************************************************************************

** 模块名称：快速傅里叶变换(FFT)

** 模块描述：本程序实现快速傅里叶变换

** 性能提升：已达到网上同类程序最高速度

** 模块版本：V6.0

** 模块作者：宋元瑞

** 最后修改：2011年5月6日

*******************************************************************************/

/*******************************************************************************

**    程序说明：

    FFT运算输入参数 source_Data(i) 是一个N大小的数组（注意是小括号）

    FFT运算输出结果 result_Data(i) 是一个N大小的数组（注意是小括号）

**    调用举例：

int main(int argc, char *argv[])

{

    int i = 0;

    ///以下为FFT运算的调用方式

    FFT_prepare();        //为FFT做好准备，此函数只需程序开始时执行一次即可，切勿写在循环中

    while(1)

    {

        for(i=0;i<N_FFT;i++)    //输入原始数据

            {source_Data(i) = sin(2*PI*i/N_FFT);}//注意inputData后面是小括号

        FFT();                //进行FFT计算

        //输出结果：XXX =  result_Data(i);

    }

    return 0;

}

*******************************************************************************/

#ifndef SYRFFT_H

//#pragma once

//#include <stdio.h>

#include <math.h>

#define FFT_prepare() CREATE_SIN_TABLE();    //创建正弦函数表

#define source_Data(i) data_FFT[i].imag        //接收输入数据的数组，大小为N

#define result_Data(i) sqrt(data_FFT[i].real*data_FFT[i].real+data_FFT[i].imag*data_FFT[i].imag)//FFT结果 

#define PI  3.14159265358979323846264338327950288419716939937510    //圆周率，50位小数

#define PI2 6.28318530717958647692528676655900576839433879875021

#define N_FFT 1024                //傅里叶变换的点数

#define M_of_N_FFT 10                //蝶形运算的级数，N = 2^M

#define Npart2_of_N_FFT 512            //创建正弦函数表时取PI的1/2

#define Npart4_of_N_FFT 256            //创建正弦函数表时取PI的1/4 

typedef float ElemType;        //原始数据序列的数据类型,可以在这里设置

typedef struct                //定义复数结构体

{

    ElemType real,imag;

}complex_of_N_FFT;

complex_of_N_FFT data_FFT[N_FFT];        //存放要进行FFT运输的数据，运算结果也存放这里

ElemType SIN_TABLE[Npart4_of_N_FFT+];

//产生模拟原始数据输入

/*

void InputData(void)

{

    int i;

    for(i=0; i<N_FFT; i++)        //制造输入序列

    {

        source_Data(i) = sin(2*PI*i/N_FFT);    //模拟输入正弦波

        //data_FFT[i].imag = sin(2*PI*i/N);    //正弦波

        //printf("%f ",data_FFT[i].imag);

    }

}*/

//创建正弦函数表

void CREATE_SIN_TABLE(void)

{

    int i=;

    for(i=; i<=Npart4_of_N_FFT; i++)

    {

        SIN_TABLE[i] = sin(PI*i/Npart2_of_N_FFT);//SIN_TABLE[i] = sin(PI2*i/N);

    }

}

ElemType Sin_find(ElemType x)

{

    int i = (int)(N_FFT*x);

    i = i>>;

    if(i>Npart4_of_N_FFT)//注意：i已经转化为0~N之间的整数了！

    {//不会超过N/2

        i = Npart2_of_N_FFT - i;//i = i - 2*(i-Npart4);

    }

    return SIN_TABLE[i];

}

ElemType Cos_find(ElemType x)

{

    int i = (int)(N_FFT*x);

    i = i>>;

    if(i<Npart4_of_N_FFT)//注意：i已经转化为0~N之间的整数了！

    {//不会超过N/2

        //i = Npart4 - i;

        return SIN_TABLE[Npart4_of_N_FFT - i];

    }

    else//i>Npart4 && i<N/2

    {

        //i = i - Npart4;

        return -SIN_TABLE[i - Npart4_of_N_FFT];

    }

}

//变址

void ChangeSeat(complex_of_N_FFT *DataInput)

{

    //变址前data_FFT[i].imag存储了输入数据，变址后data_FFT[i].real存储了输入数据

    int i,n,New_seat;

    for(i=; i<N_FFT; i++)

    {

        New_seat = ;

        for(n=;n<M_of_N_FFT;n++)

        {

            New_seat = New_seat | (((i>>n) & 0x01) << (M_of_N_FFT-n-));

        }

        DataInput[New_seat].real = DataInput[i].imag;

    }

    for(i=; i<N_FFT; i++)

    {

        DataInput[i].imag = ;

    }

}                                                      

//复数乘法

/*

complex_of_N_FFT XX_complex(complex_of_N_FFT a, complex_of_N_FFT b)

{

    complex_of_N_FFT 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=,B=,J=,K=;

    int step=, KB=;

    //ElemType P=0;

    ElemType angle;

    complex_of_N_FFT W,Temp_XX;

    ChangeSeat(data_FFT);    //变址 

    for(L=; L<=M_of_N_FFT; L++)

    {

        step = <<L;//2^L

        B = step>>;//B=2^(L-1)

        for(J=; 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_FFT; K=K+step)

            {

                KB = K + B;

                //Temp_XX = XX_complex(data_FFT[KB],W);

                //用下面两行直接计算复数乘法，省去函数调用开销

                Temp_XX.real = data_FFT[KB].real * W.real-data_FFT[KB].imag*W.imag;

                Temp_XX.imag = W.imag*data_FFT[KB].real + data_FFT[KB].imag*W.real;

                data_FFT[KB].real = data_FFT[K].real - Temp_XX.real;

                data_FFT[KB].imag = data_FFT[K].imag - Temp_XX.imag;

                data_FFT[K].real = data_FFT[K].real + Temp_XX.real;

                data_FFT[K].imag = data_FFT[K].imag + Temp_XX.imag;

            }

        }

    }

}

#define SYRFFT_H

#endif
/*******************************************************************************

** 程序名称：快速傅里叶变换(FFT)

** 程序描述：本程序实现快速傅里叶变换

** 程序版本：V6.0

** 程序作者：宋元瑞

** 最后修改：2011年5月6日

*******************************************************************************/

#include <stdio.h>

#include "syrFFT_H.h"    //包含FFT运算头文件 

//主函数

int main(int argc, char *argv[])

{

    int i = ;

    //以下3句为FFT运算的调用函数

    FFT_prepare();        //为FFT做好准备，此函数只需程序开始时执行一次即可，切勿写在循环中

    //while(1)

    {

        for(i=;i<N_FFT;i++)//模拟输入

            {source_Data(i) = sin(*PI*i/N_FFT);}//注意inputData后面是小括号

        FFT();

        for(i=; i<N_FFT; i++)//输出结果

            {printf("%lf ",result_Data(i));}

    } 

    return ;

}
#ifndef FFT_H

#include <stdlib.h>

#include <math.h>

#define RESULT(x) sqrt(data_of_N_FFT[x].real*data_of_N_FFT[x].real+data_of_N_FFT[x].imag*data_of_N_FFT[x].imag)

#define PI  3.14159265358979323846264338327950288419716939937510    //圆周率，50位小数

#define PI2 6.28318530717958647692528676655900576839433879875021

int N_FFT=;                //傅里叶变换的点数

int M_of_N_FFT=;            //蝶形运算的级数，N = 2^M

int Npart2_of_N_FFT=;        //创建正弦函数表时取PI的1/2

int Npart4_of_N_FFT=;        //创建正弦函数表时取PI的1/4 

typedef float ElemType;        //原始数据序列的数据类型,可以在这里设置

typedef struct                //定义复数结构体

{

    ElemType real,imag;

}complex_of_N_FFT,*ptr_complex_of_N_FFT;

ptr_complex_of_N_FFT data_of_N_FFT=NULL;//开辟存储单元，原始数据与负数结果均使用之

ElemType* SIN_TABLE_of_N_FFT=NULL;

//创建正弦函数表

void CREATE_SIN_TABLE(void)

{

    int i=;

    for(i=; i<=Npart4_of_N_FFT; i++)

    {

        SIN_TABLE_of_N_FFT[i] = sin(PI*i/Npart2_of_N_FFT);//SIN_TABLE[i] = sin(PI2*i/N);

    }

}

ElemType Sin_find(ElemType x)

{

    int i = (int)(N_FFT*x);

    i = i>>;

    if(i>Npart4_of_N_FFT)//注意：i已经转化为0~N之间的整数了！

    {//不会超过N/2

        i = Npart2_of_N_FFT - i;//i = i - 2*(i-Npart4);

    }

    return SIN_TABLE_of_N_FFT[i];

}

ElemType Cos_find(ElemType x)

{

    int i = (int)(N_FFT*x);

    i = i>>;

    if(i<Npart4_of_N_FFT)//注意：i已经转化为0~N之间的整数了！

    {//不会超过N/2

        //i = Npart4 - i;

        return SIN_TABLE_of_N_FFT[Npart4_of_N_FFT - i];

    }

    else//i>Npart4 && i<N/2

    {

        //i = i - Npart4;

        return -SIN_TABLE_of_N_FFT[i - Npart4_of_N_FFT];

    }

}

//变址

void ChangeSeat(complex_of_N_FFT *DataInput)

{

    int nextValue,nextM,i,k,j=;

    complex_of_N_FFT temp;

    nextValue=N_FFT/;                  //变址运算，即把自然顺序变成倒位序，采用雷德算法

    nextM=N_FFT-;

    for (i=;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/;                    //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运算函数

void FFT(void)

{

    int L=,B=,J=,K=;

    int step=, KB=;

    //ElemType P=0;

    ElemType angle;

    complex_of_N_FFT W,Temp_XX;

    ChangeSeat(data_of_N_FFT);//变址

    //CREATE_SIN_TABLE();

    for(L=; L<=M_of_N_FFT; L++)

    {

        step = <<L;//2^L

        B = step>>;//B=2^(L-1)

        for(J=; 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_FFT; K=K+step)

            {

                KB = K + B;

                //Temp_XX = XX_complex(data[KB],W);

                //用下面两行直接计算复数乘法，省去函数调用开销

                Temp_XX.real = data_of_N_FFT[KB].real * W.real-data_of_N_FFT[KB].imag*W.imag;

                Temp_XX.imag = W.imag*data_of_N_FFT[KB].real + data_of_N_FFT[KB].imag*W.real;

                data_of_N_FFT[KB].real = data_of_N_FFT[K].real - Temp_XX.real;

                data_of_N_FFT[KB].imag = data_of_N_FFT[K].imag - Temp_XX.imag;

                data_of_N_FFT[K].real = data_of_N_FFT[K].real + Temp_XX.real;

                data_of_N_FFT[K].imag = data_of_N_FFT[K].imag + Temp_XX.imag;

            }

        }

    }

}

//IFFT运算函数

void IFFT(void)

{

    int L=,B=,J=,K=;

    int step=, KB=;

    //ElemType P=0;

    ElemType angle;

    complex_of_N_FFT W,Temp_XX;

    ChangeSeat(data_of_N_FFT);//变址

    //CREATE_SIN_TABLE();

    for(L=; L<=M_of_N_FFT; L++)

    {

        step = <<L;//2^L

        B = step>>;//B=2^(L-1)

        for(J=; 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_FFT; K=K+step)

            {

                KB = K + B;

                //Temp_XX = XX_complex(data[KB],W);

                //用下面两行直接计算复数乘法，省去函数调用开销

                Temp_XX.real = data_of_N_FFT[KB].real * W.real-data_of_N_FFT[KB].imag*W.imag;

                Temp_XX.imag = W.imag*data_of_N_FFT[KB].real + data_of_N_FFT[KB].imag*W.real;

                data_of_N_FFT[KB].real = data_of_N_FFT[K].real - Temp_XX.real;

                data_of_N_FFT[KB].imag = data_of_N_FFT[K].imag - Temp_XX.imag;

                data_of_N_FFT[K].real = data_of_N_FFT[K].real + Temp_XX.real;

                data_of_N_FFT[K].imag = data_of_N_FFT[K].imag + Temp_XX.imag;

            }

        }

    }

}

//初始化FFT程序

//N_FFT是 FFT的点数，必须是2的次方

void Init_FFT(int N_of_FFT)

{

    int i=;

    int temp_N_FFT=;

    N_FFT = N_of_FFT;                    //傅里叶变换的点数 ，必须是 2的次方

    M_of_N_FFT = ;                    //蝶形运算的级数，N = 2^M

    for(i=; temp_N_FFT<N_FFT; i++)

    {

        temp_N_FFT = *temp_N_FFT;

        M_of_N_FFT++;

    }

    //printf("\n%d\n",M_of_N_FFT);

    Npart2_of_N_FFT = N_FFT/;        //创建正弦函数表时取PI的1/2

    Npart4_of_N_FFT = N_FFT/;        //创建正弦函数表时取PI的1/4 

    //ptr_complex_of_N_FFT data_of_N_FFT=NULL;//开辟存储单元，原始数据与负数结果均使用之

    data_of_N_FFT = (ptr_complex_of_N_FFT)malloc(N_FFT * sizeof(complex_of_N_FFT));

    //ptr_complex_of_N_FFT SIN_TABLE_of_N_FFT=NULL;

    SIN_TABLE_of_N_FFT = (ElemType *)malloc((Npart4_of_N_FFT+) * sizeof(ElemType));

    CREATE_SIN_TABLE();                //创建正弦函数表

}

//结束 FFT运算，释放相关内存

void Close_FFT(void)

{

    if(data_of_N_FFT != NULL)

    {

        free(data_of_N_FFT);          //释放内存

        data_of_N_FFT = NULL;

    }

    if(SIN_TABLE_of_N_FFT != NULL)

    {

        free(SIN_TABLE_of_N_FFT);     //释放内存

        SIN_TABLE_of_N_FFT = NULL;

    }

}

#define FFT_H

#endif
/*******************************************************************************

** 程序名称：快速傅里叶变换(FFT)

** 程序描述：本程序实现快速傅里叶变换及其逆变换

** 性能提升：修正了FFT的bug

** 程序版本：V6.6

** 程序作者：宋元瑞

** 最后修改：2011年5月16日

*******************************************************************************/

#include "jxust_fft6_6.h"

//产生模拟原始数据输入 ,在实际应用时替换为AD采样数据

void InputData(void)

{

    int i;

    for(i=; i<N_FFT; i++)//制造输入序列

    {

        data_of_N_FFT[i].real = sin(*PI*i/N_FFT);    //输入采样数据

        data_of_N_FFT[i].imag = ;

    }

}

//主函数 ，示例如何调用FFT运算

int main(int argc, char *argv[])

{

    int i = ;

    Init_FFT();        //①初始化各项参数，此函数只需执行一次 ; 参数为FFT的点数，必须为2的次方 

    InputData();        //②输入原始数据 ,此处在实际应用时替换为AD采样数据

    FFT();                //③进行 FFT计算

    //for(i=0; i<N_FFT; i++)//④输出FFT频谱结果  sqrt(data_of_N_FFT[i].real*data_of_N_FFT[i].real+data_of_N_FFT[i].imag*data_of_N_FFT[i].imag)

        //{printf("%f ",RESULT(i));}

    IFFT();//进行 IFFT计算

    //for(i=0; i<N_FFT; i++)//(可选步骤)⑤输出 IFFT结果 ，滤波时会用到；暂时不用

        //{printf("%f ",data_of_N_FFT[i].real/N_FFT);}

    Close_FFT();    //结束 FFT运算，释放相关内存 ;此函数在彻底结束FFT运算后调用， 

    return ;

}
#ifndef    SYRFFT_6_55_H

#include <math.h>

#define FFT_RESULT(x)     (sqrt(data_of_N_FFT[x].real*data_of_N_FFT[x].real+data_of_N_FFT[x].imag*data_of_N_FFT[x].imag))

#define IFFT_RESULT(x)    (data_of_N_FFT[x].real/N_FFT)

#define PI  3.14159265358979323846264338327950288419716939937510    //圆周率，50位小数

#define PI2 6.28318530717958647692528676655900576839433879875021

#define N_FFT            1024        //傅里叶变换的点数

#define M_of_N_FFT        10            //蝶形运算的级数，N = 2^M

#define Npart2_of_N_FFT    512            //创建正弦函数表时取PI的1/2

#define Npart4_of_N_FFT    256            //创建正弦函数表时取PI的1/4 

typedef float ElemType;        //原始数据序列的数据类型,可以在这里设置

typedef struct                //定义复数结构体

{

    ElemType real,imag;

}complex_of_N_FFT,*ptr_complex_of_N_FFT;

complex_of_N_FFT data_of_N_FFT[N_FFT];            //定义存储单元，原始数据与负数结果均使用之

ElemType SIN_TABLE_of_N_FFT[Npart4_of_N_FFT+];

//创建正弦函数表

void CREATE_SIN_TABLE(void)

{

    int i=;

    for(i=; i<=Npart4_of_N_FFT; i++)

    {

        SIN_TABLE_of_N_FFT[i] = sin(PI*i/Npart2_of_N_FFT);//SIN_TABLE[i] = sin(PI2*i/N);

    }

}

ElemType Sin_find(ElemType x)

{

    int i = (int)(N_FFT*x);//注意：i已经转化为0~N之间的整数了！

    i = i>>;// i = i / 2;

    if(i>Npart4_of_N_FFT)

    {//根据FFT相关公式，sin()参数不会超过PI， 即i不会超过N/2

        i = Npart2_of_N_FFT - i;//i = i - 2*(i-Npart4);

    }

    return SIN_TABLE_of_N_FFT[i];

}

ElemType Cos_find(ElemType x)

{

    int i = (int)(N_FFT*x);//注意：i已经转化为0~N之间的整数了！

    i = i>>;

    if(i < Npart4_of_N_FFT )

    { //不会超过N/2

        //i = Npart4 - i;

        return SIN_TABLE_of_N_FFT[Npart4_of_N_FFT - i];

    }

    else //i > Npart4 && i < N/2

    {

        //i = i - Npart4;

        return -SIN_TABLE_of_N_FFT[i - Npart4_of_N_FFT];

    }

}

//变址

void ChangeSeat(complex_of_N_FFT *DataInput)

{

    int nextValue,nextM,i,k,j=;

    complex_of_N_FFT temp;

    nextValue=N_FFT/;                  //变址运算，即把自然顺序变成倒位序，采用雷德算法

    nextM=N_FFT-;

    for (i=;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/;                    //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运算函数

void FFT(void)

{

    int L=,B=,J=,K=;

    int step=, KB=;

    //ElemType P=0;

    ElemType angle;

    complex_of_N_FFT W,Temp_XX;

    ChangeSeat(data_of_N_FFT);//变址

    //CREATE_SIN_TABLE();

    for(L=; L<=M_of_N_FFT; L++)

    {

        step = <<L;//2^L

        B = step>>;//B=2^(L-1)

        for(J=; 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_FFT; K=K+step)

            {

                KB = K + B;

                //Temp_XX = XX_complex(data[KB],W);

                //用下面两行直接计算复数乘法，省去函数调用开销

                Temp_XX.real = data_of_N_FFT[KB].real * W.real-data_of_N_FFT[KB].imag*W.imag;

                Temp_XX.imag = W.imag*data_of_N_FFT[KB].real + data_of_N_FFT[KB].imag*W.real;

                data_of_N_FFT[KB].real = data_of_N_FFT[K].real - Temp_XX.real;

                data_of_N_FFT[KB].imag = data_of_N_FFT[K].imag - Temp_XX.imag;

                data_of_N_FFT[K].real = data_of_N_FFT[K].real + Temp_XX.real;

                data_of_N_FFT[K].imag = data_of_N_FFT[K].imag + Temp_XX.imag;

            }

        }

    }

}

//IFFT运算函数

void IFFT(void)

{

    int L=,B=,J=,K=;

    int step=, KB=;

    //ElemType P=0;

    ElemType angle;

    complex_of_N_FFT W,Temp_XX;

    ChangeSeat(data_of_N_FFT);//变址

    //CREATE_SIN_TABLE();

    for(L=; L<=M_of_N_FFT; L++)

    {

        step = <<L;//2^L

        B = step>>;//B=2^(L-1)

        for(J=; 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_FFT; K=K+step)

            {

                KB = K + B;

                //Temp_XX = XX_complex(data[KB],W);

                //用下面两行直接计算复数乘法，省去函数调用开销

                Temp_XX.real = data_of_N_FFT[KB].real * W.real-data_of_N_FFT[KB].imag*W.imag;

                Temp_XX.imag = W.imag*data_of_N_FFT[KB].real + data_of_N_FFT[KB].imag*W.real;

                data_of_N_FFT[KB].real = data_of_N_FFT[K].real - Temp_XX.real;

                data_of_N_FFT[KB].imag = data_of_N_FFT[K].imag - Temp_XX.imag;

                data_of_N_FFT[K].real = data_of_N_FFT[K].real + Temp_XX.real;

                data_of_N_FFT[K].imag = data_of_N_FFT[K].imag + Temp_XX.imag;

            }

        }

    }

}

//初始化FFT程序

void Init_FFT()

{

    CREATE_SIN_TABLE();                //创建正弦函数表

}

//结束 FFT运算，释放相关内存

void Close_FFT(void)

{

}

#define SYRFFT_6_55_H

#endif
/*******************************************************************************

** 程序名称：快速傅里叶变换(FFT)

** 程序描述：本程序实现快速傅里叶变换

** 性能提升：修正了IFFT的bug,变量重新命名,使用宏定义改变N大小

** 程序版本：V6.55

** 程序作者：宋元瑞

** 最后修改：2011年5月22日

*******************************************************************************/

#include "syrFFT_6_55.h" 

//产生模拟原始数据输入 ,在实际应用时替换为AD采样数据

void InputData(void)

{

    int i;

    for(i=; i<N_FFT; i++)//制造输入序列

    {

        data_of_N_FFT[i].real = sin(*PI*i/N_FFT);    //输入采样数据

        data_of_N_FFT[i].imag = ;

    }

}

//主函数 ，示例如何调用FFT运算

int main(int argc, char *argv[])

{

    int i = ;

    Init_FFT();            //①初始化各项参数，此函数只需执行一次 ; 修改FFT的点数去头文件的宏定义处修改 

    InputData();        //②输入原始数据 ,此处在实际应用时替换为AD采样数据

    FFT();                //③进行 FFT计算

    //for(i=0; i<N_FFT; i++)//④输出FFT频谱结果  sqrt(data_of_N_FFT[i].real*data_of_N_FFT[i].real+data_of_N_FFT[i].imag*data_of_N_FFT[i].imag)

        //{printf("%f ",FFT_RESULT(i));}

    IFFT();//进行 IFFT计算

    //for(i=0; i<N_FFT; i++)//(可选步骤)⑤输出 IFFT结果 ，滤波时会用到；暂时不用

        //{printf("%f ",IFFT_RESULT(i));}

    Close_FFT();    //结束 FFT运算，此版本此句无用，可不写 

    return ;

}

回复更精彩，请到原出处关注。

相关文章