Doxygen Source Code Documentation
mri_cfft.c File Reference
#include <math.h>#include <stdlib.h>#include <stdio.h>Go to the source code of this file.
Defines | |
| #define | PI 3.1415926536 |
| #define | IDMAX 1024 |
| #define | LIDMAX 10 |
Functions | |
| void | cfft (int mode, int idim, float *xr, float *xi) |
| void | cfft2d (int mode, int nx, int ny, float *xr, float *xi) |
Define Documentation
|
|
|
|
|
|
|
|
Definition at line 12 of file mri_cfft.c. Referenced by cfft(). |
Function Documentation
|
||||||||||||||||||||
|
Definition at line 23 of file mri_cfft.c. Referenced by cfft2d().
00025 {
00026 #define IDMAX 1024
00027 #define LIDMAX 10
00028 static int idold = -999 ;
00029 register int i0,i1,i2,i4,i5,m0, id = idim;
00030 float f1,f3,f4,f5,al,co,si,md = mode;
00031 static int n,m[LIDMAX];
00032 static float f2,c[IDMAX/2],s[IDMAX/2];
00033
00034 /**** preparations if id has changed since last call ****/
00035
00036 if (idold != id) {
00037 idold = id ;
00038 /* check for power of 2 */
00039 for( i4=4 ; i4 <= IDMAX ; i4 *= 2 ){
00040 if( id == i4 )break ;
00041 }
00042 if( id != i4 ){
00043 fprintf(stderr,"\n In cfft : illegal idim=%d\n",idim);
00044 exit(1) ;
00045 }
00046 f2 = id;
00047 n = log(f2)/log(2.) + .5;
00048 m[n-1] = 1;
00049 al = 2.*PI/f2;
00050 co = cos(al);
00051 si = sin(al);
00052 c[0] = 1.;
00053 s[0] = 0.;
00054 for(i4 = 1; i4 < 512; i4++) {
00055 c[i4] = c[i4-1]*co - s[i4-1]*si;
00056 s[i4] = s[i4-1]*co + c[i4-1]*si;
00057 }
00058 for(i4 = n-1; i4 >= 1; i4--) m[i4-1] = m[i4]*2;
00059 }
00060
00061 /**** Main loop starts here ****/
00062
00063 for(i0 = 0; i0 < n; i0++) {
00064 i1 = 0;
00065 m0 = m[i0];
00066 for(i2 = 0; i2 < m[n-i0-1]; i2++) {
00067 f4 = c[i1];
00068 f5 = s[i1]*md;
00069 for(i4 = 2*m0*i2; i4 < m0*(2*i2+1); i4++) {
00070 f3 = xr[i4+m0]*f4 - xi[i4+m0]*f5;
00071 f1 = xi[i4+m0]*f4 + xr[i4+m0]*f5;
00072 xr[i4+m0] = xr[i4] - f3;
00073 xr[i4] += f3;
00074 xi[i4+m0] = xi[i4] - f1;
00075 xi[i4] += f1;
00076 }
00077 for(i4 = 1; i4 < n; i4++) {
00078 i5 = i4;
00079 if (i1 < m[i4]) goto i1_plus1;
00080 i1 -= m[i4];
00081 }
00082 i1_plus1: i1 += m[i5];
00083 }
00084 }
00085 i1 = 0;
00086 for(i4 = 0; i4 < id; i4++) {
00087 if (i1 > i4) {
00088 f3 = xr[i4];
00089 f1 = xi[i4];
00090 xr[i4] = xr[i1];
00091 xi[i4] = xi[i1];
00092 xr[i1] = f3;
00093 xi[i1] = f1;
00094 }
00095 for(i2 = 0; i2 < n; i2++) {
00096 i5 = i2;
00097 if (i1 < m[i2]) goto i1_plus2;
00098 i1 -= m[i2];
00099 }
00100 i1_plus2: i1 += m[i5];
00101 }
00102
00103 if (md > 0.) {
00104 f1 = 1./f2;
00105 for(i4 = 0; i4 < id; i4++) {
00106 xr[i4] *= f1;
00107 xi[i4] *= f1;
00108 }
00109 }
00110 return;
00111 }
|
|
||||||||||||||||||||||||
|
Definition at line 116 of file mri_cfft.c. References cfft(), free, and malloc. Referenced by mri_fft_complex().
00118 {
00119 float *rbuf , *ibuf ;
00120 register int ii , jj , jbase ;
00121
00122 rbuf = (float *)malloc( ny * sizeof(float) ) ;
00123 ibuf = (float *)malloc( ny * sizeof(float) ) ;
00124 if( rbuf == NULL || ibuf == NULL ){
00125 fprintf(stderr,"malloc error in cfft2d\n") ;
00126 exit(1) ;
00127 }
00128
00129 for( jj=0 ; jj < ny ; jj++ ){
00130 jbase = nx * jj ;
00131 cfft( mode , nx , &xr[jbase] , &xi[jbase] ) ;
00132 }
00133
00134 for( ii=0 ; ii < nx ; ii++ ){
00135 for( jj=0 ; jj < ny ; jj++ ){
00136 rbuf[jj] = xr[ii + jj*nx] ;
00137 ibuf[jj] = xi[ii + jj*nx] ;
00138 }
00139 cfft( mode , ny , rbuf , ibuf ) ;
00140 for( jj=0 ; jj < ny ; jj++ ){
00141 xr[ii+jj*nx] = rbuf[jj] ;
00142 xi[ii+jj*nx] = ibuf[jj] ;
00143 }
00144 }
00145
00146 free(rbuf) ; free(ibuf) ;
00147 return ;
00148 }
|
