DCT3回目開始！まずは、公式から

wikipedia等によると、公式はいくつかあるようなのですが、私は下記の公式を使う事にしました。
これから、以下の公式通りにプログラムを組んでみて、変換⇒逆変換と動かしてみて、
この公式が本当に正しいのかどうかを確かめてみようと思います。

===== 公式 =========================================================

N:変換する対象の要素数を表します。
π:円周率のパイです(3.1415926535・・・・)
S(x):DCT変換する前の値です。要素数N個の1次元の配列です。
S(u):DCT変換後の値です。要素数N個の1次元の配列です。
この値はDCT係数と呼ばれ、各々の波が、どの程度含まれているのか？という事を示しているそうです。

///// 変換 /////////////////////////////////////////////////////////

                        N-1                            
S(u) = √2/√N * C(u) * Σ  S(x) * cos( (2x+1)u*π/2N)  
                        x=0                            
但し、C(u)は、
    u=0  の時、1/√2
    u≠0 の時、1

///// 逆変換 ///////////////////////////////////////////////////////

                 N-1                                  
S(x) = √2/√N * Σ  C(u) * S(u) * cos( (2x+1)u*π/2N) 
                 u=0                                  
但し、C(u)は、
    u=0  の時、1/√2
    u≠0 の時、1     です。

====================================================================

↑が公式です。

今回は、要素の数を8個として実験します。
つまり、高さ1、幅8 の画像を1次元DCTするというイメージで行きたいと思います。

@	0	1	2	3	4	5	6	7
	10	20	30	40	50	60	70	80

次に、要素数を８としたときの公式と、そのプログラムを書きます。
上の公式の N が 8 になっただけです。

 /*
======変換=========================================
                     7
S(u) = 1/2 * C(u) * Σ  * S(x) * cos( (2x+1)uπ/16 )
                    x=0
但し、C(u)は、
    u=0  の時、1/√2
    u≠0 の時、1     です。
===================================================
*/
 double Cu = 0.0;
 double sigma = 0.0;
 double sigmain = 0.0;

 double cos=0.0;
 double cosin = 0.0;
 sb->AppendFormat("-----変換-----\r\n");

 for(u=0; u<8; u++){
 sigma = 0;
	 for(x=0; x<8; x++){
		 cosin = (2*x+1)*u*Math::PI/16;
		 cos = Math::Cos( cosin );
		 sigmain = Sx[x] * cos;
		 sigma += sigmain;
		 sb->AppendFormat("Sx[{4}]={5,4}, cos({0,8:F4})={1,8:f4},\tsigmain={2,10:f4},\tsigma={3,10:f4}\r\n"
			,cosin,cos,sigmain,sigma,x,Sx[x]);
	 }
	 if(u==0){
		 Cu = Math::Sqrt(1.0/2.0);
	 }else{
		 Cu = 1.0;
	 }
	 Su[u] = (int)(sigma/2.0*Cu);
	 sb->AppendFormat("Su[{0}]={1}\r\n",u,Su[u]);
 }

/*
=====逆変換========================================
              7
S(x) = 1/2 * Σ  * C(u) * S(u) * cos( (2x+1)uπ/16)
             u=0
但し、C(u)は、
	u=0  の時、1/√2
	u≠0 の時、1     です。
===================================================
*/
 sb->AppendFormat("-----逆変換-----\r\n");

 for(x=0; x<8; x++){
	 sigma = 0;

	 for(u=0; u<8; u++){
		 if(u==0){
			 Cu = Math::Sqrt(1.0/2.0);
		 }else{
			 Cu = 1.0;
		 }
		 cosin = (2*x+1)*u*Math::PI/16;
		 cos = Math::Cos(cosin);
		 sigmain = Su[u]*Cu*cos;
		 sigma += sigmain;
		 sb->AppendFormat("Su[{4}]={5,4}, cos({0,8:F4})={1,8:f4},\tsigmain={2,10:f4},\tsigma={3,10:f4}\r\n"
			,cosin,cos,sigmain,sigma,u,Su[u]);
	 }
	 Sx2[x] = (int)sigma/2;	 
	 sb->AppendFormat("Sx2[{0}]={1}\r\n",x,Sx2[x]);
 }
			 
 textBox3->Text = sb->ToString();

※計算途中の値を出力したいので、計算は、数回に分けてあります。

これを実行すると、以下のように出力され、この公式が、おそらく正しいであろうと思われるところまで来ました。

※double型からint型にキャスト(小数点切捨て)している為に、誤差が出ています。
また、DCT⇒逆DCTのみなら、プログラム上の丸め処理を考慮に入れなければ、可逆変換(完全に元に戻せる)です。

入力値	1	10	100	200	300	900	10	1
DCT後	538	-231	-417	420	-183	-119	357	-283
逆DCT後	1	10	100	199	299	899	10	0

入力値	100	100	100	100	100	100	100	100
DCT後	282	0	0	0	0	0	0	0
逆DCT後	99	99	99	99	99	99	99	99