Cosft.cin

// cosft.cin is the C++ numerical implementation of the DiscreteCos transform.

// The input is array; z_type is supposed to be defined as complex(double), although other definition (for example, just double or float) also may have sense for some applications; in the old book Numerical recipes in C, the argument is supposed to be array of float variables.

#define SWAP(a,b) tempr=(a);(a)=(b);(b)=tempr
void zfour1(z_type data[], unsigned long nn, int isign)
{
unsigned long n,mmax,m,j,istep,i;
double wtemp,wr,wpr,wpi,wi,theta;
z_type tempr,tempi;
n=nn << 1;
j=1;
for(i=1;i<n;i+=2){
                       if(j>i){SWAP(data[j],data[i]);
                               SWAP(data[j+1],data[i+1]); }
                       m=n >> 1;
                       while (m >= 2 && j > m) { j -= m; m >>= 1; }
                       j += m;
                }
mmax=2;
while(n>mmax)
{ istep=mmax << 1;
  theta=isign*(6.28318530717959/mmax);
  wtemp=sin(0.5*theta);
  wpr = -2.0*wtemp*wtemp;
  wpi=sin(theta);
  wr=1.0;
  wi=0.0;
  for(m=1;m<mmax;m+=2)
       {for (i=m;i<=n;i+=istep)
          {j=i+mmax;
           tempr=wr*data[j]-wi*data[j+1];
           tempi=wr*data[j+1]+wi*data[j];
           data[j]=data[i]-tempr;
           data[j+1]=data[i+1]-tempi;
           data[i] += tempr;
           data[i+1] += tempi;
          }
        wr=(wtemp=wr)*wpr-wi*wpi+wr;
        wi=wi*wpr+wtemp*wpi+wi;
       }
  mmax=istep;
}
}
#undef SWAP

/* #include <math.h>*/
void zrealft(z_type data[], unsigned long n, int isign)
{ /* void zfour1(z_type data[], unsigned long nn, int isign);*/
 unsigned long i,i1,i2,i3,i4,np3;
 z_type c1=0.5,c2,h1r,h1i,h2r,h2i;
 double wr,wi,wpr,wpi,wtemp,theta;
 theta=M_PI/(double) (n>>1);
 if(isign == 1){ c2 = -(double)0.5; zfour1(data,n>>1,1);} 
          else { c2 =  (double)0.5; theta = -theta; }
 wtemp=sin((double)0.5*theta);
 wpr = -(double)2.0*wtemp*wtemp;
 wpi=sin(theta);
 wr=(double)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]);
                zfour1(data,n>>1,-1); }
}

void zcosft1(z_type y[], int n)
{ /* void zrealft(z_type data[], unsigned long n, int isign);*/
       int j,n2;  z_type sum,y1,y2;
       double theta,wi=0.0,wpi,wpr,wr=1.0,wtemp;
       theta=M_PI/n;
       wtemp=sin(0.5*theta);
       wpr = -2.0*wtemp*wtemp;
       wpi=sin(theta);
       sum=0.5*(y[1]-y[n+1]);
       y[1]=0.5*(y[1]+y[n+1]);
       n2=n+2;
       for (j=2;j<=(n>>1);j++) {
               wr=(wtemp=wr)*wpr-wi*wpi+wr;
               wi=wi*wpr+wtemp*wpi+wi;
               y1=0.5*(y[j]+y[n2-j]);
               y2=(y[j]-y[n2-j]);
               y[j]=y1-wi*y2;
               y[n2-j]=y1+wi*y2;
               sum += wr*y2;
                                              }
       zrealft(y,n,1);
       y[n+1]=y[2];
       y[2]=sum;
       for(j=4;j<=n;j+=2) {sum += y[j]; y[j]=sum;}
}

void cosft(z_type a[], int N){ int n; DB d; zcosft1(a-1,N); d=sqrt(2./N);  DO(n,N) a[n]*=d; }

Keywords

DiscreteCos, Fourier operator