FFT

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

//**************************************************************************************

#ifndef M_PI
    #define M_PI 3.14159265358979323846
#endif
#define DATA_LEN 64
#define SAMPLE_FREQ 1000 // Hz

unsigned char sin_data[64] = {0x7F,0x8B,0x98,0xA4,0xB0,0xBB,0xC6,0xD0,0xD9,0xE2,
    0xE9,0xEF,0xF5,0xF9,0xFC,0xFE,0xFE,0xFE,0xFC,0xF9,0xF5,0xEF,0xE9,0xE2,0xD9,0xD0,
    0xC6,0xBB,0xB0,0xA4,0x98,0x8B,0x7F,0x73,0x66,0x5A,0x4E,0x43,0x38,0x2E,0x25,0x1C,
    0x15,0x0F,0x09,0x05,0x02,0x00,0x00,0x00,0x02,0x05,0x09,0x0F,0x15,0x1C,0x25,0x2E,
    0x38,0x43,0x4E,0x5A,0x66,0x73};

typedef struct{
    double r;
    double i;
} cplx_t;


//**************************************************************************************

void cplx_exp(cplx_t *x, cplx_t *r)
{
    double expx = exp(x->r);
    r->r = expx*cos(x->i);
    r->i = expx*sin(x->i);
}

// 复数乘法
void cplx_mul(cplx_t *x, cplx_t *y, cplx_t *r)
{
    r->r = x->r*y->r-x->i*y->i;
    r->i = x->r*y->i+x->i*y->r;
}

// 比特反置
void bit_reverse(cplx_t *x, int N)
{
    unsigned int i = 0,j = 0,k = 0;
    
    cplx_t tmp;
    
    int bit_num = log(0.0+N)/log(2.0);  // 比特位数
    
    for (i=0; i<N; i++)
    {
        int t = bit_num;
        k=i;
        j=0;
        while (t--)
        {
            j <<= 1;
            j |= k&1;
            k >>= 1;
        }
        if (j>i)
        {
            tmp = x[i];
            x[i] = x[j];
            x[j] = tmp;
        }
    }
}


void fft(cplx_t *x, int N)
{
    cplx_t u,d,p,W,tmp;
    int i=0,j=0,k=0,l=0;
    
    double M = floor(log(0.0+N)/log(2.0));  // zhiqiu 换底公式
    if (log(0.0+N)/log(2.0)-M > 0)
    {
        printf("The length of x (N) must be a power of two!!!\n");
        return;
    }
    
    bit_reverse(x,N);
    
    for (i = 0; i < M; i++)
    {
        l = 1<<i;
        for (j = 0; j < N; j += 2*l)
        {
            for (k = 0; k < l; k++)
            {
                tmp.r = 0.0;
                tmp.i = -2*M_PI*k/2/l;
                cplx_exp(&tmp,&W);
                cplx_mul(&x[j+k+l],&W,&p);
                u.r = x[j+k].r+p.r;
                u.i = x[j+k].i+p.i;
                d.r = x[j+k].r-p.r;
                d.i = x[j+k].i-p.i;
                x[j+k] = u;
                x[j+k+l] = d;
            }
        }
    }
}

void ifft(cplx_t *x, int N)
{
    unsigned int i = 0;
    for (i = 0;i < N; i++)
        x[i].i = -x[i].i;
    fft(x,N);
    for (i = 0;i < N; i++)
    {
        x[i].r = x[i].r/(N+0.0);
        x[i].i = -x[i].i/(N+0.0);
    }
}



//**************************************************************************************

int main(int argc, const char * argv[])
{
    int i;
    cplx_t x[DATA_LEN];
    for (i=0;i<DATA_LEN;i++)
    {
        x[i].r = sin_data[i];
        x[i].i = 0;
    }
    
    printf("Before...\nReal\t\tImag\n");
    for (i=0; i<DATA_LEN; i++)
        printf("%f\t%f\n",x[i].r,x[i].i);
    
    fft(x, DATA_LEN);
    
    printf("\nAfter FFT...\nReal\t\tImag\n");
    for (i=0; i<DATA_LEN; i++)
        printf("%f\t%f\n",x[i].r,x[i].i);
    
    printf("\n频率\t\t幅度\t\t相位\n");
    for (i=0; i<DATA_LEN; i++)
    {
        double fudu = sqrt(x[i].r*x[i].r+x[i].i*x[i].i);
        double xiangwei = atan(x[i].i/x[i].r)*180.0/M_PI;
        if (i == 0)
            fudu /= DATA_LEN;
        else
            fudu /= DATA_LEN/2;
        printf("%f\t%f\t%f\n", (double)SAMPLE_FREQ/DATA_LEN*i, fudu, xiangwei);
    }
    
    ifft(x, DATA_LEN);
    
    printf("\nAfter IFFT...\nReal\t\tImag\n");
    for (i=0; i<DATA_LEN; i++)
        printf("%f\t%f\n",x[i].r,x[i].i);
    
    return 0;
}