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

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

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

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

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

  • 関数のデータを求めるプログラム例
    import java.io.*;
    
    class Functiondata
    {
            public static void main(String arg[])
            {
                    double x,y;
                    x=y=0.0;
    
                    try
                    {
                            FileWriter fw=new FileWriter("function.data");
    
                            for(x=0.0;x<=10.0;x+=0.01)
                            {
                                    y=2*Math.sin(2*x*Math.PI)+5*Math.sin(4*x*Math.PI);
                                    fw.write(y + "\n");
                            }
                            fw.close();
                    }
                    catch(Exception e)
                    {
                            System.out.println(e);
                    }
            }
    }
    
  • フーリエ変換してスペクトルの強さを求めるプログラム例
    import java.io.*;
    
    class Spectrum
    {
            public static void main(String arg[])
            {
                    int k,n,N,max;
                    N=0;
                    max=10000;      //配列の容量を決める数(データの個数に応じて余裕を持たせておく)
                    double f[]=new double[max+1],ReF,ImF;
    
                    try
                    {
    //フーリエ変換したいデータをあらかじめファイル function.data に保存しておく
                            FileReader fr=new FileReader("function.data");
                            FileWriter fw=new FileWriter("fourier.data");
                            StreamTokenizer st=new StreamTokenizer(fr);
    
                            while(st.nextToken()!=StreamTokenizer.TT_EOF)
                            {
                                    if(N>max)
                                            break;
                                    f[N]=st.nval;
                                    N++;
                            }
    
                            //実数部分と虚数部分に分けてフーリエ変換
                            for(n=0;n<N;n++) 
                            {
                                    ReF=ImF=0.0;
                                    for(k=0;k<N;k++) 
                                    {
                                            ReF+=f[k]*Math.cos(2*Math.PI*k*n/N);
                                            ImF+=-f[k]*Math.sin(2*Math.PI*k*n/N);
                                    }
    
    //最後にスペクトルの強さを追加
                                    fw.write(n + " " + ReF + " " + ImF + " " + (ReF*ReF+ImF*ImF) + "\n");
                            }
    
                    }
                    catch(Exception e)
                    {
                            System.out.println(e);
                    }
            }
    }
    
  • 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 -52.420680426559009 0.0 2747.9277363834267
    1 25.40028590318969 -41.377809184938094 2357.2976169089243
    2 16.604798259233405 -24.170212729361889 859.91850861244825
    3 -0.89431323945362584 -19.100797230635198 365.64025101610326
    4 -12.478301237626592 -10.590058586519373 267.85734264286606
    5 -8.7467276530201232 -7.3471112557273743 130.48528844014277
    6 -0.503167680759154 -16.127874931378031 260.36152751713246
    7 -0.023415197660089815 -24.484401939781112 599.48648662023845
    8 -1.0955888581423205 -19.543123697553504 383.13399880396292
    9 4.4965077880292528 -8.9780747411361332 100.82440834523436
    10 55.691975823745139 -1011.6348673932472 1026506.7010969055
    11 3.8035066147768859 5.2696849762386178 42.236242317446532
    12 -4.0415864955007139 14.152421202497447 216.62544729351302
    13 -3.847447228848818 18.120003188047146 343.13736571361522
    14 -2.7761110620657519 8.5213181570811383 80.319655763124516
    15 -9.0799250966792577 -2.8928194421643 90.8134440862696
    16 -12.705564305149846 -2.2766446146921049 166.61447501390444
    17 -2.0946076472409447 4.3907612247825458 23.666165328934166
    18 14.821778598919446 4.3628670301157113 238.71972955785719
    19 22.129990500075415 -2.1232355188209988 494.24460860181119
    20 258.46814034512045 -2514.1958483233975 6387986.5433000736
    21 26.834026757978254 15.809831730238654 970.01577138635389

    以下続く

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

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