みその計算物理学
ホーム はじめに リンク集
フーリエ変換とスペクトルの強さの数値計算例(C言語)
  • スペクトルの強さとは

    スペクトルの強さとは、ある関数 f(x) をフーリエ変換したときにでてくる関数 F(ω) の絶対値をとった大きさである。また関数を f(x) をフーリエ級数で表した場合、それぞれの角振動数の成分の係数の2乗の大きさである。

    ここでは、簡単な関数をフーリエ変換してスペクトルの強さを求めてみる。この作業は自分の作ったフーリエ変換のプログラムが正しいかの確認にもなる。

  • フーリエ変換を行う関数

    ここでは、関数 f(x)=2sin(2ωx)+5sin(4ωx) をフーリエ変換する。また角振動数 ω は円周率πに等しいとする。基本的に離散的フーリエ変換の数値計算例(C言語)で行うことと一緒である。よってプログラムも主要な部分はほとんど変更はない。

  • 関数のデータを求めるプログラム例
    #include <stdio.h>
    #include <math.h>
    
    #define pi 3.1415926535
    
    int main()
    {
            double x,y,w;
    
            w=pi;
    
            FILE *function;
            
            function=fopen("function.data","w");
    
            for(x=0.0;x<=10.0;x+=0.01) {
                    y=2*sin(2*w*x)+5*sin(4*w*x);
                    fprintf(function,"%f\n",y);
            }
    
            fclose(function);
    
            return 0;
    }
  • フーリエ変換してスペクトルの強さを求めるプログラム例
    #include <stdio.h>
    #include <math.h>
    
    #define max 10000                       //Nの限界値
    #define pi 3.1415926535         //円周率
    
    int main()
    {
            int k,n,N;
            double f[max+1],ReF,ImF;
    
            FILE *function;
            FILE *fourier;
    //フーリエ変換したいデータをあらかじめファイル function.data に保存しておく
            function=fopen("function.data","r");
            fourier=fopen("fourier.data","w");
    
    //データの読み込み(ただし実数の成分のみ)
            for(N=0;N<max;N++) {
                    if(fscanf(function,"%lf",&f[N]) == EOF) {
                            N--;
                            break;
                    }
            }
    
    //実数部分と虚数部分に分けてフーリエ変換
            for(n=0;n<N;n++) {
                    ReF=ImF=0.0;
                    for(k=0;k<N;k++) {
                            ReF+=f[k]*cos(2*pi*k*n/N);
                            ImF+=-f[k]*sin(2*pi*k*n/N);
                    }
                    fprintf(fourier,"%d %f %f %f\n",n,ReF,ImF,ReF*ReF+ImF*ImF);             //最後にスペクトルの強さを追加
            }
    
            fclose(function);
            fclose(fourier);
            
            return 0;
    }
  • GNUPLOTでの出力結果

    GNUPLOTでスペクトルの強さのみ出力する。

  • スペクトルの強さの分析

    2つある画像のうち重要なのは上の方の画像である。下の画像はコンピュータでは極限をとれないことが原因で大きい角振動数が生じたので参考に載せた。なので誤差として無視してよい。

    さて、上の画像をみると、横軸が10と20のところに鋭いピークがある。つまり10番目と20番目の角振動数においてスペクトルの強さが強いということになる。具体的にこのときの角振動数を計算すると、離散的フーリエ変換(PDF形式)に記述した角振動数 2πn/Nh から、上記のプログラム例で x の刻み幅は 0.01 、得られたスペクトルの強さのデータの個数から N を約 1000 に等しいとして、

    2X3.14X10 ÷ (1000X0.01) = 6.28 = 2π
    2X3.14X20 ÷ (1000X0.01) = 12.56 = 4π

    となる。確かにフーリエ変換した関数が持つ角振動数と全て一致している。

    ではスペクトルの強さはどうだろうか?予想として f(x) のそれぞれの角振動数の係数の2乗の比から 4:25 となるはずである。GNUPLOTでの出力結果からはわかりにくいので実際の出力データを見てみると、

    0 0.000009 0.000000 0.000000
    1 -0.000001 0.000000 0.000000
    2 -0.000001 0.000001 0.000000
    3 -0.000000 0.000001 0.000000
    4 0.000000 0.000001 0.000000
    5 0.000001 0.000000 0.000000
    6 0.000001 -0.000000 0.000000
    7 0.000001 -0.000001 0.000000
    8 0.000001 -0.000001 0.000000
    9 0.000000 -0.000001 0.000000
    10 0.000004 -1000.000028 1000000.056538
    11 -0.000001 -0.000000 0.000000
    12 -0.000001 -0.000000 0.000000
    13 -0.000001 0.000000 0.000000
    14 -0.000001 0.000001 0.000000
    15 -0.000000 0.000001 0.000000
    16 0.000000 0.000001 0.000000
    17 0.000001 0.000000 0.000000
    18 0.000001 -0.000000 0.000000
    19 0.000001 -0.000002 0.000000
    20 -0.000002 -2499.999987 6249999.934255
    21 0.000000 0.000001 0.000000
    22 -0.000000 -0.000000 0.000000
    以下続く

    一番右側の列がスペクトルの大きさを表す。10と20の欄の数値から比を求めるとほぼ 4:25 である。よってこれも理論的な比とほぼ一致している。

    この一連の流れでフーリエ変換のイメージがつかめただろうか。実際に数値計算を行えばフーリエ変換の意味がつかめてくると思う。