第６章 フーリエ変換

【リスト6-1】 複素数計算用関数
typedef struct{double r; double i; } Complex;
Complex complex(double r, double i){
	Complex X; X.r=r; X.i=i; return X;
}
Complex addC(Complex A, Complex B){
	return complex(A.r+B.r, A.i+B.i);
}
Complex subC(Complex A, Complex B){
	return complex(A.r-B.r, A.i-B.i);
}
Complex multC(Complex A, Complex B){
	return complex(A.r*B.r-A.i*B.i, A.r*B.i+A.i*B.r);
}
Complex multR(Complex A, double B){
	return complex(A.r*B, A.i*B);
}
Complex expJ(double A){// exp(θj)の設定
	return complex(cos(A), sin(A));
}
Complex divR(Complex A, double B){
	return complex(A.r/B, A.i/B);
}

(C++用オペレーション定義）  
  
  Complex operator +(Complex x, Complex y){ return addC(x, y);}
  Complex operator -(Complex x, Complex y){ return subC(x, y);}
  Complex operator *(Complex x, Complex y){ return multC(x, y);}
  Complex operator *(Complex x, double  y){ return multR(x, y);}
  Complex operator /(Complex x, double  y){ return divR(x, y);}
  
 
【リスト6-2】 DFTと逆DFT
#include "myWin.h"
#include "math.h"
#define XST 420
#define YST 10
#define GWIDTH 400
#define GHIGHT 400
#define RED(X) (X & 0xFF)
#define GREEN(X) ((X>>8) & 0xFF)
#define BLUE(X) ((X>>16) & 0xFF)
static HBITMAP hBM[2];static HDC hBuff[2]; 
static BITMAP bM[2];
void clearBitmap(int i){//ビットマップのクリア
	SelectObject(hBuff[i], GetStockObject(NULL_PEN));
	PatBlt(hBuff[i], 0,0,GWIDTH,GHIGHT,WHITENESS);//黒塗り
}
void createBitmap(HWND hw){//ビットマップの設定
	HDC hdc=GetDC(hw);
	for(int i=0;i<2;i++){
		hBM[i]=CreateCompatibleBitmap(hdc,GWIDTH,GHIGHT);
		hBuff[i]=CreateCompatibleDC(hdc); SelectObject(hBuff[i], hBM[i]);
		clearBitmap(i);
	}
}
int pixel(int i, int j){//画像1からピクセルを取り出す
	return GetPixel(hBuff[1], i,j);
}
void setPixel(int i, int j, COLORREF C){//画像0にピクセルを設定する
	SetPixel(hBuff[0], i,j, C);
}

------（複素数演算用関数省略: 【リスト6-1】参照）-----

#define GWDFT 80
#define GHDFT 80
Complex beforDFT[GWDFT+1][GHDFT+1], tmDFT[GWDFT+1][GHDFT+1];
Complex afterDFT[GWDFT+1][GHDFT+1],reverceDFT[GWDFT+1][GHDFT+1];
//グレースケール画像にして１以下の値にする
void setComplexDFT(){
	for(int i=0;i<GWDFT;i++){
	for(int j=0;j<GHDFT;j++){
			COLORREF C=pixel(i+80,j+80),R=RED(C), G=GREEN(C), B=BLUE(C);
			beforDFT[i][j]=complex((double)(R+G+B)/255/3, 0);
		}
		beforDFT[i][GHDFT]=beforDFT[i][0];
	}
	beforDFT[GWDFT][GHDFT]=beforDFT[0][0];
}
//DFT
void DFT(Complex beforDFT[GWDFT+1][GHDFT+1],
				 Complex afterDFT[GWDFT+1][GHDFT+1]){
	double ARG=-2.0*PI/GHDFT;
	for(int k=0;k<=GWDFT;k++)for(int N=0;N<=GWDFT;N++){//中央への移動
		double DD=((k+N)%2)!=0 ? -1.0: 1.0;
		afterDFT[k][N]=multR(beforDFT[k][N],DD);
	}
	for(int i=0;i<GWDFT;i++){//Y方向変換
		for(int k=0;k<GHDFT;k++){
			tmDFT[i][k]=complex(0,0);
			for(int N=0;N<GHDFT;N++){
				tmDFT[i][k]=addC(tmDFT[i][k],
							 multC(afterDFT[i][N], expJ(ARG*N*k)));
			}
		}
		tmDFT[i][GHDFT]=tmDFT[i][0];
	}
	tmDFT[GWDFT][GHDFT]=tmDFT[0][0];
	ARG=-2.0*PI/GWDFT;
	for(int i=0;i<GHDFT;i++){//X方向変換
		for(int k=0;k<GWDFT;k++){
			afterDFT[k][i]=complex(0,0);
			for(int N=0;N<GWDFT;N++){
				afterDFT[k][i]=addC(afterDFT[k][i],
							 multC(tmDFT[N][i], expJ(ARG*N*k)));
			}
		}
		tmDFT[i][GHDFT]=tmDFT[i][0];
	}
}
//逆DFT
void IDFT(Complex afterDFT[GWDFT+1][GHDFT+1],
						 Complex reverceDFT[GWDFT+1][GHDFT+1]){
	double ARG=2.0*PI/GHDFT; //Y方向変換
	for(int i=0;i<GWDFT;i++){
		for(int k=0;k<GHDFT;k++){
			reverceDFT[i][k]=complex(0,0);
			for(int N=0;N<GHDFT;N++){
				reverceDFT[i][k]=addC(reverceDFT[i][k],
							 multC(afterDFT[i][N], expJ(ARG*N*k)));
			}
		}
		reverceDFT[i][GHDFT]=reverceDFT[i][0];
	}
	reverceDFT[GWDFT][GHDFT]=reverceDFT[0][0];
	for(int i=0;i<=GWDFT;i++){
		for(int k=0;k<=GHDFT;k++){
			reverceDFT[i][k]=divR(reverceDFT[i][k],GHDFT);
		}
	}
	ARG=2.0*PI/GWDFT; //X方向変換
	for(int i=0;i<GHDFT;i++){ 
		for(int k=0;k<GWDFT;k++){
			tmDFT[k][i]=complex(0,0);
			for(int N=0;N<GWDFT;N++){
				tmDFT[k][i]=addC(tmDFT[k][i],
							 multC(reverceDFT[N][i], expJ(ARG*N*k)));
			}
		}
		tmDFT[GWDFT][i]=tmDFT[0][i];
	}
	for(int i=0;i<=GWDFT;i++){
		for(int k=0;k<=GHDFT;k++){
			tmDFT[i][k]=divR(tmDFT[i][k],GHDFT);
		}
	}
	for(int k=0;k<=GWDFT;k++)for(int N=0;N<=GWDFT;N++){//周辺への移動
		double DD=((k+N)%2)!=0 ? -1.0: 1.0;
		reverceDFT[k][N]=multR(tmDFT[k][N],DD);
	}
}
//DFT結果のスペクトル表示（急激な変化を対数値で緩和する）
double spec[GWDFT][GHDFT];
void dspSpec(Complex DT[GWDFT+1][GHDFT+1], 
						int startX, int startY, int flag ){
	for (int i = 0; i < GWDFT; ++i) for (int j = 0; j < GHDFT; ++j){
		double DR=DT[i][j].r,  DI=DT[i][j].i;
		double IR=(DR>=0 ? 1: -1), II=(DI>=0 ? 1: -1);
		DR=abs(DR); DI=abs(DI); 
		switch(flag){
		case 0: spec[i][j] = DI*log(DR+1)/log(10.0);break;//実数部
		case 1: spec[i][j] = II*log(DI+1)/log(10.0);break;//虚数部
		default: spec[i][j] = DI*II*log(DR+DI+1)/log(10.0);//絶対値
		}
	}
	double max=spec[0][0], min=spec[0][0]; //最大値と最小値を求める
	for (int i = 0; i < GWDFT; ++i) for (int j = 0; j < GHDFT; ++j) {
		if(spec[i][j] > max) max = spec[i][j];
     if(spec[i][j] < min) min = spec[i][j];
	}
	for (int i = 0; i < GWDFT; ++i) for (int j = 0; j < GHDFT; ++j){
		int C = (int)(((spec[i][j] - min)/(max - min)*255)+0.5);//比例配分
		if(C > 255) C = 255; if(C<0)C=0;
        setPixel(startX+i, startY+j, RGB(C,C,C));
	}
}
//0〜1の値を4種類のいずれかで表示(RED, GREEN, BLUE, グレースケール)
void dspValue(Complex DT[GWDFT+1][GHDFT+1], 
						int startX, int startY, int flag){
	for (int i = 0; i < GWDFT; i++) {
     for (int j = 0; j < GHDFT; j++) {
			int C=(int)(DT[i][j].r*255+0.5);if(C>255)C=255;
			switch(flag){
			case 1 : C=RGB(C, 0, 0);break;//RED表示
			case 2 : C=RGB(0, C, 0);break;//GREEN表示
			case 3 : C=RGB(0, 0, C);break;//BLUE表示
			default: C=RGB(C, C, C);break;//グレースケール表示
			}
			setPixel(startX+i, startY+j, C);
		}
	}
}
//テスト用メイン
void procDFT(){
	setComplexDFT();  dspValue(beforDFT, 0,0, 0);
	DFT(beforDFT,afterDFT);
	dspSpec(afterDFT, 100, 0, 0 ); dspSpec(afterDFT, 200, 0, 1 );
	IDFT(afterDFT,reverceDFT);  dspValue(reverceDFT, 300,0,0);
}
void procLButtonDown(HWND hw, WPARAM wp,LPARAM lp){
   procDFT();
   InvalidateRect(hw,NULL,TRUE);
}
void procRButtonDown(HWND hw, WPARAM wp,LPARAM lp){
	clearBitmap(0);//画像0をクリア
    InvalidateRect(hw,NULL,TRUE);
}
void procCreate(HWND hw, WPARAM wp,LPARAM lp){//初期処理
	createBitmap(hw);MoveWindow(hw,0,0,850,460,TRUE);
	hBM[1]=LoadBitmap(hInstance, TEXT("suzume"));//画像1に表示
	GetObject(hBM[1],sizeof(BITMAP), &bM[1]);
    InvalidateRect(hw,NULL,TRUE);
}
void procPaint(HWND hw, WPARAM wp,LPARAM lp){//画面表示
	PAINTSTRUCT ps;	
	HDC hdc  = BeginPaint(hw,&ps);
	SelectObject(hBuff[0],hBM[0]);SelectObject(hBuff[1],hBM[1]);
	BitBlt(hdc,10,YST,GWIDTH, GHIGHT, hBuff[0], 0,0,SRCCOPY);
	BitBlt(hdc,XST,YST,GWIDTH, GHIGHT, hBuff[1], 0,0,SRCCOPY);
	EndPaint(hw,&ps);
}
LRESULT CALLBACK WndProc(HWND hw, UINT msg, WPARAM wp,LPARAM lp){
	switch(msg){
 	case WM_DESTROY : PostQuitMessage(0)  ; return 0;
	case WM_CREATE  : procCreate(hw,wp,lp); return 0;
	case WM_LBUTTONDOWN  : procLButtonDown(hw,wp,lp); return 0;
	case WM_RBUTTONDOWN  : procRButtonDown(hw,wp,lp); return 0;
	case WM_PAINT   : procPaint (hw,wp,lp); return 0;
	}
	return DefWindowProc(hw,msg,wp,lp);
}


【リスト6-3】 R/G/B別の配列の用意とデータの設定
Complex beforR[GWDFT+1][GHDFT+1],afterR[GWDFT+1][GHDFT+1];
Complex reverceR[GWDFT+1][GHDFT+1];
Complex beforG[GWDFT+1][GHDFT+1],afterG[GWDFT+1][GHDFT+1];
Complex reverceG[GWDFT+1][GHDFT+1];
Complex beforB[GWDFT+1][GHDFT+1],afterB[GWDFT+1][GHDFT+1];
Complex reverceB[GWDFT+1][GHDFT+1];
void setComplex2DFT(){
	for(int i=0;i<GWDFT;i++){
		for(int j=0;j<GHDFT;j++){
			COLORREF C=pixel(i+80,j+80);
			beforR[i][j]=complex((double)RED(C)/255.0, 0);
			beforG[i][j]=complex((double)GREEN(C)/255.0, 0);
			beforB[i][j]=complex((double)BLUE(C)/255.0, 0);
		}
		beforR[i][GHDFT]=beforR[i][0];
		beforG[i][GHDFT]=beforG[i][0];
		beforB[i][GHDFT]=beforB[i][0];
	}
	beforR[GWDFT][GHDFT]=beforR[0][0];
	beforG[GWDFT][GHDFT]=beforG[0][0];
	beforB[GWDFT][GHDFT]=beforB[0][0];
}
  
【リスト6-4】 R/G/B別のフーリエ変換と逆変換
void procDFT(){
	setComplex2DFT();
	dspValue(beforR, 0,  0, 1);
	dspValue(beforG, 0,100, 2);
	dspValue(beforB, 0,200, 3);
	dspRGB(beforR,beforG,beforB, 0,300);

	DFT(beforR,afterR); DFT(beforG,afterG); DFT(beforB,afterB);
	dspSpec(afterR, 100,   0, 0 );
	dspSpec(afterR, 200,   0, 1 );
	dspSpec(afterG, 100, 100, 0 );
	dspSpec(afterG, 200, 100, 1 );
	dspSpec(afterB, 100, 200, 0 );
	dspSpec(afterB, 200, 200, 1 );

	IDFT(afterR,reverceR); IDFT(afterG,reverceG);	IDFT(afterB,reverceB);
	dspValue(reverceR, 300,  0, 1);
	dspValue(reverceG, 300,100, 2);
	dspValue(reverceB, 300,200, 3);
	dspRGB(reverceR,reverceG,reverceB, 300,300);
}


【リスト6-4】 単純なローパスフィルタ
void LPF(int r){// (方法１) 円の外側をカット
	int rr=r*r, DX, DY;Complex C0=complex(0,0);
	int X0=GWDFT/2, Y0=GHDFT/2;
	for(int i=0;i<GWDFT; i++)for(int j=0; j<GHDFT;j++){
		DX=i-X0;DY=j-Y0;if((DX*DX+DY*DY)>rr) afterDFT[i][j]=C0;	
	}
}
void LPF2(int r){// (方法2) 四隅から半径R分をカット
	int rr=r*r;
	int MX=GWDFT-1, MY=GHDFT-1;Complex C0=complex(0,0);
	for(int i=0;i<r; i++) for(int j=0; j<r;j++){
		if((i*i+j*j)<rr)
			afterDFT[i][j]= afterDFT[MX-i][j]
				= afterDFT[i][MY-j] = afterDFT[MX-i][MY-j] = C0;
	}
}

【リスト6-5】 中心からの距離によって強調度を変化させるローパスフィルタ
void LPFX(int r){
	double R=r, DX, DY;
	double X0=GWDFT/2, Y0=GHDFT/2,XR, FF, alfa=0.1;
	for(int i=0;i<GWDFT; i++)	for(int j=0; j<GHDFT;j++){
		DX=i-X0;DY=j-Y0;XR=sqrt(DX*DX+DY*DY);
		if(XR>R) FF=0.5*exp(alfa*(R-XR)); else FF=1-0.5*exp(alfa*(XR-R));
		afterDFT[i][j]=multR(afterDFT[i][j], FF);	
	}
}
  
【リスト6-6】 単純なハイパスフィルタ
void HPF(int r){
	int rr=r*r, DX, DY;Complex C0=complex(0,0);
	int X0=GWDFT/2, Y0=GHDFT/2;
	for(int i=0;i<GWDFT; i++)for(int j=0; j<GHDFT;j++){
		DX=i-X0;DY=j-Y0;if((DX*DX+DY*DY)<=rr)afterDFT[i][j]=C0;	
	}
}

【リスト6-7】 中心からの距離によるハイパスフィルタ
void HPFX(int r){
	double R=r, DX, DY;
	double X0=GWDFT/2, Y0=GHDFT/2,XR, FF, alfa=0.1;
	for(int i=0;i<GWDFT; i++)	for(int j=0; j<GHDFT;j++){
		DX=i-X0;DY=j-Y0;XR=sqrt(DX*DX+DY*DY);
		if(XR>R) FF=1-0.5*exp(alfa*(R-XR)); else FF=0.5*exp(alfa*(XR-R));
		afterDFT[i][j]=multR(afterDFT[i][j],FF);	
	}
}

【リスト6-8】  1次元FFT
Complex w_tbl[GWDFT/2];	//回転因子計算用テーブル
int b_tbl[GWDFT]; 		//ビット逆順交換用テーブル
static int n_old=0;
int fft_tbl(Complex * x, int n_FFT){
	Complex xtmp; int pow2, j,jnh,jp,jx,num, ne, n_half, n_half2;
	num=labs(n_FFT);
	if(num !=n_old){   //1回目の呼び出しのときのみテーブル作成
		pow2=1;			//要素数が2のべき乗かどうかを判定
		for(j=0;j<32;j++) if(num==pow2) break; else pow2<<=1;
		if(j>=32){
			MessageBox(NULL,TEXT("\n データ数が2のべき乗になっていません"),
							TEXT("Error"), MB_OK);
			return 1;
		}
		double argV=-2.0*PI/num; //回転因子計算用テーブル作成
		for(j=0;j<num/2;j++) w_tbl[j]=expJ(argV*j);
		n_half=num/2; b_tbl[0]=0; //ビット逆順交換用テーブル作成
		for(ne=1; ne<num;ne<<=1){
			for(jp=0;jp<ne;jp++)b_tbl[jp+ne]=b_tbl[jp]+n_half;
			n_half>>=1;
		}
		n_old=num;
	}
	n_half=num/2;//FFT計算開始
	for(ne=1; ne<num; ne<<=1){
		n_half2=n_half*2;
		for(jp=0; jp<num; jp += n_half2){
			jx=0;
			for(j=jp; j<jp+n_half; j++){
				jnh=j + n_half; xtmp = x[j]; x[j] = addC(xtmp,x[jnh]);
				x[jnh]=multC(subC(xtmp,x[jnh]), w_tbl[jx]);
				jx += ne;
			}
		}
		n_half>>=1;
	}
	//ビット逆順入れ換え
	for(j=0; j<num;j++) if(j<b_tbl[j]) swap(Complex, x[j],x[b_tbl[j]]);
	if(n_FFT<0) for(j=0;j<num; j++)	x[j]=divR(x[j],num);//逆FFT
	return 0;
}

【リスト6-9】  2次元FFT
void fft_2d(Complex beforDFT[][GHDFT+1],Complex afterDFT[][GHDFT+1],
					int n_fft){
	int n1, n2,num=labs(n_fft);  Complex u[GWDFT+1];
	for(n1=0; n1<num; n1++){
		for(n2=0; n2<num; n2++){
			double DD=((n1+n2) & 1) ? -1.0: 1.0;
			afterDFT[n1][n2]=multR(beforDFT[n1][n2],DD);
		}
		if(fft_tbl(afterDFT[n1], n_fft))return;
	} 
	for(n2=0; n2<num; n2++){
		for(n1=0; n1<num; n1++) u[n1]=afterDFT[n1][n2];
		fft_tbl(u, n_fft);
		for(n1=0; n1<num; n1++) afterDFT[n1][n2]=u[n1];
	}
} 
void procDFT(){//以下はテスト用プログラムです。
	setComplexDFT(); dspValue(beforDFT, 0,0, 0);//データ設定
	fft_2d(beforDFT,afterDFT, GWDFT);           //FFT呼び出し
	dspSpec(afterDFT, 100, 0, 0 ); dspSpec(afterDFT, 200, 0, 1 );
	IDFT(afterDFT,reverceDFT);      //通常の逆フーリエ変換で結果を確認
	dspValue(reverceDFT, 300,0,0);
}

【リスト6-10】  高速化2次元FFT
void fft_2d(Complex beforDFT[][GHDFT+1],Complex afterDFT[][GHDFT+1],
						 int n_fft){
	int n1, n2,num=labs(n_fft),num1=num+1; Complex u[GWDFT+1];
	memmove(afterDFT, beforDFT,sizeof(Complex)*num1*num1);
	for(n1=0; n1<num; n1++) if(fft_tbl(afterDFT[n1], n_fft))return;
	for(n2=0; n2<num; n2++){
		for(n1=0; n1<num; n1++) u[n1]=afterDFT[n1][n2];
		fft_tbl(u, n_fft);
		for(n1=0; n1<num; n1++) afterDFT[n1][n2]=u[n1];
	}
}

【リスト6-11】  スペクトル位置の移動(GWDFT=64, GHDFT=64)
void setDspSpec(Complex TM[][GHDFT+1], Complex DT[][GHDFT+1]){
	int w2=GWDFT/2, h2=GHDFT/2;
	for(int x=0; x<GWDFT; x++)for(int y=0;y<GHDFT;y++){
		int x2=x+w2, y2=y+h2;
		if(x2>GWDFT) x2= x2-GWDFT; if(y2>GHDFT) y2= y2-GHDFT;
		TM[x2][y2]=DT[x][y];
	}
} 
Complex temp[GWDFT+1][GHDFT+1];
void procDFT(){// テスト用
	setComplexDFT(); dspValue(beforDFT, 0,0, 0);
	fft_2d(beforDFT,afterDFT, GWDFT); setDspSpec(temp, afterDFT);
	memmove(afterDFT, temp,sizeof(Complex)*(GWDFT+1)*(GHDFT+1));
	dspSpec(afterDFT, 100, 0, 0 ); dspSpec(afterDFT   , 200, 0, 1 );
	IDFT(afterDFT,reverceDFT)    ; dspValue(reverceDFT, 300,0,0);
}

【リスト6-12】  2次元高速逆フーリエ変換(GWDFT=64, GHDFT=64)
void rev_fft_2d(Complex afterDFT[][GHDFT+1],Complex reverceDFT[][GHDFT+1], int n_fft){
	int num=labs(n_fft);
	fft_2d(afterDFT,reverceDFT, -num);

	for(int k=0;k<=num;k++)for(int N=0, N2=num-1;N<num/2;N++,N2--)
			swap(Complex, reverceDFT[k][N],reverceDFT[k][N2]);
	for(int k=0;k<=num;k++)for(int N=0, N2=num-1;N<GWDFT/2;N++,N2--)
				swap(Complex, reverceDFT[N][k],reverceDFT[N2][k]);
	}

void procDFT(){
	setComplexDFT(); dspValue(beforDFT, 0,0, 0);
	fft_2d(beforDFT,afterDFT, GWDFT); setDspSpec(temp, afterDFT);
	dspSpec(afterDFT, 100, 0, 0 ); dspSpec(afterDFT, 200, 0, 1 );
	rev_fft_2d(afterDFT,reverceDFT, GWDFT);
	dspValue(reverceDFT, 300,0,0);
}

【リスト6-13】  2次元高速逆フーリエ変換(GWDFT=64, GHDFT=64)
Complex temp[GWDFT+1][GHDFT+1];
void procDFT(){
	setComplexDFT(); dspValue(beforDFT, 0,0, 0);
	fft_2d(beforDFT,afterDFT, GWDFT);
	setDspSpec(temp, afterDFT); 
	dspSpec(temp, 100, 0, 0 ); dspSpec(temp, 200, 0, 1 );  
	rev_fft_2d(afterDFT,reverceDFT, GWDFT);
	dspValue(reverceDFT, 300,0,0);
}

【リスト6-14】  2次元DCT
double red0[GWDFT][GHDFT],green0[GWDFT][GHDFT],blue0[GWDFT][GHDFT];
double red1[GWDFT][GHDFT],green1[GWDFT][GHDFT],blue1[GWDFT][GHDFT];
double red2[GWDFT][GHDFT],green2[GWDFT][GHDFT],blue2[GWDFT][GHDFT];
double trt[GWDFT][GHDFT],wrt[GWDFT][GWDFT],hrt[GHDFT][GHDFT]; //作業用配列
oid setCosDT(){
	for(int i=0;i<GWDFT;i++)for(int j=0;j<GHDFT;j++){
		COLORREF C=pixel(i+80,j+80);
		red0  [i][j]=(double)RED  (C)/255.0;
		green0[i][j]=(double)GREEN(C)/255.0;
		blue0 [i][j]=(double)BLUE (C)/255.0;
	}
}	
void DCT(double dt[][GHDFT],double out[][GHDFT],int inv ){
	double nm1, nm2; int stID;
	if(inv){ nm1 =  0.5; nm2 = 2.0/GWDFT; stID   = 1;}
	else   { nm1 =  0.0; nm2 = 1.0      ; stID   = 0;}
	double M2    = PI/(double)GWDFT;
	for(int i=0; i<GWDFT; i++){// x方向変換のためのテーブル設定
		for(int j=stID; j<GWDFT; j++){
			if(inv) wrt[i][j]=cos(M2*(i+0.5)*j);
			else    wrt[i][j]=cos(M2*(j+0.5)*i); 
		}
	}
	for(int y=0;y<GHDFT; y++)for(int x=0;x<GWDFT; x++){// x方向DCT
		trt[x][y]=0.0;
		for(int xx=stID; xx<GWDFT; xx++) trt[x][y]+=dt[xx][y]*wrt[x][xx];
		trt[x][y] += dt[0][y]*nm1; trt[x][y] *= nm2;
	}
	M2 = PI/(double)GHDFT; nm2=(inv)?2.0/(double)GHDFT:1;
	for(int i=0; i<GHDFT; i++){// y方向変換のためのテーブル設定
		for(int j=stID; j<GHDFT; j++){
			if(inv) hrt[i][j]=cos(M2*(i+0.5)*j);
			else    hrt[i][j]=cos(M2*(j+0.5)*i); 
		}
	}
	for(int x=0;x<GWDFT; x++)for(int y=0;y<GWDFT; y++){// y方向DCT
		out[x][y]=0.0;
		for(int yy=stID; yy<GHDFT; yy++)	
		    out[x][y]+=trt[x][yy]*hrt[y][yy];
		out[x][y] += trt[x][0]*nm1;	out[x][y] *= nm2;
	}
}
void dspDBL(double DT[][GHDFT],int xs, int ys, int flag){
	for(int i=0;i<GWDFT;i++)for(int j=0;j<GHDFT;j++){
		int C=(int)(DT[i][j]*255+0.5);
		switch(flag){
		case 1: C=RGB(C, 0, 0);break;
		case 2: C=RGB(0, C, 0);break;
		case 3: C=RGB(0, 0, C);break;
		default:C=RGB(C,C,C);break;
		}
		setPixel(i+xs,j+ys, C);
	}
}
void cosTransform(){
	setCosDT();
	DCT(red0, red1,0);	DCT(green0, green1,0); DCT(blue0, blue1,0);
	DCT(red1, red2,true);DCT(green1, green2,true); DCT(blue1, blue2,true);
	dspDBL(red0,100,0,1); dspDBL(green0,200,0, 2); dspDBL(blue0,300,0,3);
	dspDBL(red1,100,100,1);dspDBL(green1,200,100,2);dspDBL(blue1,300,100,3);
	dspDBL(red2,100,200,1);dspDBL(green2,200,200,2);dspDBL(blue2,300,200,3);
	for(int i=0;i<GWDFT;i++)for(int j=0;j<GHDFT;j++){
		setPixel(i,j, pixel(i+80,j+80));
		setPixel(i,j+200, RGB((int)(red2  [i][j]*255+0.5),
			    (int)(green2[i][j]*255+0.5), (int)(blue2 [i][j]*255+0.5)));
	}
}