Doxygen Source Code Documentation
csfft_AJ.c File Reference
#include "mrilib.h"
Go to the source code of this file.
Defines | |
#define | PI (3.141592653589793238462643) |
Functions | |
void | csfft_trigconsts (int idim) |
void | csfft (int mode, int idim, complex *xc) |
Variables | |
complex * | csplus = NULL |
complex * | csminus = NULL |
int | nold = -666 |
Define Documentation
|
Definition at line 12 of file csfft_AJ.c. Referenced by csfft_trigconsts(). |
Function Documentation
|
Definition at line 86 of file csfft_AJ.c. References csfft_trigconsts(), complex::i, i1, i2, nold, complex::r, and xc.
00087 { 00088 register unsigned int m, n, i0, i1, i2, i3, k; 00089 register complex *r0, *r1, *csp; 00090 register float co, si, f0, f1, f2, f3, f4; 00091 double al; 00092 00093 /** perhaps initialize **/ 00094 00095 if( nold != idim ) csfft_trigconsts( idim ) ; 00096 00097 /** Main loop starts here **/ 00098 00099 n = idim; 00100 i2 = idim >> 1; 00101 i1 = 0; 00102 csp = (mode > 0) ? csplus : csminus ; /* choose const array */ 00103 00104 for (i0=0; i0 < n; i0 ++) { 00105 if ( i1 > i0 ) { 00106 r0 = xc + i0; 00107 r1 = xc + i1; 00108 f1 = r0->r; 00109 f2 = r0->i; 00110 r0->r = r1->r; 00111 r0->i = r1->i; 00112 r1->r = f1; 00113 r1->i = f2; 00114 } 00115 m = i2; 00116 while ( m && !(i1 < m) ) { 00117 i1 -= m; 00118 m >>= 1; 00119 } 00120 i1 += m; 00121 } 00122 00123 m = 1; 00124 k = 0; 00125 while (n > m) { 00126 i3 = m << 1; 00127 for (i0=0; i0 < m; i0 ++) { 00128 co = (csp + k)->r; 00129 si = (csp + k)->i; 00130 for (i1=i0; i1 < n; i1 += i3) { 00131 r0 = xc + i1; 00132 r1 = r0 + m; 00133 f1 = r1->r * co; 00134 f2 = r1->i * si; 00135 f3 = r1->r * si; 00136 f4 = r1->i * co; 00137 f1 -= f2; 00138 f3 += f4; 00139 f2 = r0->r; 00140 f4 = r0->i; 00141 r1->r = f2 - f1; 00142 r1->i = f4 - f3; 00143 r0->r = f2 + f1; 00144 r0->i = f4 + f3; 00145 } 00146 k++; 00147 } 00148 m = i3; 00149 } 00150 00151 #ifdef SCALE_INVERSE 00152 if (mode > 0) { 00153 f0 = 1.0 / idim ; 00154 i0 = 0; 00155 i1 = 1; 00156 while (i0 < n) { 00157 r0 = xc + i0; 00158 r1 = xc + i1; 00159 f1 = r0->r; 00160 f2 = r0->i; 00161 f3 = r1->r; 00162 f4 = r1->i; 00163 f1 *= f0; 00164 f2 *= f0; 00165 f3 *= f0; 00166 f4 *= f0; 00167 r0->r = f1; 00168 r0->i = f2; 00169 r1->r = f3; 00170 r1->i = f4; 00171 i0 += 2; 00172 i1 += 2; 00173 } 00174 } 00175 #endif 00176 00177 return ; 00178 } |
|
Definition at line 14 of file csfft_AJ.c. References free, complex::i, i1, i2, malloc, nold, PI, and complex::r.
00015 { 00016 register unsigned int m, n, i0, i1, i2, i3, k; 00017 register complex *r0, *csp; 00018 register float co, si, f0, f1, f2, f3, f4; 00019 double al; 00020 00021 if( idim == nold ) return ; 00022 00023 if( idim > nold ){ 00024 if( csplus != NULL ){ free(csplus) ; free(csminus) ; } 00025 csplus = (complex *) malloc( sizeof(complex) * idim ) ; 00026 csminus = (complex *) malloc( sizeof(complex) * idim ) ; 00027 if( csplus == NULL || csminus == NULL ){ 00028 fprintf(stderr,"\n*** fft cannot malloc space! ***\n"); exit(-1) ; 00029 } 00030 } 00031 nold = n = idim ; 00032 00033 f1 = 1.0 ; /* csplus init */ 00034 m = 1; 00035 k = 0; 00036 while (n > m) { 00037 i3 = m << 1; 00038 f2 = m; 00039 al = f1*PI/f2; 00040 co = cos(al); 00041 si = sin(al); 00042 (csplus + k)->r = 1.; 00043 (csplus + k)->i = 0.; 00044 for (i0=0; i0 < m; i0++) { 00045 k++; 00046 csp = csplus + k; 00047 r0 = csp - 1; 00048 csp->r = r0->r * co - r0->i * si; 00049 csp->i = r0->i * co + r0->r * si; 00050 } 00051 m = i3; 00052 } 00053 00054 f1 = -1.0 ; /* csminus init */ 00055 m = 1; 00056 k = 0; 00057 while (n > m) { 00058 i3 = m << 1; 00059 f2 = m; 00060 al = f1*PI/f2; 00061 co = cos(al); 00062 si = sin(al); 00063 (csminus + k)->r = 1.; 00064 (csminus + k)->i = 0.; 00065 for (i0=0; i0 < m; i0++) { 00066 k++; 00067 csp = csminus + k; 00068 r0 = csp - 1; 00069 csp->r = r0->r * co - r0->i * si; 00070 csp->i = r0->i * co + r0->r * si; 00071 } 00072 m = i3; 00073 } 00074 return ; 00075 } |
Variable Documentation
|
Definition at line 8 of file csfft_AJ.c. |
|
Definition at line 8 of file csfft_AJ.c. |
|
Definition at line 9 of file csfft_AJ.c. Referenced by csfft(), and csfft_trigconsts(). |