C语言FFT

/************FFT***********/



#include < stdio.h >

    #include < math.h >

    #include < stdlib.h >

    #define N 1000

typedef struct



{



    double real;



    double img;



}

complex;



void fft(); /*快速傅里叶变换*/



void ifft(); /*快速傅里叶逆变换*/



void initW();



void change();



void add(complex, complex, complex * ); /*复数加法*/



void mul(complex, complex, complex * ); /*复数乘法*/



void sub(complex, complex, complex * ); /*复数减法*/



void divi(complex, complex, complex * ); /*复数除法*/



void output(); /*输出结果*/



complex x[N], * W; /*输出序列的值*/



int size_x = 0; /*输入序列的长度，只限2的N次方*/



double PI;



int main()



{



    int i, method;



    system("cls");



    PI = atan(1) * 4;



    printf("Please input the size of x:\n");



    /*输入序列的长度*/



    scanf("%d", & size_x);



    printf("Please input the data in x[N]:(such as:5 6)\n");



    /*输入序列对应的值*/



    for (i = 0; i < size_x; i++)



        scanf("%lf %lf", & x[i].real, & x[i].img);



    initW();



    /*选择FFT或逆FFT运算*/



    printf("Use FFT(0) or IFFT(1)?\n");



    scanf("%d", & method);



    if (method == 0)



        fft();



    else



        ifft();



    output();



    return 0;



}



/*进行基-2 FFT运算*/



void fft()



{



    int i = 0, j = 0, k = 0, l = 0;



    complex up, down, product;



    change();



    for (i = 0; i < log(size_x) / log(2); i++)



    {



        l = 1 << i;



        for (j = 0; j < size_x; j += 2 * l)



        {



            for (k = 0; k < l; k++)



            {



                mul(x[j + k + l], W[size_x * k / 2 / l], & product);



                add(x[j + k], product, & up);



                sub(x[j + k], product, & down);



                x[j + k] = up;



                x[j + k + l] = down;



            }



        }



    }



}



void ifft()



{



    int i = 0, j = 0, k = 0, l = size_x;



    complex up, down;



    for (i = 0; i < (int)(log(size_x) / log(2)); i++) /*蝶形运算*/



    {



        l /= 2;



        for (j = 0; j < size_x; j += 2 * l)



        {



            for (k = 0; k < l; k++)



            {



                add(x[j + k], x[j + k + l], & up);



                up.real /= 2;

                up.img /= 2;



                sub(x[j + k], x[j + k + l], & down);



                down.real /= 2;

                down.img /= 2;



                divi(down, W[size_x * k / 2 / l], & down);



                x[j + k] = up;



                x[j + k + l] = down;



            }



        }



    }



    change();



}



void initW()



{



    int i;



    W = (complex * ) malloc(sizeof(complex) * size_x);



    for (i = 0; i < size_x; i++)



    {



        W[i].real = cos(2 * PI / size_x * i);



        W[i].img = -1 * sin(2 * PI / size_x * i);



    }



}



void change()



{



    complex temp;



    unsigned short i = 0, j = 0, k = 0;



    double t;



    for (i = 0; i < size_x; i++)



    {



        k = i;

        j = 0;



        t = (log(size_x) / log(2));



        while ((t--) > 0)



        {



            j = j << 1;



            j |= (k & 1);



            k = k >> 1;



        }



        if (j > i)



        {



            temp = x[i];



            x[i] = x[j];



            x[j] = temp;



        }



    }



}



void output() /*输出结果*/



{



    int i;



    printf("The result are as follows\n");



    for (i = 0; i < size_x; i++)



    {



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



        if (x[i].img >= 0.0001)



            printf("+%.4fj\n", x[i].img);



        else if (fabs(x[i].img) < 0.0001)



            printf("\n");



        else



            printf("%.4fj\n", x[i].img);



    }



}



void add(complex a, complex b, complex * c)



{



    c - > real = a.real + b.real;



    c - > img = a.img + b.img;



}



void mul(complex a, complex b, complex * c)



{



    c - > real = a.real * b.real - a.img * b.img;



    c - > img = a.real * b.img + a.img * b.real;



}



void sub(complex a, complex b, complex * c)



{



    c - > real = a.real - b.real;



    c - > img = a.img - b.img;



}



void divi(complex a, complex b, complex * c)



{



    c - > real = (a.real * b.real + a.img * b.img) / (



        b.real * b.real + b.img * b.img);



    c - > img = (a.img * b.real - a.real * b.img) / (b.real * b.real + b.img * b.img);



}