二次元離散フーリエ変換、逆変換を実装したい

Question

二次元離散フーリエ変換および逆変換のプログラムを実装したいです

二次元離散フーリエ変換および逆変換のプログラムを書いたのですが、変換、逆変換後に結果が戻りません。

入力信号は画像ファイルから１画素あたり８ビットの輝度値が格納された配列とします。

該当のソースコード

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

#define N1 128
#define N2 128
#define H 54 //54

int main(int argc, char **argv)
{
	FILE *fp;
	double x[N1][N2],ReH[N1][N2],ImH[N1][N2],ReX[N1][N2],ImX[N1][N2];
	int k1,k2,n1,n2;//Re~実部　Im~虚部　

	
	unsigned char header[H];
	unsigned char img[N1][N2][3];
 	unsigned char in[N1][N2];

	
	int i,j,k;
	int height_original=N1;	//画像の縦のサイズ 
	int width_original =N2;	//画像の横のサイズ
 	
	
	
	if( argc != 2 ){
        fprintf( stderr, "Usage: ./a input.bmp \n" );
        exit(-1);
    }
    if( (fp = fopen(argv[1], "rb")) == NULL ){
        fprintf( stderr, "Cannot open %s\n", argv[1] );
        exit(-2);
    }
	fread(header,1,H,fp); /* ヘッダ(54バイト)を読み込み */
	fread(img,1, height_original * width_original * 3,fp); /* 残りはデータ(最下行から順に入る) */
	fclose(fp);
	
	fp=fopen("af6x.txt","w");
 	for (i=0;i<height_original;i++){//縦の長さ
 		for (j=0;j<width_original;j++){//横の長さ
 			in[i][j]  = ( img[i][j][0] + img[i][j][1] + img[i][j][2])/3;
 			x[i][j]=(double)in[i][j];
 			fprintf(fp,"%5.5f\n",x[i][j]);
 		}
 		fprintf(fp,"\n\n");
 	}
 	fclose(fp);


	
	//フーリエ変換　行と列交互に変換

	fp=fopen("af6fourier.txt","w");
	for (n2=0;n2<N2;n2++){
		for (k1=0;k1<N1;k1++){
			ReH[k1][n2]=ImH[k1][n2]=0.0;
			for (n1=0;n1<N1;n1++){
				ReH[k1][n2]+=x[n1][n2]*cos((n1*2*M_PI*k1)/(double)N1);
				ImH[k1][n2]+=-x[n1][n2]*sin((n1*2*M_PI*k1)/(double)N1);
			}
		}
	}
	
	for (k1=0;k1<N1;k1++){
		for (k2=0;k2<N2;k2++){
			ReX[k1][k2]=ImX[k1][k2]=0.0;
			for (n2=0;n2<N2;n2++){
				ReX[k1][k2]+=ReH[k1][n2]*cos((n2*2*M_PI*k2)/(double)N2)+ImH[k1][n2]*sin((n2*2*M_PI*k2)/(double)N2);
				ImX[k1][k2]+=-ReH[k1][n2]*sin((n2*2*M_PI*k2)/(double)N2)+ImH[k1][n2]*cos((n2*2*M_PI*k2)/(double)N2);
			}
			fprintf(fp,"%5.5f\t%5.5f\n",ReX[k1][k2],ImX[k1][k2]);
		}
		fprintf(fp,"\n");
		
	}
	fclose(fp);
	
	
	
	//逆フーリエ変換
	double Rex[N1][N2],Imx[N1][N2],Reh[N1][N2],Imh[N1][N2];

	fp=fopen("af6gyakufourier.txt","w");
	for (k2=0;k2<N2;k2++){
		for (n1=0;n1<N1;n1++){
			Reh[n1][k2]=Imh[n1][k2]=0.0;
			for (k1=0;k1<N1;k1++){
				Reh[n1][k2]+=(ReX[k1][k2]*cos(n1*2*M_PI*k1/(double)N1)-ImX[k1][k2]*sin(n1*2*M_PI*k1/(double)N1))/(double)N1;
				Imh[n1][k2]+=(ReX[k1][k2]*sin(n1*2*M_PI*k1/(double)N1)+ImX[k1][k2]*cos(n1*2*M_PI*k1/(double)N1))/(double)N1;		
			}
			
		}
	}
	
	
	for (n1=0;n1<N1;n1++){
		for (n2=0;n2<N2;n2++){
			Rex[n1][n2]=Imx[n1][n2]=0.0;
			for (k2=0;k2<N2;k2++){
				Rex[n1][n2]+=(Reh[n1][k1]*cos(n2*2*M_PI*k2/(double)N2)-Imh[n1][k1]*sin(n2*2*M_PI*k2/(double)N2))/(double)N2;
				Imx[n1][n2]+=(Reh[n1][k1]*sin(n2*2*M_PI*k2/(double)N2)+Imh[n1][k1]*cos(n2*2*M_PI*k2/(double)N2))/(double)N2;	

			}
			fprintf(fp,"%5.5f\t%5.5f\n",Rex[n1][n2],Imx[n1][n2]);
		}
		
	}

 	fclose(fp); 
	
	fclose(fp);


	return 0;
}

試したこと

定義通りにプログラムを組んだつもりですがうまくできません...

補足情報（FW/ツールのバージョンなど）

ここにより詳細な情報を記載してください。

Comments

[y_waiwai]

なにがどういうふうにうまくいかないんでしょうか

[VanS]

入力信号をフーリエ変換した後に逆フーリエ変換をして元に戻ることを確認したいです。説明不足ですみません。。。

[jbpb0]

> 結果が戻りません。

結果を何かで割ったら元に戻りませんか？

[jbpb0]

> 変換、逆変換後に

いきなりそこまで行かず、変換だけの結果が、既存のコードでの結果と合うか、確認したらいかがでしょうか？

[VanS]

既存のコードでの結果、正しい変換後の値はどのようにしてもとめればいいのでしょうか。。？

[VanS]

xを求める際にあるiのとき0<=j<N2-1の範囲で配列の和を求めると逆変換後に同じ値が出てきます。。

[jbpb0]

もしPythonとかMATLABとかGnu Octaveとか、数値計算系でよく使われるツールがあるなら、それでFFTしたのと比べるとか

あるいは、正弦波(Sin、-1～1)を10周期くらい、1°間隔でデータを作って、それをFFTしてみるとか
その場合、データの開始と終了がきれいにつながる(終了の次が開始と同じになる)ように作れば、FFT後は二カ所にだけデータが現れます
たとえば、元データが10周期なら、FFT後は±10Hzにだけデータが現れるため、11番目と、後ろから10番目に結果が現れます
周波数が0Hzから始まるため、10Hzは11番目になります
また、後半は負の周波数で、一番最後が-1Hz、その前が-2Hz、と並びます

なお、元データの平均が0ではない場合は、FFT後は0Hzにもデータが現れます

[VanS]

なるほど、、やってみます

[jbpb0]

> データの開始と終了がきれいにつながる(終了の次が開始と同じになる)ように作れば

を補足します

たとえば、0°から始まり1°間隔でSinを1周期分のデータを作る場合、0, 1...358, 359°のSinを計算する、ということです
最後は360°ではなく、その一つ前の359°になります
Sin(0°)=Sin(360°)だから、最後をSin(359°)にすれば、最初と最後がきれいにつながります(最後の次に最初が来ても大丈夫)

1周期じゃなくても、何周期でも同じです
最後は、最初と同じものの一つ前

Answer 1

最後のforループ内における変数が間違っていました・・・
Rex[n1][n2]+=(Reh[n1][k1]cos(n22M_PIk2/(double)N2)-Imh[n1][k1]sin(n22M_PIk2/(double)N2))/(double)N2;
Imx[n1][n2]+=(Reh[n1][k1]sin(n22M_PIk2/(double)N2)+Imh[n1][k1]cos(n22M_PIk2/(double)N2))/(double)N2;

正しくはk1→k2です
ご協力いただいた方々ありがとうございました。