Realft.cin

// realft.cin is conversion of routine zrealft.cin from complex double to float, // as I could not access the original. I hope, it is close to the original from the Numerical recipes and does exactly the same.

void realft(float data[], unsigned long n, int isign)
{ /* void zfour1(z_type data[], unsigned long nn, int isign);*/
        unsigned long i,i1,i2,i3,i4,np3;
        float c1=0.5,c2,h1r,h1i,h2r,h2i; //z_type
        float wr,wi,wpr,wpi,wtemp,theta; //double

        theta=M_PI/(float) (n>>1);
        if (isign == 1) {
                c2 = -(float)0.5;
                four1(data,n>>1,1);
        } else {
                c2=(float)0.5;
                theta = -theta;
        }
        wtemp=sin((float)0.5*theta);
        wpr = -(float)2.0*wtemp*wtemp;
        wpi=sin(theta);
        wr=(float)1.0+wpr;
        wi=wpi;
        np3=n+3;
        for (i=2;i<=(n>>2);i++) {
                i4=1+(i3=np3-(i2=1+(i1=i+i-1)));
                h1r=c1*(data[i1]+data[i3]);
                h1i=c1*(data[i2]-data[i4]);
                h2r = -c2*(data[i2]+data[i4]);
                h2i=c2*(data[i1]-data[i3]);
                data[i1]=h1r+wr*h2r-wi*h2i;
                data[i2]=h1i+wr*h2i+wi*h2r;
                data[i3]=h1r-wr*h2r+wi*h2i;
                data[i4] = -h1i+wr*h2i+wi*h2r;
                wr=(wtemp=wr)*wpr-wi*wpi+wr;
                wi=wi*wpr+wtemp*wpi+wi;
        }
        if (isign == 1) {
                data[1] = (h1r=data[1])+data[2];
                data[2] = h1r-data[2];
        } else {
                data[1]=c1*((h1r=data[1])+data[2]);
                data[2]=c1*(h1r-data[2]);
                four1(data,n>>1,-1);
        }
}

Keywords

Fourier transform, sinft, DSFT]], cosft, DCFT]]