/*******
 * dct_transform.c
 *
 * Looking at the discrete cosine transform
 *
 *******/

#include <stdio.h>
#include <math.h>

// Number of points.  JPEG's typically used 8x8 2-dim DCT's.
#define N 8

// ------------------------------------------------------
// Define the 64 terms in the DCT and DCT_inverse matrix.
void defineDCT(double dct[N][N]){
  double pi = 4.0*atan(1.0);	// typical way to get pi, using tan(pi/4)=1
  double piOver2N = pi/(2*N);
  double invSqrt2 = 1.0/sqrt(2.0);
  double sqrt2overN = sqrt(2.0/N);
  double C;
  int u,x;
  for (u=0; u<N; u++){
    C = u==0 ? invSqrt2 : 1;
    for (x=0; x<N; x++) {
      dct[u][x] = sqrt2overN*C*cos( (2*x+1)*u*piOver2N );
    }
  }
}



// ------------------------------------------------------
// The DCT matrix is orthogonal, so its transpose is also its inverse.
void transposeMatrix(double a[N][N], double aTrans[N][N]){
  int row,col;
  for (row=0; row<N; row++){
    for (col=0; col<N; col++){
      aTrans[col][row] = a[row][col];
    }
  }
}


// ------------------------------------------------------
void matrixTimesVector(double m[N][N], double v[N], double answer[N]){
  int row,col;
  double tmp;
  for (row=0; row<N; row++){
    tmp=0;
    for (col=0; col<N; col++){
      tmp += m[row][col]*v[col];
    }
    answer[row]=tmp;
  }
}

// ------------------------------------------------------
int main(){

  double dct[N][N];		// All these are [row][column]
  double dctInv[N][N];

  double f[N] = { -.2, -.1, 0.0, 0.1, 0.2, 0.2, 0.2, 0.2 };
  double S[N];
  double ff[N];
  int i;

  defineDCT(dct);			
  transposeMatrix(dct, dctInv);		

  matrixTimesVector(dct, f, S);		// transform f -> S
  matrixTimesVector(dctInv, S, ff);	// inverse transform S -> ff

  printf("# Here's an array and its discrete cosine transform. \n");
  printf("# Using N=%d.\n",N);
  printf("#\n");
  printf("#  f[x] original,   S[u] transform,   ff[x] restored original \n");
  for (i=0; i<N; i++){
    printf("  % 7.5f     % 7.5f     % 7.5f \n", f[i], S[i], ff[i]);
  }
  

}