C语言FFT算法

matlab测试结果

FFT.c

#include "fft.h"

void swap(Complex *a, Complex *b)
{
    Complex temp = *a;
    *a = *b;
    *b = temp;
}

void change(Complex y[], int len)
{
    int i, j, k;
    for (i = 1, j = len / 2; i < len - 1; i++)
    {
        if (i < j)
            swap(&y[i], &y[j]);

        k = len / 2;

        while (j >= k)
        {
            j = j - k;
            k = k / 2;
        }

        if (j < k)
            j += k;
    }
}

void fft(Complex y[], int len, int on)
{
    change(y, len);
    for (int h = 2; h <= len; h <<= 1)
    {
        Complex wn = {cos(2 * PI / h), -sin(on * 2 * PI / h)};

        for (int j = 0; j < len; j += h)
        {
            Complex w = {1, 0};

            for (int k = j; k < j + h / 2; k++)
            {
                /*Complex u=y[k];
                Complex t=w*y[k+h/2];
                y[k]=u+t;
                y[k+h/2]=u-t;
                w=w*wn;*/

                Complex u = y[k];
                Complex t;
                t.real = w.real * y[k + h / 2].real - w.imag * y[k + h / 2].imag;
                t.imag = w.real * y[k + h / 2].imag + w.imag * y[k + h / 2].real;

                y[k].real = u.real + t.real;
                y[k].imag = u.imag + t.imag;

                y[k + h / 2].real = u.real - t.real;
                y[k + h / 2].imag = u.imag - t.imag;

                Complex temp;
                temp.real = w.real * wn.real - w.imag * wn.imag;
                temp.imag = w.real * wn.imag + w.imag * wn.real;
                w.real = temp.real;
                w.imag = temp.imag;
            }
        }
    }

    if (on == -1)
    {
        for (int i = 0; i < len; i++)
        {
            y[i].real /= len;
            y[i].imag /= len;
        }
    }
}

FFT.h

#ifndef __FFT_H
#define __FFT_H

#include <math.h>

//#define PI acos(-1.0)
#define PI 3.1415926

typedef struct  {
    float real;
    float imag;
} Complex;

void swap(Complex *a, Complex *b);
void change(Complex y[],int len);
void fft(Complex y[],int len,int on);

#endif