CでDCTと逆DCTをしてみた。

DCTⅡのDCTと逆DCTを素直にCで実装。

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

#define PI     3.14159265358979323846
#define K_MAX 64
#define N_MAX 64

void dct(int * x, int * X, const int N, const int K)
{
    int n, k;
    double C[K_MAX][N_MAX];
    double Scale;

    Scale = sqrt(2.0 / N);
    for (k = 0; k < K; k++) {
        for (n = 0; n < N; n++) {
            C[k][n] = Scale * cos((n + 0.5) * k * PI / N);
            if (k == 0) {
                C[0][n] *= 1.0 / sqrt(2.0);
            }
        }
    }
        
    for (k = 0; k < K; k++) {
        X[k] = 0;
        for (n = 0; n < N; n++) {
            X[k] += C[k][n] * x[n];     
        }
    }
}

void idct(int * x, int * X, const int N, const int K)
{
    int n, k;
    double C[K_MAX][N_MAX];
    double Scale;

    Scale = sqrt(2.0 / N);
    for (k = 0; k < K; k++) {
        for (n = 0; n < N; n++) {
            C[k][n] = Scale * cos((2*n + 1) * k * PI / (2 * N));
        }
    }

    for (n = 0; n < N; n++) {
        X[n] = 0;
        for (k = 0; k < K; k++) {
            X[n] += C[k][n] * x[k];     
        }
    }
}

適当な三角波を生成して、N=64 K=32でDCTと逆DCTをしてみた。

int main(void)
{
   int i,j,k,n;
   const int N = 64;
   const int K = 32;
   int a_time[64];
   int a_spec[64];
   int a_out[64];

   for (i = 0; i < 64; i++) {
       if( i < 31)      
           a_time[i] = i * 256  ;
       else
           a_time[i] = 8192 - (i-31) * 256 ;
   }

   printf("\n a_time");
   for (i = 0; i < N; i++) {
       printf(" %d",a_time[i]);
   }
   printf("\n\n");

   dct(&a_time[0], a_spec, N, K);

   printf("\n a_spec");
   for (i = 0; i < K; i++) {
       printf(" %d",a_spec[i]);
   }
   printf("\n\n");

   idct(a_spec, a_out, N, K);

   printf("\n out");
   for (i = 0; i < N; i++) {
       printf(" %d",a_out[i]);
   }
   printf("\n\n");

   return 0;
}

実行結果

DCTをするとN次元ベクトルがK次元ベクトルになる。

逆DCTをするとK次元ベクトルがN次元ベクトルになる。

なんとなく、それっぽい波形が復元できている。

端の方は誤差が多めだね。

 a_time 0 256 512 768 1024 1280 1536 1792 2048 2304 2560 2816 3072 3328 3584 3840 4096 4352 4608 4864 5120 5376 5632 5888 6144 6400 6656 6912 7168 7424 7680 8192 7936 7680 7424 7168 6912 6656 6400 6144 5888 5632 5376 5120 4864 4608 4352 4096 3840 3584 3328 3072 2816 2560 2304 2048 1792 1536 1280 1024 768 512 256 0

 a_spec 31776 -29 -18811 0 43 7 -2115 -3 41 13 -778 -9 45 15 -406 -12 35 21 -260 -17 39 24 -180 -19 35 28 -138 -24 37 30 -102 -27

 out 1655 1863 2146 2401 2636 2900 3167 3417 3663 3929 4195 4445 4690 4953 5217 5463 5711 5978 6247 6493 6734 7001 7275 7515 7749 8024 8313 8544 8757 9050 9413 9663 9620 9338 9021 8784 8565 8298 8019 7766 7534 7272 7000 6749 6508 6250 5981 5722 5478 5222 4960 4703 4462 4202 3934 3678 3434 3176 2903 2653 2420 2153 1864 1652

Kを4にしてもそれなりの波形が復元できるのが、すごいね。

 a_time 0 256 512 768 1024 1280 1536 1792 2048 2304 2560 2816 3072 3328 3584 3840 4096 4352 4608 4864 5120 5376 5632 5888 6144 6400 6656 6912 7168 7424 7680 8192 7936 7680 7424 7168 6912 6656 6400 6144 5888 5632 5376 5120 4864 4608 4352 4096 3840 3584 3328 3072 2816 2560 2304 2048 1792 1536 1280 1024 768 512 256 0

 a_spec 31776 -29 -18811 0

 out 2289 2321 2385 2480 2604 2759 2941 3148 3378 3631 3902 4190 4491 4804 5125 5449 5776 6100 6420 6734 7035 7323 7594 7847 8078 8285 8467 8621 8746 8841 8905 8937 8938 8906 8842 8747 8624 8470 8288 8081 7852 7599 7328 7040 6739 6427 6107 5783 5456 5132 4813 4500 4199 3911 3640 3387 3157 2950 2768 2615 2491 2396 2332 2300

参考

貴家仁志,村松正吾マルチメディア技術の基礎 DCT(離散コサイン変換)入門 (Amazon.co.jp)