質問をすることでしか得られない、回答やアドバイスがある。

15分調べてもわからないことは、質問しよう!

新規登録して質問してみよう
ただいま回答率
85.48%
C

C言語は、1972年にAT&Tベル研究所の、デニス・リッチーが主体となって作成したプログラミング言語です。 B言語の後継言語として開発されたことからC言語と命名。そのため、表記法などはB言語やALGOLに近いとされています。 Cの拡張版であるC++言語とともに、現在世界中でもっとも普及されているプログラミング言語です。

Q&A

解決済

1回答

3889閲覧

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

VanS

総合スコア6

C

C言語は、1972年にAT&Tベル研究所の、デニス・リッチーが主体となって作成したプログラミング言語です。 B言語の後継言語として開発されたことからC言語と命名。そのため、表記法などはB言語やALGOLに近いとされています。 Cの拡張版であるC++言語とともに、現在世界中でもっとも普及されているプログラミング言語です。

0グッド

1クリップ

投稿2021/01/17 07:02

編集2021/01/17 07:03

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

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

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

該当のソースコード

C

1#include <stdio.h> 2#include <math.h> 3#include <stdlib.h> 4 5#define N1 128 6#define N2 128 7#define H 54 //54 8 9int main(int argc, char **argv) 10{ 11 FILE *fp; 12 double x[N1][N2],ReH[N1][N2],ImH[N1][N2],ReX[N1][N2],ImX[N1][N2]; 13 int k1,k2,n1,n2;//Re~実部 Im~虚部  14 15 16 unsigned char header[H]; 17 unsigned char img[N1][N2][3]; 18 unsigned char in[N1][N2]; 19 20 21 int i,j,k; 22 int height_original=N1; //画像の縦のサイズ 23 int width_original =N2; //画像の横のサイズ 24 25 26 27 if( argc != 2 ){ 28 fprintf( stderr, "Usage: ./a input.bmp \n" ); 29 exit(-1); 30 } 31 if( (fp = fopen(argv[1], "rb")) == NULL ){ 32 fprintf( stderr, "Cannot open %s\n", argv[1] ); 33 exit(-2); 34 } 35 fread(header,1,H,fp); /* ヘッダ(54バイト)を読み込み */ 36 fread(img,1, height_original * width_original * 3,fp); /* 残りはデータ(最下行から順に入る) */ 37 fclose(fp); 38 39 fp=fopen("af6x.txt","w"); 40 for (i=0;i<height_original;i++){//縦の長さ 41 for (j=0;j<width_original;j++){//横の長さ 42 in[i][j] = ( img[i][j][0] + img[i][j][1] + img[i][j][2])/3; 43 x[i][j]=(double)in[i][j]; 44 fprintf(fp,"%5.5f\n",x[i][j]); 45 } 46 fprintf(fp,"\n\n"); 47 } 48 fclose(fp); 49 50 51 52 //フーリエ変換 行と列交互に変換 53 54 fp=fopen("af6fourier.txt","w"); 55 for (n2=0;n2<N2;n2++){ 56 for (k1=0;k1<N1;k1++){ 57 ReH[k1][n2]=ImH[k1][n2]=0.0; 58 for (n1=0;n1<N1;n1++){ 59 ReH[k1][n2]+=x[n1][n2]*cos((n1*2*M_PI*k1)/(double)N1); 60 ImH[k1][n2]+=-x[n1][n2]*sin((n1*2*M_PI*k1)/(double)N1); 61 } 62 } 63 } 64 65 for (k1=0;k1<N1;k1++){ 66 for (k2=0;k2<N2;k2++){ 67 ReX[k1][k2]=ImX[k1][k2]=0.0; 68 for (n2=0;n2<N2;n2++){ 69 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); 70 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); 71 } 72 fprintf(fp,"%5.5f\t%5.5f\n",ReX[k1][k2],ImX[k1][k2]); 73 } 74 fprintf(fp,"\n"); 75 76 } 77 fclose(fp); 78 79 80 81 //逆フーリエ変換 82 double Rex[N1][N2],Imx[N1][N2],Reh[N1][N2],Imh[N1][N2]; 83 84 fp=fopen("af6gyakufourier.txt","w"); 85 for (k2=0;k2<N2;k2++){ 86 for (n1=0;n1<N1;n1++){ 87 Reh[n1][k2]=Imh[n1][k2]=0.0; 88 for (k1=0;k1<N1;k1++){ 89 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; 90 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; 91 } 92 93 } 94 } 95 96 97 for (n1=0;n1<N1;n1++){ 98 for (n2=0;n2<N2;n2++){ 99 Rex[n1][n2]=Imx[n1][n2]=0.0; 100 for (k2=0;k2<N2;k2++){ 101 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; 102 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; 103 104 } 105 fprintf(fp,"%5.5f\t%5.5f\n",Rex[n1][n2],Imx[n1][n2]); 106 } 107 108 } 109 110 fclose(fp); 111 112 fclose(fp); 113 114 115 return 0; 116}

試したこと

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

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

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

気になる質問をクリップする

クリップした質問は、後からいつでもMYページで確認できます。

またクリップした質問に回答があった際、通知やメールを受け取ることができます。

バッドをするには、ログインかつ

こちらの条件を満たす必要があります。

y_waiwai

2021/01/17 07:12

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

2021/01/17 07:22

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

2021/01/17 09:33

> 結果が戻りません。 結果を何かで割ったら元に戻りませんか?
jbpb0

2021/01/17 09:35

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

2021/01/17 13:25

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

2021/01/17 13:27

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

2021/01/17 13:56 編集

もし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

2021/01/18 01:25

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

2021/01/18 01:34

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

回答1

0

自己解決

最後の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です
ご協力いただいた方々ありがとうございました。

投稿2021/01/20 15:00

VanS

総合スコア6

バッドをするには、ログインかつ

こちらの条件を満たす必要があります。

あなたの回答

tips

太字

斜体

打ち消し線

見出し

引用テキストの挿入

コードの挿入

リンクの挿入

リストの挿入

番号リストの挿入

表の挿入

水平線の挿入

プレビュー

15分調べてもわからないことは
teratailで質問しよう!

ただいまの回答率
85.48%

質問をまとめることで
思考を整理して素早く解決

テンプレート機能で
簡単に質問をまとめる

質問する

関連した質問