/*價4 fft*/
#include <stdio.h>
#include <math.h>
#if defined(__GNUC__)
typedef long long int64_t;
typedef unsigned long long uint64_t;
#elif defined(WIN32)
typedef __int64 int64_t;
typedef unsigned __int64 uint64_t;
#endif
#define pi 3.14159265358979
#define TEST
typedef struct
{
float re;
float im;
}comp;
int BitReverse(int src, int size)
{
int tmp = src;
int des = 0;
int i;
for (i=size-1; i>=0; i--)
{
des = ((tmp & 0x3) << (i*2)) | des;
tmp = tmp >> 2;
}
return des;
}
void reverse_idx(comp *in,int log4_N)
{
int i;
int N = 1 << (log4_N*2);
comp * temp;
temp = malloc(N * sizeof(comp));
if(!temp)
{
printf("malloc failed!\n");
exit(0);
}
for(i = 0; i < N;i++)
{
int idx;
idx = BitReverse(i,log4_N);
temp[idx].re = in[i].re;
temp[idx].im = in[i].im;
}
for(i = 0; i < N; i++)
{
in[i].re = temp[i].re;
in[i].im = temp[i].im;
}
free(temp);
}
void fft_ifft_4_common(comp *in,comp * win,int log4_N,int reverse)
{
int N = (1 << log4_N * 2);
int i,j,k;
int span = 1;
int n = N >> 2;
float wr,wi;
int widx;
comp temp1,temp2,temp3,temp4;
int idx1,idx2,idx3,idx4;
int **cof;
for(i = 0;i < log4_N; i++)
{
for(j = 0; j < n; j++)
{
widx = 0;
idx1 = j * span * 4;
idx2 = idx1 + span;
idx3 = idx2 + span;
idx4 = idx3 + span;
for(k = 0; k < span; k++)
{
temp1.re = in[idx1 + k].re;
temp1.im = in[idx1 + k].im;
temp2.re = win[widx].re * in[idx2 + k].re - win[widx].im * in[idx2 + k].im;
temp2.im = win[widx].im * in[idx2 + k].re + win[widx].re * in[idx2 + k].im;
temp3.re = win[widx * 2].re * in[idx3 + k].re -win[widx * 2].im * in[idx3 + k].im;
temp3.im = win[widx * 2].im * in[idx3 + k].re + win[widx * 2].re * in[idx3 + k].im;
temp4.re = win[widx * 3].re * in[idx4 + k].re -win[widx * 3].im * in[idx4 + k].im;
temp4.im = win[widx * 3].im * in[idx4 + k].re + win[widx * 3].re * in[idx4 + k].im;
in[idx1 + k].re = temp1.re + temp3.re;
in[idx1 + k].im = temp1.im + temp3.im;
in[idx2 + k].re = temp1.re - temp3.re;
in[idx2 + k].im = temp1.im - temp3.im;
in[idx3 + k].re = temp2.re + temp4.re;
in[idx3 + k].im = temp2.im + temp4.im;
in[idx4 + k].re = temp2.re - temp4.re;
in[idx4 + k].im = temp2.im - temp4.im;
temp1.re = in[idx1 + k].re + in[idx3 + k].re;
temp1.im = in[idx1 + k].im + in[idx3 + k].im;
if(reverse == 0)
{
temp2.re = in[idx2 + k].re + in[idx4 + k].im;
temp2.im = in[idx2 + k].im - in[idx4 + k].re;
}
else
{
temp2.re = in[idx2 + k].re - in[idx4 + k].im;
temp2.im = in[idx2 + k].im + in[idx4 + k].re;
}
temp3.re = in[idx1 +k].re - in[idx3 + k].re;
temp3.im = in[idx1 + k].im - in[idx3 + k].im;
if(reverse == 0)
{
temp4.re = in[idx2 + k].re - in[idx4 + k].im;
temp4.im = in[idx2 + k].im + in[idx4 + k].re;
}
else
{
temp4.re = in[idx2 + k].re + in[idx4 + k].im;
temp4.im = in[idx2 + k].im - in[idx4 + k].re;
}
in[idx1 + k].re = temp1.re;
in[idx1 + k].im = temp1.im;
in[idx2 + k].re = temp2.re;
in[idx2 + k].im = temp2.im;
in[idx3 + k].re = temp3.re;
in[idx3 + k].im = temp3.im;
in[idx4 + k].re = temp4.re;
in[idx4 + k].im = temp4.im;
widx += n;
}
}
n >>= 2;
span <<= 2;
}
}
void fft4(comp *in,int log4_N)
{
comp *win;
int N = 1 << (log4_N*2);
int i;
win = malloc((3*N/4 - 2) * sizeof(comp));
if(!win)
{
printf("malloc failed!\n");
exit(0);
}
reverse_idx(in,log4_N);
for(i = 0; i < (3*N/4-2); i++)
{
win[i].re = cos(2*pi*i/(float)N);
win[i].im = -sin(2*pi*i/(float)N);
}
fft_ifft_4_common(in,win,log4_N,0);
free(win);
}
void ifft4(comp *in,int log4_N)
{
int N = 1 << (log4_N * 2);
comp *win;
int i;
win = malloc((3*N/4 - 2) * sizeof(comp));
if(!win)
{
printf("malloc failed!\n");
exit(0);
}
reverse_idx(in,log4_N);
for(i = 0; i < (3*N/4-2); i++)
{
win[i].re = cos(2*pi*i/(float)N);
win[i].im = sin(2*pi*i/(float)N);
}
fft_ifft_4_common(in,win,log4_N,1);
for(i = 0; i < N; i++)
{
in[i].re /= (float)N;
in[i].im /= (float)N;
}
free(win);
}
#ifdef TEST
int main()
{
int i;
int N = 64;
int log4_N = 3;
comp in[64];
for(i = 0;i < 32; i++)
{
in[i].re = 1;
in[i].im = 0;
}
for(i = 32;i < 64;i++)
{
in[i].re = 2;
in[i].im = 1;
}
for(i = 0; i < N; i++)
{
printf("in[%d]: %f + %fi\n",i,in[i].re,in[i].im);
}
fft4(in,log4_N);
printf("fft:\n");
for(i = 0; i < N; i++)
{
printf("out[%d]: %f + (%f)i\n",i,in[i].re,in[i].im);
}
ifft4(in,log4_N);
printf("ifft:\n");
for(i = 0; i < N; i++)
{
printf("out[%d]: %f + (%f)i\n",i,in[i].re,in[i].im);
}
return 0;
}
#endif
评论2