bpmnr/nr_realft.c

Go to the documentation of this file.
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     

Generated on Thu Jan 10 10:18:04 2008 for libbpm by  doxygen 1.5.1