#define N 9
#define M 7
#define PI 3.141592654
#define COEFF 0.707106781

/***************************************************************************
 * Func: dct                                                               *
 *                                                                         *
 * Desc: forward Discrete Cosine Transform on MxN array                    *
 *                                                                         *
 * Params: forward - MxN array of data to transform                        *
 *         DCT - MxN array of trasformed data                              *
 ***************************************************************************/
void dct(double forward[M][N], double DCT[M][N])
    {
    int u, v;         /* frequency domain variables */
    int x, y;         /* spatial domain variables */
    double accum;     /* accumulator */
    double twoN;      /* 2 times N */
    double twoM;      /* 2 times M */
    double scale;     /* 2 / square_root(MN) */

    twoN = 2.0 * N;
    twoM = 2.0 * M;
    scale = (2.0 / sqrt((double)M*N));

    for(u=0; u<M; u++)
	for(v=0; v<N; v++)
	    {
	    accum = 0.0;
	    for(x=0; x<M; x++)
		for(y=0; y<N; y++)
		    {
		    accum += (cos((PI*u*(2*x+1))/twoM) * 
		    cos((PI*v*(2*y+1))/twoN) * forward[x][y]);
		    }
		accum *= scale;
		if(u == 0)
		    accum *= COEFF;
		if(v == 0)
		    accum *= COEFF;
		DCT[u][v] = accum;
	 }
    }