/*--------------------------------------------------------------------------------
高速フーリエ変換と逆高速フーリエ変換
Fast Fourier transform and inverse Fast Fourier transform
calc_FFT.c
*--------------------------------------------------------------------------------*/
#include "FFT.h"
// 高速フーリエ変換
// Fast Fourier transform
void FFT( int num, complex *x, complex *F )
{
int i, j, k, r, p, j1, j2, k1, bit;
double a;
double b;
double angle;
for( i = 0; i < num; i++ ) {
F[i] = x[i];
}
k = 0;
j2 = num;
r = log2( (double)num );
j1 = r - 1;
for( j = 0; j < r; j++ ) {
j2 /= (double)2.0;
while(1) {
for( i = 0; i < j2; i++ ) {
p = k/(double)pow( 2.0, (double)j1 );
k1 = k + j2;
angle = 2.0 * M_PI * bitr(p, r)/(double)num;
a = F[k1].real * cos( angle ) + F[k1].imag * sin( angle );
b = - F[k1].real * sin( angle ) + F[k1].imag * cos( angle );
/*--------------------------------------------------------*/
// バタフライ演算処理
// butterfly operation
F[k1].real = F[k].real - a;
F[k1].imag = F[k].imag - b;
F[k].real = F[k].real + a;
F[k].imag = F[k].imag + b;
/*--------------------------------------------------------*/
k++;
}
k += j2;
if( k >= num ) {
break;
}
}
k = 0;
j1--;
}
// ビット反転処理
// bit inversion processing
for( k = 0; k < num; k++ ) {
bit = bitr( k, r );
if( bit > k ) {
swap( &F[k], &F[bit]);
}
}
}
/*-------------------------------------------------------------------------------*/
// 逆高速フーリエ変換
// inverse Fast Fourier transform
void inv_FFT( int num, complex *x, complex *F )
{
int i, j, k, r, p, j1, j2, k1, bit;
double a;
double b;
double angle;
for( i = 0; i < num; i++ ) {
x[i] = F[i];
}
k = 0;
j2 = num;
r = log2( (double)num );
j1 = r - 1;
for( j = 0; j < r; j++) {
j2 /= (double)2.0;
while(1) {
for( i = 0; i < j2; i++ ) {
p = k/(double)pow( 2.0, (double)j1 );
k1 = k + j2;
angle = 2.0 * M_PI * bitr(p, r)/(double)num;
a = x[k1].real * cos( angle ) - x[k1].imag * sin( angle );
b = x[k1].real * sin( angle ) + x[k1].imag * cos( angle );
/*--------------------------------------------------------*/
// バタフライ演算処理
// butterfly operation
x[k1].real = x[k].real - a;
x[k1].imag = x[k].imag - b;
x[k].real = x[k].real + a;
x[k].imag = x[k].imag + b;
/*--------------------------------------------------------*/
k++;
}
k += j2;
if( k >= num ) {
break;
}
}
k = 0;
j1--;
}
// ビット反転処理
// bit inversion processing
for( k = 0; k < num; k++ ) {
bit = bitr( k, r );
if( bit > k ) {
swap( &x[k], &x[bit] );
}
}
for( k = 0; k < num; k++ ) {
x[k].real = x[k].real/(double)num;
x[k].imag = x[k].imag/(double)num;
}
}
/*-------------------------------------------------------------------------------*/
int bitr( int bit, int r )
{
int i, bitr;
bitr = 0;
for( i = 0; i < r; i++ ) {
bitr *= 2.0;
bitr |= ( bit & 1 );
bit /= (double)2.0;
}
return ( bitr );
}
/*-------------------------------------------------------------------------------*/
void swap( complex *a, complex *b )
{
complex c;
c = (*b);
(*b) = (*a);
(*a) = c;
}
/*-------------------------------------------------------------------------------*/