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