c语言做快速傅里叶变换和快速逆傅里叶变换

C语言做快速傅里叶变换和快速逆傅里叶变换

  快速傅里叶变换（FFT）和快速逆傅里叶变换(IFFT)要求做傅里叶变换的数据点数只能是2的整数次幂，比如2，4，8，16，32，64，128，256，512，1024，2048，.......如果是2000个数据，那么用快速傅里叶变换（FFT）的结果就不对了，就需要使用对数据点数不限制的离散傅里叶变化（DFT）了，据说也可以使用补0的方法凑够2048个点来使用FFT。
   这里只分享C语言做快速傅里叶变化的代码，以后再分享C语言做傅里叶变换的代码。C语言做快速傅里叶变化（FFT）的时间复杂度是O(nlogn),傅里叶变化（DFT）的时间复杂度是O(n*n).
以下是做快速傅里叶变换和快速逆傅里叶变换的代码（包含测试代码）：

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

//定义复数结构体
typedef struct
{
    double real;
    double img;
}complex;

void initW(complex W[], int size_x, double PI); //初始化变化核
void change(int size_x, complex x[]);
void add(complex a, complex b, complex *c);
void mul(complex a, complex b, complex *c);
void sub(complex a, complex b, complex *c);
void divi(complex a, complex b, complex *c);
void output(int size_x, complex x[]);  //输出结果
void fft(int size_x, complex x[], complex W[]);
void ifft(int size_x, complex x[], complex W[]);

#define    N     1000

int main()
{
    int size_x = 0;     //输入序列的长度，只限2的N次方
    complex x[N], W[N]; //输出序列的值
    double PI;
    PI = atan(1) * 4;                   //pi等于4乘以1.0的正切值

    size_x = 16;
//    x[0].real = 136, x[0].img=0;
//    x[1].real = -8, x[1].img=40.2187;
//    x[2].real = -8, x[2].img=19.3137;
//    x[3].real = -8, x[3].img=11.9728;
//    x[4].real = -8, x[4].img=8;
//    x[5].real = -8, x[5].img=5.3454;
//    x[6].real = -8, x[6].img=3.3137;
//    x[7].real = -8, x[7].img=1.5913;
//    x[8].real = -8, x[8].img=0;
//    x[9].real = -8, x[9].img=-1.5913;
//    x[10].real = -8, x[10].img=-3.3137;
//    x[11].real = -8 x[11].img=-5.3454;
//    x[12].real = -8, x[12].img=-8;
//    x[13].real = -8, x[13].img=-11.9728;
//    x[14].real = -8, x[14].img=-19.3137;
//    x[15].real = -8, x[15].img=-40.2187;

    x[0].real = 1,   x[0].img=0;
    x[1].real = 2,   x[1].img=0;
    x[2].real = 3,   x[2].img=0;
    x[3].real = 4,   x[3].img=0;
    x[4].real = 5,   x[4].img=0;
    x[5].real = 6,   x[5].img=0;
    x[6].real = 7,   x[6].img=0;
    x[7].real = 8,   x[7].img=0;
    x[8].real = 9,   x[8].img=0;
    x[9].real = 10,  x[9].img=0;
    x[10].real = 11, x[10].img=0;
    x[11].real = 12, x[11].img=0;
    x[12].real = 13, x[12].img=0;
    x[13].real = 14, x[13].img=0;
    x[14].real = 15, x[14].img=0;
    x[15].real = 16, x[15].img=0;

    initW(W, size_x, PI);
    fft(size_x, x, W);         //做傅里叶变换
    //ifft(size_x, x, W);     //做逆傅里叶变换
    output(size_x, x);
    return 0;
}

void initW(complex W[], int size_x, double PI)
{
    for(int 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(int size_x, complex x[])
{
    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 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);
}

void output(int size_x, complex x[])
{
    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 fft(int size_x, complex x[], complex W[])
{
    int i = 0, j = 0, k = 0, l = 0;
    complex up, down, product;
    change(size_x, x);

    for(i = 0; i < (int)(log10(size_x)/log10(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 size_x, complex x[], complex W[])
{
    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(size_x, x);
}

 分享一个在线做傅里叶变换的网站，可以验证自己做傅里叶变换的结果。http://www.yunsuan.info/signalprocessing/compute_fft_2d.htm