/* @(#)fft_pow2.c 17.1.1.1 (ES0-DMD) 01/25/02 17:21:21 */ /*=========================================================================== Copyright (C) 1995 European Southern Observatory (ESO) This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version. This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software Foundation, Inc., 675 Massachusetss Ave, Cambridge, MA 02139, USA. Corresponding concerning ESO-MIDAS should be addressed as follows: Internet e-mail: midas@eso.org Postal address: European Southern Observatory Data Management Division Karl-Schwarzschild-Strasse 2 D 85748 Garching bei Muenchen GERMANY ===========================================================================*/ /************************************************************************** ** Copyright (C) 1993 by European Southern Observatory *************************************************************************** ** ** UNIT ** ** Version: 17.1.1.1 ** ** Author: Jean-Luc Starck ** ** Date: 02/01/25 ** ** File: fft_pow2.c ** **************************************************************************** ** ** DESCRIPTION ** ----------- ** ** this module contains the routines converning the Fourier transform ** **************************************************************************** ** ** ft_cf_any_power_of_2(dat,direction,length) ** complex_float *dat; ** int direction; ** int length; ** ** Takes the array pointed to by dat which is assumed to be of side-length ** length and performs a fourier transform of its contents ** ** External function calls ** None ** ** Return codes ** 0 : No problems encountered. ** 1 : Length was not a power of 2 FATAL_ERROR ** ** INPUT dat = complex_floats images ** direction = 1 ==> directe Fourier transform ** direction = -1 ==> inverse Fourier transform ** length = number of lines X number of columns ** ************************************************************************** ** ** prepare_fft_real (Pict, FFT_Dat, N) ** complex_float *FFT_Dat; ** float *Pict; ** int N; ** ** copy a real image to a complex image with a imaginary part equal to ** zero. ** **************************************************************************/ static char sccsid[] = "@(#)fft_pow2.c 17.1.1.1 02/01/25 ESO @(#)"; #include #include #include "Def_Math.h" #include "Def_Mem.h" /* extern char *malloc(); */ typedef struct { float cos; float sin; } trig_table_struct; typedef struct { int table_len; int length; int len_exp; trig_table_struct *table; } trig_table_type; typedef struct { int source; int destin; } rev_map_struct; typedef struct { int map_len; int length; rev_map_struct *map; } rev_map_type; /*************************2222222222222222222****************************/ static complex_float *b_s_row_to_vector(dat,index,length) complex_float *dat; int index; int length; { complex_float transfer; register int loop; complex_float *vector,*vrem; vrem = vector = (dat + index*length); length /= 2; for ( loop = 0 ; loop < length ; ++loop, ++vector){ transfer = *(vector); *(vector) = *(vector+ length); *(vector+length) = transfer; } return (vrem); } /******************************22222222222*********************************/ static int b_s_row_from_vector(vector,length) complex_float *vector; int length; { /* this routine assumes that the correct address of thr row vector */ /* within the data array is correctly specified by the pointer vector*/ /* ie vector has not been changed since it was set by b_s_row_to_vector*/ complex_float transfer; register int loop; length /= 2; for ( loop = 0 ; loop < length ; ++loop, ++vector){ transfer = *(vector); *(vector) = *(vector+ length); *(vector+length) = transfer; } } /******************************22222222222*********************************/ static int b_s_col_to_vector(dat,vector,index,length) complex_float *dat; complex_float *vector; int index; int length; { register int ofst; register int loop; register int len_b_2; len_b_2 = length/2; ofst = len_b_2*length; dat += index; for ( loop = 0 ; loop < len_b_2 ; ++loop, ++vector, dat += length){ *(vector) = *(dat+ ofst); *(vector+len_b_2) = *dat; } } /******************************22222222222*********************************/ static int b_s_col_from_vector(dat,vector,index,length) complex_float *dat; complex_float *vector; int index; int length; { register int ofst; register int loop; register int len_b_2; len_b_2 = length/2; ofst = len_b_2*length; dat += index; for ( loop = 0 ; loop < len_b_2 ; ++loop , ++vector , dat += length){ *(dat + ofst) = *vector; *dat = *(vector+len_b_2); } } /******************************22222222222*********************************/ static int bitreverse(vector,rev_map) rev_map_type rev_map; complex_float *vector; { register rev_map_struct *map; complex_float transfer; register int loop; map = rev_map.map; for ( loop = 0 ; loop < rev_map.map_len ; ++loop, ++map){ transfer = *(vector+map->destin); *(vector+map->destin) = *(vector+map->source); *(vector+map->source) = transfer; } } /******************************22222222222*********************************/ static int ft_1d(vector,trig_table) trig_table_type trig_table; complex_float *vector; { register float c,s; register int butterfly = 1 ; int rank; int l2,end_of_loop2; register complex_float *top_btrfly_ptr,*mid_btrfly_ptr; register int trig_table_pos; register float re_transfer_buffer,im_transfer_buffer; trig_table_struct *table; table = trig_table.table; for (rank = 1 ; rank <= trig_table.len_exp ; ++rank){ top_btrfly_ptr = vector +0; mid_btrfly_ptr =vector+ butterfly; end_of_loop2 = trig_table.length >> rank;/* 2^(len_exp - rank) */ for ( l2 = 1 ; l2 <= end_of_loop2 ; ++l2){ for ( trig_table_pos = 0 ; trig_table_pos < trig_table.table_len ; trig_table_pos += end_of_loop2 ){ c = (table+trig_table_pos)->cos; s = (table+trig_table_pos)->sin; re_transfer_buffer = (*mid_btrfly_ptr).re*c - (*mid_btrfly_ptr).im*s; im_transfer_buffer = (*mid_btrfly_ptr).im*c + (*mid_btrfly_ptr).re*s; (*mid_btrfly_ptr).re = (*top_btrfly_ptr).re - re_transfer_buffer; (*mid_btrfly_ptr).im = (*top_btrfly_ptr).im - im_transfer_buffer; (*top_btrfly_ptr).re += re_transfer_buffer; (*top_btrfly_ptr).im += im_transfer_buffer; top_btrfly_ptr += 1; mid_btrfly_ptr += 1; } top_btrfly_ptr += butterfly; mid_btrfly_ptr += butterfly; } butterfly *= 2; } } /******************************22222222222*********************************/ static int set_trig_table(trig_table,direction) trig_table_type trig_table; int direction; { register int l; double arg; trig_table_struct *table; table = trig_table.table; for ( l = 0 ; (l < trig_table.table_len) ; ++l, ++table){ arg =direction*PI*l/trig_table.table_len; table->cos = cos(arg); table->sin = sin(arg); } } /******************************22222222222*********************************/ static int set_rev_map(rev_map) rev_map_type rev_map; { register int source,destin,s_mask,d_mask ; rev_map_struct *map; int adrs = 0; map = rev_map.map; for ( source = 0 ; source < rev_map.length -1 ; ++source ) { destin = 0; d_mask = (int) rev_map.length/2; for ( s_mask = 1 ; (s_mask < rev_map.length) ; ) { if ((s_mask & source) > 0) { destin = destin | d_mask; } d_mask = d_mask >> 1; s_mask = s_mask << 1; } if ( destin > source) { (map+adrs)->source = source; (map+adrs)->destin = destin; ++adrs; } } } /******************************22222222222*********************************/ static int print_array(ary,length) complex_float *ary; int length; { int l1,l2; for( l1 = 0 ; l1 < length ; ++l1){ for( l2 = 0 ; l2 < length ; ++l2, ++ary){ printf("%6f ",sqrt(ary->re*ary->re+ary->im*ary->im)); } printf("\n"); } } /******************************22222222222*********************************/ static normalisation (ary,length) complex_float *ary; int length; { int l1,l2; for( l1 = 0 ; l1 < length ; ++l1) { for( l2 = 0 ; l2 < length ; ++l2, ++ary) { ary->re /= (float)(length * length); ary->im /= (float)(length * length); } } } /*************************11111111111111118*******************************/ ft_cf_any_power_of_2(dat,direction,length) complex_float *dat; int direction; int length; /* Takes the array pointed to by dat which is assumed to be of side-length length and performs a fourier transform of its contents External function calls None Return codes 0 : No problems encountered. 1 : Length was not a power of 2 FATAL_ERROR */ { int error = 0; complex_float *vector; complex_float *vect_space_holder; int index; rev_map_type rev_map; trig_table_type trig_table; int len_exp; int temp; /* First we must check length to make sure it is a power of 2 */ len_exp = (int)(0.3+log((double)(length))/(log(2.0))); INT_POW(2,len_exp,temp); if ( length != temp ){ io_err_message_exit (ERR_POWER_OF_2, " "); } else{ /* Now we must reserve space for the various working arrays */ INT_POW(2,len_exp-1,rev_map.map_len); INT_POW(2,(len_exp-1)/2,temp); rev_map.map_len -= temp; rev_map.length = length; rev_map.map = (rev_map_struct*) calloc((unsigned)sizeof(rev_map_struct)*rev_map.map_len,1); INT_POW(2,len_exp-1,trig_table.table_len); trig_table.length = length; trig_table.len_exp = len_exp; trig_table.table = (trig_table_struct*) calloc((unsigned)sizeof(trig_table_struct)*trig_table.table_len,1); vect_space_holder = (complex_float*) calloc((unsigned)sizeof(complex_float)*length,1); set_trig_table(trig_table,direction); set_rev_map(rev_map); /* row vectors actually occur in the data array so all we need to do*/ /* is referance their start position */ for (index = 0 ; index < length ; ++index){ vector = b_s_row_to_vector(dat,index,length); bitreverse(vector,rev_map); ft_1d(vector,trig_table); b_s_row_from_vector(vector,length); } /* column vectors are dispersed in the data array and thus have to */ /* collected together the variable vect_space_holder reserves enough*/ /* space in thi routine to hold the colledted vector */ vector = &(vect_space_holder[0]); for (index = 0 ; index < length ; ++index){ b_s_col_to_vector(dat,vector,index,length); bitreverse(vector,rev_map); ft_1d(vector,trig_table); b_s_col_from_vector(dat,vector,index,length); } /*Now free-up space used by the working arrays */ free((char*)vect_space_holder); free((char*)trig_table.table); free((char*)rev_map.map); } if (direction == -1) normalisation (dat,length); return(error); } /******************************22222222222*********************************/ prepare_fft_real (Pict, FFT_Dat, N) complex_float *FFT_Dat; float *Pict; int N; { int ind = 0; int l1,l2; for ( l1 = 0 ; l1 < N ; l1++) { for ( l2 = 0 ; l2 < N ; l2++) { FFT_Dat [ind].re = Pict [ind]; FFT_Dat [ind].im = 0.; ind ++; } } } /************************************************************************/