00001 00005 #include <bpm/bpm_messages.h> 00006 #include <bpm/bpm_nr.h> 00007 00008 00032 int nr_four1( double data[], unsigned long nn, int isign ) { 00033 00034 unsigned long n,mmax,m,j,istep,i; 00035 double wtemp,wr,wpr,wpi,wi,theta; 00036 double tempr,tempi; 00037 00038 if ( ! nr_is_pow2( nn ) ) { 00039 bpm_error( "Data length is not power of 2 in nr_four1(...)", 00040 __FILE__, __LINE__ ); 00041 return BPM_FAILURE; 00042 } 00043 00044 n=nn<<1; 00045 j=1; 00046 for (i=1;i<n;i+=2) 00047 { 00048 if (j > i) 00049 { 00050 SWAP(data[j],data[i]); 00051 SWAP(data[j+1],data[i+1]); 00052 } 00053 00054 m=nn; 00055 while (m >= 2 && j > m) 00056 { 00057 j -= m; 00058 m >>= 1; 00059 } 00060 j +=m; 00061 } 00062 00063 mmax=2; 00064 while(n>mmax) 00065 { 00066 istep=mmax << 1; 00067 theta=isign*(2*PI/mmax); 00068 wtemp=sin(0.5*theta); 00069 wpr = -2.0*wtemp*wtemp; 00070 00071 wpi=sin(theta); 00072 wr=1.0; 00073 wi=0.0; 00074 00075 for (m=1;m<mmax;m+=2) 00076 { 00077 for (i=m;i<=n;i+=istep) 00078 { 00079 j=i+mmax; 00080 tempr=wr*data[j]-wi*data[j+1]; 00081 tempi=wr*data[j+1]+wi*data[j]; 00082 data[j]=data[i]-tempr; 00083 data[j+1]=data[i+1]-tempi; 00084 data[i] += tempr; 00085 data[i+1] += tempi; 00086 } 00087 wr=(wtemp=wr)*wpr-wi*wpi+wr; 00088 wi=wi*wpr+wtemp*wpi+wi; 00089 } 00090 mmax=istep; 00091 } 00092 00093 return BPM_SUCCESS; 00094 }