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と関連する良く見られている質問

Q離散フーリエ変換

今、離散フーリエ変換の値が求まっています。
これから、振幅の値を出すのは、どうしたらいいのでしょうか?

自分で調べたところ、離散フーリエ変換の値に標本化関数のフーリエ変換をかけて、サンプリングの間隔で割ればいいのでは、と考えているのですが、標本化関数というものがよくわかりません。このやり方で良いのかもわかりません。
アドバイスお願いします。

Aベストアンサー

元の信号の式っていう事ですよね?:

離散フーリエ変換式と逆変換式を調べて並べて補足してください

Qフーリエ変換のC言語プログラムについて

正弦波(およびガウス性雑音)をフーリエ変換(離散)→逆フーリエ変換するというプログラムを組みました。正弦波をフーリエ変換すると実部は2回ピークがくるはずなのですが、すべて「0.000000」または「-0.000000」と表示されてしまいます。虚部は正常なようで実装の仕方もさほど違わないので、何が問題なのかわからずにいます。念のためコードはすべて載せますが、該当箇所は関数Fourierの
fp = fopen("reX.txt", "w");//書き込む
あたりです。問題点を教えていただけないでしょうか。お願いします。


//gauss.txt, sin.txt:発生させたガウス性雑音&正弦波
//reX, imX:フーリエ変換の実部と虚部
//re-1, im-1:逆フーリエ変換の実部と虚部

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

#define PI3.14159265358979323846
#define N256

//n:長さ, w:角周波数, p:位相(phase), a:振幅
void SinCurve(int n, double w, double p, double a)
{
FILE *fp;
double x;
int t;

fp = fopen("sin.txt", "w");//書き込むので"w"
if(fp == NULL)
{
printf("file open error\n");
} else
{
for(t = 0; t < n; t++)
{
x = a * sin( w*(double)t + p );
fprintf(fp, "%f\n", x);
}
}

fclose(fp);
}

//n:長さ, s:分散, m:平均
void Gauss(int n, double s, double m)
{
FILE *fp;
double x, x1, x2, y1;
int t;

srand((unsigned) time(NULL));

fp = fopen("gauss.txt", "w");//書き込むので"w"
if(fp == NULL)
{
printf("file open error\n");
} else
{
for(t = 0; t < n; t++)
{
x1 = ( (double)rand() + 1.0 ) / ( (double)RAND_MAX + 2.0);
x2 = ( (double)rand() + 1.0 ) / ( (double)RAND_MAX + 2.0);
y1 = pow(-2.0*log(x1), 0.5) * cos(2.0*PI*x2);
y1 = s * y1 + m;
fprintf(fp, "%f\n", y1);
}
}

fclose(fp);
}

//ファイル名sのデータをフーリエ変換し、実部と虚部をreX, imXに保存
void Fourier(int num, char *s)
{

FILE *fp;
int k, n;
double largeX, x[N+100], t;

fp = fopen(s, "r"); //読み込み
if(fp == NULL)
{
printf("file open error\n");
} else
{
//printf("%s\n", s);
for(k = 0; k < num; k++)
{
fscanf(fp, "%lf", &x[k]);
printf("x[%d]=%f\n", k, x[k]);
}
}

fp = fopen("reX.txt", "w");//書き込む
if(fp == NULL)
{
printf("file open error\n");
} else
{
for(k = 0; k < num; k++)
{
largeX = 0.0;
t = 2.0*PI*(double)k / (double)N;
for(n = 0; n < num; n++)
{
largeX += x[n] * cos((double)n*t);
//printf("%f\n", largeX);
}
fprintf(fp, "%f\n", largeX);
printf("reX[%d]=%f\n", k, largeX);
}
}

fp = fopen("imX.txt", "w");//書き込む
if(fp == NULL)
{
printf("file open error\n");
} else
{
for(k = 0; k < num; k++)
{
largeX = 0.0;
t = 2.0*PI*k / (double)N;
for(n = 0; n < num; n++)
{
largeX -= x[n] * sin(n*t);
}
fprintf(fp, "%f\n", largeX);
}
}

fclose(fp);
}

void InverseFourier(int num)
{
FILE *fp;
int k, n;
double a[N+100], b[N+100], x, t;
//a:reX, b:imX

fp = fopen("reX.txt", "r"); //読み込み
if(fp == NULL)
{
printf("file open error\n");
} else
{
for(k = 0; k < num; k++)
{
fscanf(fp, "%lf", &a[k]);
//printf("a[%d]=%f\n", k, a[k]);
}
}

fp = fopen("imX.txt", "r"); //読み込み
if(fp == NULL)
{
printf("file open error\n");
} else
{
for(k = 0; k < num; k++)
{
fscanf(fp, "%lf", &b[k]);
//printf("b[%d]=%f\n", k, b[k]);
}
}

fp = fopen("re-1.txt", "w"); //読み込み
if(fp == NULL)
{
printf("file open error\n");
} else
{
for(n = 0; n < num; n++)
{
x = 0.0;
t = 2.0*PI*(double)n / (double)N;
for(k = 0; k < num; k++)
{
x +=a[k] *cos(k*t) - b[k] *sin(k*t);
}
x /= (double)N;
fprintf(fp, "%f\n", x);
//printf("x[%d]=%f\n", n, x);
}
}
/*
fp = fopen("im-1.txt", "w"); //読み込み
if(fp == NULL)
{
printf("file open error\n");
} else
{
for(n = 0; n < num; n++)
{
x = 0.0;
for(k = 0; k < num; k++)
{
t = 2.0*PI*(double)k / (double)N;
x = x + a[k] *sin(n*t) + b[k] *cos(n*t);
}
x /= (double)N;
fprintf(fp, "%f\n", x);
}
}
*/

fclose(fp);
}

int main(void)
{
SinCurve(N, PI/8.0, 0.0, 1.0);

//Gauss(N, 1.0, 0.0);

Fourier(N, "sin.txt");

//Fourier(N, "gauss.txt");

InverseFourier(N);

return 0;
}

正弦波(およびガウス性雑音)をフーリエ変換(離散)→逆フーリエ変換するというプログラムを組みました。正弦波をフーリエ変換すると実部は2回ピークがくるはずなのですが、すべて「0.000000」または「-0.000000」と表示されてしまいます。虚部は正常なようで実装の仕方もさほど違わないので、何が問題なのかわからずにいます。念のためコードはすべて載せますが、該当箇所は関数Fourierの
fp = fopen("reX.txt", "w");//書き込む
あたりです。問題点を教えていただけないでしょうか。お願いします。


//gauss.txt...続きを読む

Aベストアンサー

元のデータがsin関数なら、虚部だけしか出ないので正常では?

fopen()に対応したfclose()が必要なのはその通りなので直すべきですが。

Q離散フーリエ変換の周期とサンプリング間隔と周波数

離散フーリエ変換で求まるスペクトルの各点の周波数について質問があります。
離散フーリエ変換で時間軸上の各点(0~T[s]でΔt刻みにN個の点を取った)を周波数軸上の各点に変換したときの周波数の換算式を調べると、
Δf=1/Tとなっていたり、Δf=1/(N*Δt)
となっていました。
意味上はどちらでも良さそうな気がしたのですが、実際に計算してみると両者の式で周波数軸上の各点での周波数がずれていました。
たとえば0秒から0.1秒刻みで10点とると一周期T=0.9秒になるのですが、N*Δtで計算すると一周期1秒になってしまいます。0.9秒しか見ていないのに一秒周期の関数としてフーリエ変換していることになると思いますが、周波数間隔はどちらの式で計算すべきでしょうか?それとも用いるフーリエ変換の式によって異なるのでしょうか?
教えていただければ助かります。よろしくお願いします。

Aベストアンサー

0~T[s]でΔt刻みにN個の点を取った
の意味ですが、
0からT秒までをΔtの幅で分割します。
分割した区間内のある時刻に計測するわけです。
たとえば、すべての区間で最初に測るときもあれば、
すべての区間で最後に測るときもあるでしょう。
 N回測ったときにかかる時間は、
それぞれの区間で計測した瞬間だけではなく、
各区間の残りの時間を含めて考えるべきです。
 そうしないと、
(N+1)回目の測定は、2*Δt秒間で1回計測することになってしまいます。
均等にするには0.9秒の瞬間に測ったとしても残りの0.1秒をつけた区間の中で1回測ったとして扱いますので、10回測るのに必要な時間は1秒となります。

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離散フーリエ変換(DFT)について。

離散フーリエ変換(DFT)について。
次の有限長N=4のディジタル信号の離散フーリエ変換(DFT)の周波数スペクトルを求めよ。[F[0],F[1],F[2],F[3]]=[-1,1,-1,2]
について。
算出した所、
DFTは
F[0]=1
F[1]=j
F[2]=-5
F[3]=-jと算出できましたが正解でしょうか。

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

Aベストアンサー

逆変換して検算すれば良いですね。DFTの定義によっては定数倍の係数がつきますんで、ご確認を。係数のことを無視すれば、合ってるっぽいです。

QC言語でフーリエ変換を作成したいです。

(d^n/dx^n)f(x)=1/2π∫(∞to-∞)(ik)^nF(x)*e^ikxdk
をC言語で表したいのですが、上手くいきません。
nは実数でF(x)はf(x)の微分です。

アイデアだけでもいいので教えてください。

Aベストアンサー

フリーのライブラリを一つ紹介します。

汎用 FFT (高速 フーリエ/コサイン/サイン 変換) パッケージ
Copyright Takuya OOURA, 1996-2001

設計方法も説明されていますので、参考にされてはいかがでしょうか。

参考URL:http://momonga.t.u-tokyo.ac.jp/~ooura/fft-j.html

Q離散フーリエ変換について

離散フーリエ変換によって得られた値についての質問です。
多くのサイトでその値は
Σ(k=0~N-1) f(k)exp(-2πkni/N) という式から求められるとあります。
離散フーリエ変換は本来、ある周期関数が、どのくらいの振幅でどのくらいの周波数の波からできているかを調べるために行うものだと思います。
しかし上記の公式から得られるスペクトル(sqrt(Re^2+Im^2))では振幅の値は得られません。振幅を得るには刻み幅Δ(関数をサンプリングした際の幅)を乗じて
Σ(k=0~N-1) f(k)exp(-2πkni/N)*Δ
とすれば得られることが分かりました。
最初の公式から得られるスペクトルはなにを表しているのでしょうか?またなぜ刻み幅Δを乗じることで、振幅が求まるのでしょうか?
よろしくお願いします。m(__)m

Aベストアンサー

パッと見だけど, Δを掛けた式で N → ∞ の極限をとると連続の Fourier 変換にならないかなぁ?

QC言語プログラムの質問です。 実数をxを読み込み次の計算をするCプログラムを作成し、そのプログラムリ

C言語プログラムの質問です。
実数をxを読み込み次の計算をするCプログラムを作成し、そのプログラムリストを記しなさい。
2sin(x)cos(x) および sin(2x)
次にこのプログラムを用いて、x=0.785を計算しなさい。

画像のプログラムを作成し、計算をしたのですが、計算結果が全て0.00000となってしまいます。
どこが間違っているか教えてください!

Aベストアンサー

scanfを以下のように変えてください。
scanf("%lf", &x);

Q離散フーリエ変換

連続関数f(x)をフーリエ変換してときの係数をcnとすると、
連続関数からサンプル点での値をとってf(xj)の離散フーリエ変換したときの係数c'nとすると両者は同じものなのでしょうか?サンプル点の数とかによるのかもしれませんが、c'nというのは厳密なものなので同じになる気もします。
よくわからないので教えてください.お願いします。

Aベストアンサー

「ナイキストのサンプリング定理」の話ですね。

以下、めんどくさいから定数倍を無視して書きますと、

fをフーリエ変換したときに得られるのは(係数cnなんかじゃなくて)関数F(ω)であり、
f(x) ~ ∫ F(ω) exp(-iωx) dω (積分はω=-∞ ~ ∞)
となる。特にfが周期関数の場合には、δ関数の列
F(ω)=Σc[n] δ(ω-nΔω) (総和はn=-∞ ~ ∞)
が得られ
f(x) ~ ∫ F(ω) exp(-iωx) dω = Σc[n] exp(-inΔωx) (総和はn=-∞ ~ ∞)
となります。

 一方、fが周期関数であるとき、これをフーリエ級数展開すると係数c[n]が得られ、
f(x) ~ Σ c[n]exp(inωx) (総和はn=-∞ ~ ∞)
です。この場合、fのフーリエ変換はδ関数の列
F(ω) = Σc[n]δ(ω-nW)(総和はn=-∞ ~ ∞)
となります。
 
 ところで、離散フーリエ変換ってのは、「f(x)が周期関数であって、それをサンプリングしたものをフーリエ変換する」という場合の話です。

 サンプリングをするというのは、サンプリング関数をs(x)と書いて
f(x)s(x)
を作るということ。ここに、
s(x) = Σδ(x-mΔx) (総和はm=-∞ ~ ∞)
です。sのフーリエ変換は
S(ω) = Σδ(ω-nW) (総和はn=-∞ ~ ∞)
ここにWはサンプリング周波数。
従って、f(x)のフーリエ変換をF(ω)とすると、f(x)s(x)のフーリエ変換は
F(ω)*S(ω) (*は畳み込み積分)
です。Sの式から分かるように、これは周期Wの周期関数であり、すなわち、
F(ω+kW)*S(ω+kW) = F(ω)*S(ω) (kは整数)
を満たしています。

 特に、もしf(x)がサンプリング周波数の半分(W/2)(これを「ナイキスト周波数」と言います)以上の周波数成分を全く含まない場合に限ると、
|ω|≦(W/2)のとき、F(ω)*S(ω) = F(ω)
となりますから、この場合には
G(ω) = |ω|≦(W/2)のとき1、それ以外の時0
という関数を掛け算してやると、全てのωについて
(F(ω)*S(ω))G(ω) = F(ω)
となり、f(x)のフーリエ変換と一致する。
両辺を逆変換すれば、
Σf(nΔx)δ(x-nΔx) *sinc(x/(πΔx)) = f(x) (総和はn=-∞ ~ ∞, sinc(x)=sin(x)/x )
である。つまり、fをサンプリングしたものから、f全体を再現できるわけです。

 さて、以上の話では、f(x)が周期関数かどうかは関係ない。なので、離散フーリエ変換(「f(x)が周期関数であって、それをサンプリングしたものをフーリエ変換する」という場合の話)でも同じように成り立ちます。
 ですから、「もしf(x)がナイキスト周波数(W/2)以上の周波数成分を全く含まない」ならば、ご質問の答はyesである。さもなくばnoである。

「ナイキストのサンプリング定理」の話ですね。

以下、めんどくさいから定数倍を無視して書きますと、

fをフーリエ変換したときに得られるのは(係数cnなんかじゃなくて)関数F(ω)であり、
f(x) ~ ∫ F(ω) exp(-iωx) dω (積分はω=-∞ ~ ∞)
となる。特にfが周期関数の場合には、δ関数の列
F(ω)=Σc[n] δ(ω-nΔω) (総和はn=-∞ ~ ∞)
が得られ
f(x) ~ ∫ F(ω) exp(-iωx) dω = Σc[n] exp(-inΔωx) (総和はn=-∞ ~ ∞)
となります。

 一方、fが周期関数であるとき、これをフーリエ級数展開すると...続きを読む

Q離散フーリエ変換について

閲覧して頂きありがとうございます。

現在、課題で
http://www.cryptrec.go.jp/estimation/rep_ID0211.pdf
の3.6の検定をC言語で作成しようとしていますが、
Step2とStep5に記載されている数式が理解できずCのコードに書き直せません。

何か参考になるホームページやソースコードなどあれば教えて頂けないでしょうか。
漠然とした内容で申し訳ありませんが、よろしくお願いします。

Aベストアンサー

Step3???
むしろ Step2 よりはるかに簡単なはずなんだけど... T と N0 を計算するだけ, ですよ.


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

人気Q&Aランキング

おすすめ情報