fourier.cpp

复数类

#include "fourier.h"
#include <stdlib.h>
#include <stdio.h>
#include <math.h>

double fourier::PI = 3.14159265358979323846264338327950288 ;
	
fourier::fourier(){
	w = NULL ;
}

fourier::~fourier(){
	if ( w != NULL ) free(w);
}

//exp(-j * ( 2 * pi * i / N ) ) 
//size W
//dir: 1 , increase
//dir: -1 , decrease
void fourier::initw(int size , int dir){
	if ( w != NULL ) free(w);
	w = (complex *) malloc(size*sizeof(complex)) ;
	int i ;
	if ( dir == 1){
		for ( i = 0 ; i < size ; i ++ ){
			w[i].real = cos( 2 * PI * i / size ) ;
			w[i].imag = -sin( 2 * PI * i / size ) ;
		}
	}
	else{
		for ( i = 0 ; i < size ; i ++ ){
			w[i].real = cos( 2 * PI * (-i) / size ) ;
			w[i].imag = -sin( 2 * PI * (-i) / size ) ;
		}	
	}
}

int fourier::rbit(int b , int m){
	int t = 0 ;
	int i = 1 ;
	while ( m -- ) {
		t <<= 1 ;
		t |= b & 1 ;
		b >>= 1;
	}
	return t ;
}

//bit reverse order for fft caculating
void fourier::bitreverse(complex *x , int size , int m){
	int i , j ;
	complex tmp ;
	for ( i = 0 ; i < size ; i ++ ){
		j = rbit(i , m);
		if ( j > i ) {
			tmp = x[i] ;
			x[i] = x[j] ;
			x[j] = tmp ;
		}
	}
}

void fourier::output(complex *x , int size){
	int i;
	for ( i = 0 ; i < size ; i ++ ){
		x[i].output() ;
		putchar('\n');
	}
}

bool fourier::dft(complex *x , int size){
	initw(size , 1);
	complex *tmp = ( complex * ) malloc( size * sizeof ( complex )) ;
	int i , j ;
	for ( i = 0 ; i < size ; i ++ ) {
		tmp[i].real = 0 ;
		tmp[i].imag = 0 ;
		for ( j = 0 ; j < size ; j ++ ){
			tmp[i] += x[j] * w[ i * j % size ] ;
		}
	}
	for ( i = 0 ; i < size ; i ++ ) {
		x[i] = tmp[i] ;
	}
	free(tmp) ;
	return true ;
}

bool fourier::idft(complex *x , int size){
	initw(size , -1) ;
	complex *tmp = ( complex * ) malloc( size * sizeof ( complex )) ;
	int i , j ;
	for ( i = 0 ; i < size ; i ++ ) {
		tmp[i].real = 0 ;
		tmp[i].imag = 0 ;
		for ( j = 0 ; j < size ; j ++ ){
			tmp[i] += x[j] * w[ i * j % size ] ;
		}
	}
	for ( i = 0 ; i < size ; i ++ ) {
		x[i].real = tmp[i].real / size ;
		x[i].imag = tmp[i].imag / size ;
	}
	free(tmp) ;
	return true ;
}

void fourier::fft_ifft_common(complex *x , int size){
	int i , j ;
	int m ;
	complex high , low ;

	//caculate power m 
	i = size ;
	m = 0 ;
	while ( i ) {
		m ++ ;
		i >>= 1 ;
	}

	bitreverse(x , size , m - 1 ) ;

	int group_num = size / 2 ;
	int group_length = 2 ;
	int group_step = 1 ;
	while ( m -- ){
		for ( i = 0 ; i < group_num ; i ++ ){
			int group_start = i * group_length ;
			for ( j = 0 ; j < group_step ; j ++ ){
				int a = group_start + j ;
				int b = group_start + j + group_step ;
				high = x[a] ;
				low = x[b] * w[j * group_num];
				x[a] = high + low ;
				x[b] = high - low ;
			}
		}
		group_num >>= 1 ;
		group_length <<= 1 ;
		group_step <<= 1 ;
	}
}

//pre-condition: size = 2 ^ m
//post-condition: x be the FT of x
bool fourier::fft(complex *x , int size){

	initw(size , 1) ;
	fft_ifft_common(x , size);
	
	return true ;
}

bool fourier::ifft(complex *x , int size){
	int i;

	initw(size , -1) ;
	fft_ifft_common(x , size);

	for ( i = 0 ; i < size ; i ++ ) {
		x[i].imag /= size ;
		x[i].real /= size ;
	}

	return true ;
}