00001 00005 #include <bpm/bpm_messages.h> 00006 #include <bpm/bpm_nr.h> 00007 00008 00027 int nr_realft( double data[], unsigned long n, int isign ) { 00028 00029 00030 unsigned long i, i1,i2,i3,i4,np3; 00031 double c1=0.5,c2,h1r,h1i,h2r,h2i; 00032 double wr,wi,wpr,wpi,wtemp,theta; 00033 00034 if ( ! nr_is_pow2( n ) ) { 00035 bpm_error( "Data length is not power of 2 in nr_realft(...)", 00036 __FILE__, __LINE__ ); 00037 return BPM_FAILURE; 00038 } 00039 00040 theta = PI/(double)(n>>1); 00041 00042 if (isign == 1){ 00043 c2 = -0.5; 00044 nr_four1(data,n>>1,1); 00045 } else { 00046 c2=0.5; 00047 theta = -theta; 00048 } 00049 00050 wtemp=sin(0.5*theta); 00051 wpr = -2.0*wtemp*wtemp; 00052 wpi=sin(theta); 00053 wr=1.0+wpr; 00054 wi=wpi; 00055 np3=n+3; 00056 for(i=2;i<=(n>>2);i++) 00057 { 00058 i4=1+(i3=np3-(i2=1+(i1=i+i-1))); 00059 h1r=c1*(data[i1]+data[i3]); 00060 h1i=c1*(data[i2]-data[i4]); 00061 h2r = -c2*(data[i2]+data[i4]); 00062 h2i=c2*(data[i1]-data[i3]); 00063 data[i1]=h1r+wr*h2r-wi*h2i; 00064 data[i2]=h1i+wr*h2i+wi*h2r; 00065 data[i3]=h1r-wr*h2r+wi*h2i; 00066 data[i4] = -h1i+wr*h2i+wi*h2r; 00067 wr=(wtemp=wr)*wpr-wi*wpi+wr; 00068 wi=wi*wpr+wtemp*wpi+wi; 00069 } 00070 if (isign == 1) 00071 { 00072 data[1] = (h1r=data[1])+data[2]; 00073 data[2] = h1r-data[2]; 00074 } 00075 else 00076 { 00077 data[1]=c1*((h1r=data[1])+data[2]); 00078 data[2]=c1*(h1r-data[2]); 00079 00080 nr_four1(data,n>>1,-1); 00081 } 00082 00083 return BPM_SUCCESS; 00084 } 00085