C言語プログラムの離散フーリエ変換について教えてください。「C言語による画像再構成の基礎」という本のプログラムをもとに二次元画像をDFT(通常の離散フーリエ変換)→InveresFFT(逆高速フーリエ変換)すると画像が左右反転、上下反転してしまいます。DFT→InverseDFTやFFT→InverseFFTだとそのようにはなりません。通常のDFTとFFTのアルゴリズムの違いからしかたがないのでしょうか?それともプログラムの変更で修正できるのでしょうか?どうしてもDFT→InverseFFTでがぞうをもとに戻したいのです。
サンプルページ http://www.iryokagaku.co.jp/frame/03-honwosagasu …
P4-14fourier2d1d.c (離散フーリエ変換DFT)   P4-15fft.c(高速フーリエ変換)プログラムです

このQ&Aに関連する最新のQ&A

A 回答 (1件)

なぜ「どうしてもDFT→InverseFFTでがぞうをもとに戻したい」んだろう.



言い換えれば「なぜ FFT→IFFT じゃだめ」なんだろう.
    • good
    • 0

このQ&Aに関連する人気のQ&A

お探しのQ&Aが見つからない時は、教えて!gooで質問しましょう!

このQ&Aを見た人はこんなQ&Aも見ています

このQ&Aを見た人が検索しているワード

このQ&Aと関連する良く見られている質問

Qc言語でDFTのプログラムを作成したのですが

c言語でDFTのプログラムを作成しました。
以下にソースを載せます。
#include<stdio.h>
#include<math.h>
#include<stdlib.h>
#include<time.h>
#define PI 3.141592653589793
#define N 64 //データ数


DFT(double result[]){
int i,k;
double A[N],B[N],T=0; //A[N]:実数部,B[N]:虚数部
double a,b;
for(k=0;k<N;++k){
a=b=0;
for(i=0;i<N;++i){
a=a+result[i]*cos(2*PI*i*k/N);     
b=b+(-1.0)*result[i]*sin(2*PI*i*k/N);
}
A[k]=a/N;
B[k]=b/N;
}

for(i=0;i<N;++i){
printf("[%f秒]:Re:%f,Im:%f\n",T,A[i],B[i]); //変換後の値を表示
T=T+(0.1/N);
}

}

main(){
int i;
double T=0;
double Sampdata;
double result[N];

for(i=0;i<N;++i){
Sampdata=5*sin(20*PI*T);      //0~0.1秒間をN個にサンプリング
result[i]=Sampdata; //サンプリングデータを代入
T=T+(0.1/N);
}
clock_t start,end; //処理時間計測開始
start=clock();
DFT(result);
end=clock();
printf("%.2f秒かかりました\n",(double)(end-start)/CLOCKS_PER_SEC); //処理時間表示
}

元信号には5sin(20πt)の値を入れています。この信号は周期は0.1secです。
これでフーリエ変換を行うとデータ数N/2を中心に対称なデータが出てくるのですが、処理が終わるのが早い気がするんです。
例えば2^15個のデータで実行しても2分もかからずに処理が終わってしまいます。一応、対称性が出てるとはいえ、終わるのが早すぎる気がするのですが、おかしい所があれば教えていただけると嬉しいです。
よろしくお願いします。

c言語でDFTのプログラムを作成しました。
以下にソースを載せます。
#include<stdio.h>
#include<math.h>
#include<stdlib.h>
#include<time.h>
#define PI 3.141592653589793
#define N 64 //データ数


DFT(double result[]){
int i,k;
double A[N],B[N],T=0; //A[N]:実数部,B[N]:虚数部
double a,b;
for(k=0;k<N;++k){
a=b=0;
for(i=0;i<N;++i){
a=a+result[i]*cos(2*PI*i*k/N);     
b=b+(-1.0)*result[i]*sin(2*PI*i*k/N);
}
A[k]=a/N;
B[k]...続きを読む

Aベストアンサー

計算に使っているコンピュータの性能はどれくらいなのでしょう?

このプログラムで一番時間がかかると思われるのは
> for(i=0;i<N;++i){
> printf("[%f秒]:Re:%f,Im:%f\n",T,A[i],B[i]); //変換後の値を表示>
> T=T+(0.1/N);
> }

この出力の部分です。

手許のcore i7 2.66GHz
gcc -O0 (最適化無し)
DFT(result); を2^15回ループ
結果出力部分をコメントアウト
で22秒程でした。

CPUの速度、コンパイラの最適化等でもっと速くなります。
2分かからないくらいなら、別に変では無いと思います。

心配なら、この計算結果を、別の方法で計算したもの(実績のあるDFT計算ライブラリ、Excel等の別ソフト等)と比較してみては?

Q離散フーリエ変換をC言語でどの様に書けばいいですか?

C言語でDFT離散フーリエ変換を書くにはどの様に書けばよろしいですか?

Googleで検索すれば書き方は出てくるのですが、使ってる関数がいまいちよく分かりません。


・データの入力

・フーリエ変換の計算

・結果の出力というのをやればいいのは理解できるのですが、C言語でどの様に書けばいいか分からなくて…

Aベストアンサー

学習中である場合とする解答です。使用パソコンは Linuxまたは Mac OSX などの UNIX系OSです。

<・データの入力>
1)エディタで専用データファイルを作ることから始めます。
 マトリックスは、それに対応した数値をファイル内配置します。

  N
  W0 W0 W0 W0 ....
  W0 W1 W2 W3 ....
  W0 W2 W4 W6 ....
  ....
  N0
  N1
  N2
   .
   .

 上記定義されたマトリックスを ↓のようにファイルに書き込む
 仮に整数値としていますが、実数ならば実数値を書き込んで下さい。

5
  0 0 0 0 ....
0 1 2 3 ....
0 2 4 6 ....
  ....
  0
  1
  2
  .
  .


 各数値間は半角スペースで区切ってそれぞれの数式定義対応マトリックスを作ります。
 ファイル名は半角英字がエラーなく行えるので英字ファイル名を使うことを勧めます。

2)データの読み込み
 データの入力は scanf() を使います。
 最初にNを読み込めば、マトリックス行数がわかるため for() を使ってプログラムしますが、慣れない場合は腕力で scanf() 関数を連発するのもひとつの方法です。
 データファイル読み込みに際しての疑問は、作成したプログラムを起動する際、ターミナルのシェルプロンプトとから
   ./a.out<データファイル名
 と打鍵すればデータファイルを取り込むことができます。
 以上をCで表すと次のようになります。scanf() はそちらで勉強して下さい。
 プログラム作成に際しては、一気に書き込まず。途中で printf() を入れてデータ取り込みを確認されることを勧めます。

3)Cプログラム
 /* magatai.c DFT program
  * file name: magatai.c
  * compile: gcc magatai.c
  * execution: ./a.out<data_file
  */

 #define N 10 取り込みデータ数+α

 double n;
 double w[N][N];
 double x[N];

 int main(void)
 {
  int i, j;
  scanf(%d, &n);
  for (i = 0; i < n; i++)
   scanf(%......
   .... ↓のフーリエ式が入る。
   .... ↓↓の結果の出力が入る。
  return 0;
 }


<・フーリエ変換>
 これは for() の入れ子になります。
  for (i = 0; i < n; i++) {
   for (j = 0; j < n; j++ {
    フーリエ変換式;
   }
  }


<・結果の出力>
 printf() を使います。↑のフーリエ・プログラムに続いて書きます。
  for (i = 0; i < n; i++ )
   printf("%f?n", x[i]);

 計算結果を特定のファイルに残したいという場合は、プログラム起動の際、リダイレクトを使います。
   ./a.out<データファイル名>書き込むファイル名
 後は、 cat などを使ってファイルをリストすれば良いでしょう。

参考URL:http://laputa.cs.shinshu-u.ac.jp/~yizawa/InfSys1/basic/chap6/index.htm

学習中である場合とする解答です。使用パソコンは Linuxまたは Mac OSX などの UNIX系OSです。

<・データの入力>
1)エディタで専用データファイルを作ることから始めます。
 マトリックスは、それに対応した数値をファイル内配置します。

  N
  W0 W0 W0 W0 ....
  W0 W1 W2 W3 ....
  W0 W2 W4 W6 ....
  ....
  N0
  N1
  N2
   .
   .

 上記定義されたマトリックスを ↓のようにファイルに書き込む
 仮に整数値としていますが、実数ならば実数値を書き込んで下...続きを読む

Q離散フーリエ変換のプログラムについて

今DFTのプログラムをC言語で書いているんですが、うまく動いてくれません。
DFTの式はX(k)=1/N{Σx(k)*e^(-j2πkn/N)}のシグマの中をn=0からN-1まで足し合わせればいいと思っているのですがちがいますか。
下にプログラムを書きますのでお願いします。

void DFTcore(double x[],int N,double A[],double fai[],double a_rl[],double a_im[]){

    x[]はデータ
    Nはデータ数
    A[]は振幅
    fai[]は位相差
    a_rl[]は実数部、a_im[]は虚数部

double w,a;
int k,n,t;

for(k=0;k<N;k++){  //初期化
   a_rl[k]=0.0;  //実数部
   a_im[k]=0.0;  //虚数部
  }
    時間間隔は1秒
for(k=0;k<N;k++){
      for(n=0;n<N;n++){
  w=((2*PI)/N)*k*n;
        a_rl[k]=a_rl[k]+x[k]*cos(w);
  a_im[k]=-a_im[k]+x[k]*sin(w);
 }
  a_rl[k]=a_rl[k]/(double)N;実数部
  a_im[k]=-a_im[k]/(double)N; 虚数

  A[k]=sqrt(a_rl[k]*a_rl[k]+a_im[k]*a_im[k]);//振幅
   fai[k]=atan(a_im[k]/a_rl[k]);//位相

}
}

今DFTのプログラムをC言語で書いているんですが、うまく動いてくれません。
DFTの式はX(k)=1/N{Σx(k)*e^(-j2πkn/N)}のシグマの中をn=0からN-1まで足し合わせればいいと思っているのですがちがいますか。
下にプログラムを書きますのでお願いします。

void DFTcore(double x[],int N,double A[],double fai[],double a_rl[],double a_im[]){

    x[]はデータ
    Nはデータ数
    A[]は振幅
    fai[]は位相差
    a_rl[]は実数部、a_im[]は虚数部

double w,a;
...続きを読む

Aベストアンサー

正しくは、X(k)=1/N{Σx(n)*e^(-j2πkn/N)}
(ただし、k=f*N/f_sample , fは調べたい周波数、f_sampleはサンプリング周波数)
です。
とりあえずは、
a_rl[k]=a_rl[k]+x[n]*cos(w);
a_im[k]=a_im[k]+x[n]*sin(w);
とすれば動作確認はできます。
あとはご自分のプログラムの中でのkの意味(f*N/f_sampleは整数でなくてもよい)、f<(f_sample/2)、に注意し、プログラムを書きなおしてみてください。またきれいなスペクトルを得るには窓関数処理をお忘れなく。

数学、数式が大の苦手の私ですが、DFTはおぼろげに理解できたので、以前N88BASICという言語で(今は知っている人が少ないらしい)、DFTのプログラムを書いたことがあります。なにせクロック周波数12MHzのパソコンだったので、サインテーブルを作って積和演算の部分はマシン語で書いたり、いろいろやりました。
きれいなスペクトルが出ると感動しますよ。ご健闘を祈ります。
 

正しくは、X(k)=1/N{Σx(n)*e^(-j2πkn/N)}
(ただし、k=f*N/f_sample , fは調べたい周波数、f_sampleはサンプリング周波数)
です。
とりあえずは、
a_rl[k]=a_rl[k]+x[n]*cos(w);
a_im[k]=a_im[k]+x[n]*sin(w);
とすれば動作確認はできます。
あとはご自分のプログラムの中でのkの意味(f*N/f_sampleは整数でなくてもよい)、f<(f_sample/2)、に注意し、プログラムを書きなおしてみてください。またきれいなスペクトルを得るには窓関数処理をお忘れなく。

数学、数式が大の苦手の私ですが、DFTはおぼろ...続きを読む

Qフーリエ変換やFFTのプログラム

C言語で書かれているFFTやフーリエ変換のプログラムのあるお勧めのサイトがあればおしえていただけないでしょうか?

よろしくお願いいたします。

Aベストアンサー

No1のepistemeさんが紹介してくださっているライブラリは、それなりに有名だと思うので、
検索すれば情報が得られると思います。

ちなみにライブラリ本体は、
http://www.fftw.org/

http://www.fftw.org/download.html
からダウンロードできます…。


参考:
http://www32.atwiki.jp/amaeda/

参考URL:http://www32.atwiki.jp/amaeda/

QC言語 配列の長さの上限

C言語で配列Array[N]の長さNの上限っていくらなんでしょうか?
もし可能なのであれば上限を2147483647にしたいのですが、方法を教えてください。

Aベストアンサー

そもそもWindowsの32bit版はアプリが仮想メモリ空間を2GBしか使えません。2GBを超えるには64bit版が必要です。
たとえ64bit版OSだとしても添え字が2147483647って、単純なintの配列だとしても4x2147483647=8GB必要ですね。実メモリ16GBとかのPCを用意しますか?
そもそも配列で2147483647個必要なアルゴリズムに問題ありだと思います。

Qフーリエ変換について教えてください

フーリエ変換をすると横軸が時間から周波数になるのはわかったのですが、縦軸が何になるのかわかりません。

一般的に縦軸はなにになるのでしょうか?

また横軸が時間で、縦軸が距離をフーリエ変換したら縦軸は何になるのでしょうか?

よろしくお願いします。

Aベストアンサー

時間関数をフーリエ変換すると結果は、その時間関数の周波数成分が
得られます。スペクトルとも言います。従って、縦軸は、周波数成分です。一般に複素数です。
大きさと偏角による表現もできます。
大きさの方は振幅特性、位相角の方は位相特性と呼ばれます。
画像のように空間座標の上の関数の場合には、フーリエ変換すると
空間周波数成分が得られます。横軸は、空間周波数(2次元)となります。
対象とする関数により結果はそれぞれ意味が異なります。
「一般に何になる」とは言えません。

>横軸が時間で縦軸が距離の場合・・・
フーリエ変換の結果は、距離を表す時間関数の周波数成分です。

フーリエ変換の対象の関数は別に時間関数でなければならないということは
ありません。従って、フーリエ変換の結果は適用する人が解釈(定義)すれば
よいと思います。
たとえば、
時間関数をフーリエ変換し、その結果の絶対値の対数のフーリエ変換を
することもあります。これの結果には、発明者らがケプストラムという名前
をつけています。Cepstrum は Spectrum から作った造語です。

時間関数をフーリエ変換すると結果は、その時間関数の周波数成分が
得られます。スペクトルとも言います。従って、縦軸は、周波数成分です。一般に複素数です。
大きさと偏角による表現もできます。
大きさの方は振幅特性、位相角の方は位相特性と呼ばれます。
画像のように空間座標の上の関数の場合には、フーリエ変換すると
空間周波数成分が得られます。横軸は、空間周波数(2次元)となります。
対象とする関数により結果はそれぞれ意味が異なります。
「一般に何になる」とは言えません。

>横軸が...続きを読む

Q離散フーリエ変換(DFT)の公式について

離散フーリエ変換の公式は、参考書等によりますと色々な記述があります。
今回は2例の違いを教えていただきたいのです。

1)F(u)=1/NΣf(x)・・・・
2)F(u)=Σf(x)・・・・

との式があります(両式とも詳細は省いて記述しました)。この規格化定数(1/N)がある公式1)と、ない公式2)があります。
本来の周波数スペクトル(振幅)を表しているのは1)式であると考えていますが
いかがでしょうか?

公式2)を使用して算出した場合には、その値をサンプリング数で除すれば公式1)同じ結果となるのですが、なぜ公式2)が記述してあるのでしょうか?

Aベストアンサー

戻って来ました。
No.7の補足に対してです。
答は、「そのとおり」です。(^^)

蛇足になりますが、この時の1/Nは形式的には
サンプル数で割っていますが、考えるときは、
F[0]=1が周波数内にしめる割合が1/8と考える方が
分かり易いかもしれません。

何度も書くようですが、周波数が標本化周波数
(標本化間隔の逆数)の1/2を超えるときは注意しましょう。
例の場合、周波数7/8は-1/8と同様に表せるだけでなく
15/8, 23/8,… ,-9/8,-17/8…
としても同様に表せます。
この任意に表せる周波数の内どれが一番妥当かも考えてみましょう。
(参考) 折り返し雑音、エイリアシング(aliasing)

以上です。

Q三角関数の記述の仕方

タイトルそのまんまなんですが、三角関数はC言語ではどのように記述すればいいでしょうか?
角度にラジアン表記でπ(パイ)を使いたいんですが、その表記方法もわかりません。
僕の持っている本に載ってなかったので質問させていただきました。
よろしくお願いします。

Aベストアンサー

C言語で三角関数を使うためには、math.h をインクルードする必要があります。使い方は例えば、こんな感じです。

#define M_PI 3.14159265358979 /* 円周率 */

double x, y, theta;

theta = M_PI / 4.0;
x = cos(theta); /* sin,cos,tanの引数は弧度法の角度です。*/
y = sin(theta);

πは上記の例のように自分で定義して使ってください。

Qファイル出力の場所を指定

現在C++にてhtmlファイルを出力するプログラムを作っているのですが、出力場所を指定することはできるのでしょうか?(現在はそのプログラムソースが保存されている場所と同じファイル内に出力されますが、それをデスクトップに出力するなど。)
もし、方法がありましたら、教えてください。
ソースや参考HPのURLなどのせていただけたらありがたいです。
環境はVisualStudio.NET2003です。
よろしくお願いします。

Aベストアンサー

単にファイル名の前にパスを指定する。

絶対パス指定
fp=fopen("c:/temp/test.txt","w");

相対パス指定
fp=fopen("./hoge/test.txt","w");


デスクトップはOSやユーザによって場所が異なるので、少し面倒です。
XPの場合環境変数を利用してこんな感じで出来ると思います。

例:
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

void main(void)
{
FILE *fp;
char fname[1024];
strcpy(fname,getenv("USERPROFILE"));
strcat(fname,"/デスクトップ/test.txt");
fp=fopen(fname,"w");
//処理
fclose(fp);
}

Q高速フーリエ変換とフーリエ変換の違い

高速フーリエ変換とフーリエ変換の違いについて教えて下さい。
高速フーリエ変換は何か近似を行うことによって、計算速度を速くしているのでしょうか?
もし、何かの極限で出てくる結果が違う場合などがあれば教えて下さい。

Aベストアンサー

>出てくる結果は全く同じだということなのでしょうか?
その通りです。


このQ&Aを見た人がよく見るQ&A

人気Q&Aランキング