ＤＣＴ（離散コサイン変換）プログラム

/*		dct.c													*/
/*																*/
/*	DCT : Discrete Cosine Transformation （離散コサイン変換）	*/
/*																*/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#define		N		32

int dct1(double x[], double y[], int n, int flag)
{
	int i, j, k, ks, n2, n4;
	double sigma, sq2, nrm, t, data;
	double *cosTable;

	cosTable = (double *)malloc(4 * n * sizeof(double));
	if(cosTable == NULL)
	{
		fprintf(stderr, "Error : out of memory in dct1.\n");
		return 999;
	}

	n2 = n * 2;
	n4 = n * 4;
	t = M_PI / 2. / n;
	cosTable[0] = 1.;
	for(i = 1; i < n; i++)	cosTable[i] = cos(t * i);
	cosTable[i++] = 0.;
	for(j = n - 1; i < n2; i++)	cosTable[i] = -cosTable[j--];
	cosTable[i++] = -1.;
	for(j = n2 - 1; i < n4; i++)	cosTable[i] = cosTable[j--]; 

	nrm = sqrt(2.0 / n);
	sq2 = sqrt(2.0) / 2.;

	if(flag == 0)
	{
		ks = 0;
		for(j = 0; j < n; j++, ks += 2)
		{
			sigma = 0.;
			k = j;

			for(i = 0; i < n; i++)
			{
				sigma += x[i] * cosTable[k] ;
				k = (k + ks) % n4;
			}
			y[j] = nrm * sigma ;
		}
		y[0] *= sq2;
	}
	else
	{
		ks = 1;
		for(j = 0; j < n; j++, ks += 2)
		{
			sigma = x[0] * cosTable[0] * sq2;
			k = ks;

			for(i = 1; i < n; i++)
			{
				sigma += x[i] * cosTable[k];
				k = (k + ks) % n4;
			}
			y[j] = nrm * sigma ;
		}
	}

	free((char *)cosTable);
	return 1;
}

int main(void)
{
	double data[N];
	double *a1, *a2, *b, x, y;
	int i;

	a1 = (double *)calloc(N, sizeof(double));
	a2 = (double *)calloc(N, sizeof(double));
	b = (double *)calloc(N, sizeof(double));
	if(a1 == NULL || a2 == NULL || b == NULL)
	{
		fprintf(stderr, "Error : out of memory.\n");
		return 0;
	}

	for(i = 0; i < N; i++)	a1[i] = sin(2. * M_PI * i / N);

	dct1(a1, b, N, 0);
	dct1(b, a2, N, 1);

	for(i = 0; i < N; i++)
	{
		x = 20. * log10(fabs(b[i]));
		y = a1[i] - a2[i];
		printf("%2d : %8.4f %e\n", i, x, y);
	}

	free((char *)a1);
	free((char *)a2);
	free((char *)b);
	return 1;
}

 0 : -308.0148 -2.164118e-16
 1 :  10.5861 -1.942890e-16
 2 :  -8.1328 -2.220446e-16
 3 :   6.1211 2.220446e-16
 4 : -323.4714 -1.110223e-16
 5 :  -6.4015 -1.110223e-16
 6 : -332.4222 -3.330669e-16
 7 : -13.1115 -1.110223e-16
 8 : -327.0282 0.000000e+00
 9 : -17.9046 -1.110223e-16
10 : -315.9273 -1.110223e-16
11 : -21.7105 0.000000e+00
12 : -337.8625 -2.220446e-16
13 : -24.9232 -2.220446e-16
14 : -317.5464 -1.110223e-16
15 : -27.7576 -2.498002e-16
16 : -325.1496 -3.002430e-17
17 : -30.3536 -1.110223e-16
18 : -313.0131 -5.551115e-17
19 : -32.8193 -2.220446e-16
20 : -313.2597 0.000000e+00
21 : -35.2559 -2.220446e-16
22 : -314.4256 0.000000e+00
23 : -37.7811 -2.220446e-16
24 : -316.1316 2.220446e-16
25 : -40.5663 1.110223e-16
26 : -319.2635 2.220446e-16
27 : -43.9280 -1.110223e-16
28 : -321.2751 -1.110223e-16
29 : -48.6519 1.110223e-16
30 : -327.5641 0.000000e+00
31 : -58.3362 -5.551115e-17