Zcosft2.cin

// zcosft2.cin is the C++ implementation of the DCTII and DCTIII

// For the use, files zfour1.cin and zrealft.cin also should be loaded.

// y is the input array, that will be replaced with the DCTII or DCTIII transform; in the array y-1, the numeration of elements begins with zero.

// n is number of elements; \(n=2^q\) for some natural \(q\).

// isign should be -1 for DCTII and 1 for  DCTIII

void zcosft2(z_type y[], int n, int isign) { // void zrealft(z_type data[], unsigned long n, int isign); int i; z_type sum,sum1,y1,y2,ytemp; double theta,wi=0.0,wi1,wpi,wpr,wr=1.0,wr1,wtemp; theta=0.5*M_PI/n; wr1=cos(theta); wi1=sin(theta); wpr = -2.0*wi1*wi1; wpi=sin(2.0*theta); if(isign == 1) { for (i=1;i<=n/2;i++) { y1=0.5*(y[i]+y[n-i+1]); y2=wi1*(y[i]-y[n-i+1]); y[i]=y1+y2; y[n-i+1]=y1-y2; wr1=(wtemp=wr1)*wpr-wi1*wpi+wr1; wi1=wi1*wpr+wtemp*wpi+wi1; }  zrealft(y,n,1); for (i=3;i<=n;i+=2) {    wr=(wtemp=wr)*wpr-wi*wpi+wr; wi=wi*wpr+wtemp*wpi+wi; y1=y[i]*wr-y[i+1]*wi; y2=y[i+1]*wr+y[i]*wi; y[i]=y1; y[i+1]=y2; } sum=0.5*y[2]; for (i=n;i>=2;i-=2) {   sum1=sum; sum += y[i]; y[i]=sum1; } } else if (isign == -1) {    ytemp=y[n]; for (i=n;i>=4;i-=2) y[i]=y[i-2]-y[i]; y[2]=2.0*ytemp; for (i=3;i<=n;i+=2) { wr=(wtemp=wr)*wpr-wi*wpi+wr; wi=wi*wpr+wtemp*wpi+wi; y1=y[i]*wr+y[i+1]*wi; y2=y[i+1]*wr-y[i]*wi; y[i]=y1; y[i+1]=y2; }  zrealft(y,n,-1); for(i=1;i<=n/2;i++) { y1=y[i]+y[n-i+1]; y2=(0.5/wi1)*(y[i]-y[n-i+1]); y[i]=0.5*(y1+y2); y[n-i+1]=0.5*(y1-y2); wr1=(wtemp=wr1)*wpr-wi1*wpi+wr1; wi1=wi1*wpr+wtemp*wpi+wi1; }   } }

Comments
"z_type" should be defined as "double" in order to deal with real objects, and as "complex(double)" in odder to transform the complex arrays. However, the transformation of more complicated objects can be performed with the same routines, at the appropriate definition of "z_type".

Articles DCTII and DCTIII provide the examples of call zcosft2.

Keywords
DCT, DCTII, DCTIII, C++