#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 // Hzunsigned 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;}