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.

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;ii){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*/ 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