傅里叶变换

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

typedef char int8_t;
typedef unsigned char uint8_t;
typedef unsigned short uint16_t;
typedef short int16_t;

#define N    64    //傅里叶变换的点数
#define M     6    //蝶形运算的级数，N = 2^M
#define N2    32    //N/2
#define N4    16    //N/4

#define PI  3.14159    //圆周率
#define PI2 6.28318

//定义复数结构体
typedef struct                
{
    float real;    //实部
    float imag;    //虚部
}complex;

//傅里叶变换序列，一维
float pr[64] = 
{
    1234,1285,1151,1086,896,671,541,392,368,565,762,905,987,1103,1352,1601,1593,1565,1576,1565,1379,1152,1208,1150,1086,1124,1092,945,791,561,393,291,352,596,950,1307,1490,1707,1760,1494,1351,1510,1427,1191,1156,1090,953,917,944,847,655,580,457,495,679,877,1015,1084,1256,1519,1469,1411,1423,1509
};

//正弦函数表
float sin_table[N4 + 1];


float find_sin(float x)
{
    int i = ((int)(N * x)) >> 1;

    if (i > N4)        //不会超过N/2 
        i = N2 - i;

    return sin_table[i];
}

float find_cos(float x)
{
    int i = ((int)(N * x)) >> 1;

    if (i < N4)
        return sin_table[N4 - i];
    else                //i>Npart4 && i<N/2
        return -sin_table[i - N4];
}

//dir:1 - 傅里叶变换， -1 - 傅里叶逆变换
void fft(complex data[], uint8_t num, uint8_t n, int8_t dir)
{
    uint8_t i = 0, j = 0, k = 0, m = 0, t = num - 1;
    float angle;
    complex w, temp;

    //变址，把自然顺序变为倒位序
    for (i = 0; i < t; i++)
    {
        if (i < j)
        {
            temp = data[j];
            data[j] = data[i];
            data[i] = temp;
        }
        k = N2;
        while (k <= j)
        {
            j -= k;
            k >>= 1;
        }

        j += k;    
    }

    for (i = 1; i <= n; i++)
    {
        uint8_t km, step = 1 << i;
        m = step >> 1;

        for (j = 0; j < m; j++)
        {
            angle = ((double)j) / m;

            if (dir == 1)
                w.imag = -find_sin(angle);
            else if (dir == -1)
                w.imag = find_sin(angle);
            w.real = find_cos(angle);

            for (k = j; k < num; k += step)
            {
                km = k + m;
                //用下面两行直接计算复数乘法，省去函数调用开销
                temp.real = data[km].real * w.real - data[km].imag * w.imag;
                temp.imag = w.imag * data[km].real + w.real * data[km].imag;

                data[km].real = data[k].real - temp.real;
                data[km].imag = data[k].imag - temp.imag;

                data[k].real = data[k].real + temp.real;
                data[k].imag = data[k].imag + temp.imag;            
            }
        }    
    }

    if (dir == -1)
    {
        for (i = 0; i < num; i++)
        {
            data[i].real /= num;
        }
    }
}


int main()
{
    int i;
    complex data[N];

    //初始化输入数据
    for (i = 0; i < N; i++)
    {
        data[i].real = pr[i];
        data[i].imag = 0;
    }

    //创建正弦表
    for (i = 0; i <= N4; i++)
    {
        sin_table[i] = (float)sin(PI * i /N2);
    }

    //正变换
    fft(data, N, M, 1);

    for (i = 0; i < N; i++)
    {
        printf("%.f    %.f
", data[i].real, data[i].imag);
    }
    printf("
");

    //波形处理
    for (i = 5; i < N; i++)    //10:3.9Hz,11:4.3Hz
    {
        data[i].imag= 0;
        data[i].real = 0;
    }

    //逆变换
    fft(data, N, M, -1);

    for (i = 0; i < N; i++)
    {
        printf("%.f    %.f
", data[i].real, data[i].imag);
    }
    printf("
");

    return 0;
}