[N32G45x] 快速傅里叶变换

#include <stdio.h>

#include <math.h>

#define PI 3.1415926

#define FFT_LENGTH                8

double input[FFT_LENGTH]={1,1,1,1,1,1,1,1};

struct complex1{                //定义一个复数结构体

        double real;        //实部

        double image;        //虚部

};        

//将input的实数结果存放为复数

struct complex1 result_dat[8];

/*

        虚数的乘法

*/

struct complex1 con_complex(struct complex1 a,struct complex1 b){

        struct complex1 temp;

        temp.real=(a.real*b.real)-(a.image*b.image);

        temp.image=(a.image*b.real)+(a.real*b.image);

        return temp;

}

/*

        简单的a的b次方

*/

int mypow(int a,int b){

        int i,sum=a;

        if(b==0)return 1;

        for(i=1;i<b;i++){

                sum*=a;        

        }

        return sum;

}

/*

        简单的求以2为底的正整数

*/

int log2(int n){

        unsigned i=1;

        int sum=1;

        for(i;;i++){

                sum*=2;

                if(sum>=n)break;

        }

        return i;

}

/*

        简单的交换数据的函数

*/

void swap(struct complex1 *a,struct complex1 *b){

        struct complex1 temp;

    temp=*a;

    *a=*b;

    *b=temp;

}

/*

        dat为输入数据的数组

        N为抽样次数  也代表周期  必须是2^N次方

*/

void fft(struct complex1 dat[],unsigned char N){        

        /*最终  dat_buf计算出 当前蝶形运算奇数项与W  乘积

                        dat_org存放上一个偶数项的值

        */

        struct complex1 dat_buf,dat_org;

        /*        L为几级蝶形运算    也代表了2进制的位数

                n为当前级蝶形的需要次数  n最初为N/2 每级蝶形运算后都要/2

                i j为倒位时要用到的自增符号  同时  i也用到了L碟级数   j是计算当前碟级的计算次数

                re_i i_copy均是倒位时用到的变量

                k为当前碟级  cos(2*pi/N*k)的  k   也是e^(-j2*pi/N)*k  的  k

        */

        unsigned char L,i,j,re_i=0,i_copy=0,k=0,fft_flag=1;

        //经过观察，发现每级蝶形运算需要N/2次运算，共运算N/2*log2N  次  

        unsigned char fft_counter=0;

        //在此要进行补2   N必须是2^n   在此略

        //蝶形级数  (L级)

        L=log2(N);        

        //计算每级蝶形计算的次数(这里只是一个初始值)  之后每次要/2

        //n=N/2;

        //对dat的顺序进行倒位

        for(i=1;i<N/2;i++){

                i_copy=i;

                re_i=0;

                for(j=L-1;j>0;j--){

                        //判断i的副本最低位的数字  并且移动到最高位  次高位  ..

                        //re_i为交换的数   每次它的数字是不能移动的 并且循环之后要清0

                        re_i|=((i_copy&0x01)<<j);               

                        i_copy>>=1;

                }

                swap(&dat[i],&dat[re_i]);

        }

        //进行fft计算

        for(i=0;i<L;i++){

               

                fft_flag=1;

                fft_counter=0;

                for(j=0;j<N;j++){

                        if(fft_counter==mypow(2,i)){                //控制隔几次，运算几次,

                                fft_flag=0;

                        }else if(fft_counter==0){                //休止结束，继续运算

                                fft_flag=1;

                        }

                        //当不判断这个语句的时候  fft_flag保持  这样就可以持续运算了

                        if(fft_flag){

                                dat_buf.real=cos((2*PI*k)/(N/mypow(2,L-i-1)));

                                dat_buf.image=-sin((2*PI*k)/(N/mypow(2,L-i-1)));

                                dat_buf=con_complex(dat[j+mypow(2,i)],dat_buf);

                                                                                                //计算 当前蝶形运算奇数项与W  乘积

                                dat_org.real=dat[j].real;

                                dat_org.image=dat[j].image;                //暂存

                                dat[j].real=dat_org.real+dat_buf.real;

                                dat[j].image=dat_org.image+dat_buf.image;               

                                                                                                        //实部加实部   虚部加虚部

                                dat[j+mypow(2,i)].real=dat_org.real-dat_buf.real;

                                dat[j+mypow(2,i)].image=dat_org.image-dat_buf.image;

                                                                                                        //实部减实部        虚部减虚部

                                k++;

                                fft_counter++;

                        }else{

                                fft_counter--;                                //运算几次，就休止几次

                                k=0;

                        }

                }

        }

        

}

void main(){

        int i;

        //先将输入信号转换成复数

        for(i=0;i<FFT_LENGTH;i++){

                result_dat[i].image=0;        

                                                                //输入信号是二维的，暂时不存在复数

                result_dat[i].real=input[i];

                //result_dat[i].real=10;

                                                                //输入信号都为实数

        }

        fft(result_dat,FFT_LENGTH);

        for(i=0;i<FFT_LENGTH;i++){

                input[i]=sqrt(result_dat[i].real*result_dat[i].real+result_dat[i].image*result_dat[i].image);

                //取模

                printf("%lf\n",input[i]);

        }

        while(1);

}