資料集：Igorの使い方

9. 実際的な使用例 2 (Power spectrum)

周波数解析は、いろいろな分野で使用される解析方法です。フーリエ変換してどのような周波数の成分が多く含まれるかを調べること、という程度の理解はあっても、実際に自分で行うにはなかなか閾が高そうに思えます。しかし実際に行ってみるととても簡単です。

ここでは、まず単純な周波数のデータを作成し、Igor Proによりフーリエ変換を行い、どのような結果が出てくるかを検討します。

以下の条件でテスト用データを作成します。
sampling rate 128 Hz
data length 4096
3 Hz、振幅3のcosine波、5 Hz、振幅5のsine波、 10 Hz、振幅10のcosine波の合わさった波を作成。
Igor Proを起動して
make/o/d/n=4096 a
variable f = 128
variable dt = 1/f
SetScale/P X 0, dt, "s", a
display a
時間軸が0から32までの直線が現れます。
テストデータを作成
variable dw = 2 * pi
a = 3*cos(3*dw*x) + 5*sin(5*dw*x) + 10 * cos(10*dw*x)
少し拡大すると、このような波形が出来ているのがわかります。
さて、いよいよフーリエ変換です。
通常のフーリエ変換では、実数部分と虚数部分に結果が出てきます。虚数!？と思うかもしれませんが、位相がずれたものという程度に理解しておけばよいと思います。
（三角関数の加法定理によると、位相がaずれたcosine波は、cos(x-a) = cos(a)*cos(x)+sin(a)*sin(x)で表されます。ここでcos(a)とsin(a)は定数と考えれば、cos(x)とsin(x)を適当に混ぜ合わせると、任意の位相の波をつくことが出来ることを示しています。）

FFT /OUT=1 /DEST=aout a
オプションの/OUTはFFTのmodeを指定するもので、1は通常の実数+虚数、2は実数のみ、3はmagnitude、4はmagnitudeの2乗を出力します。出力先は、/DESTで指定します。
早速、結果を表示してみます。
display aout

周波数が、3, 5, 10のところにピークが見られます。これは期待通り。しかし振幅が桁違いに大きいです。またcosine波とsine波の区別がつきません。
フーリエ変換の注意点の一つに、scalingがあります。このIgorのFFTの場合では、データ数の半分の数（4096/2 = 2048）で割ります。
aout = aout / 2048
そうすると、それぞれのピーク値は、3, 5, 10となっています。これでフーリエ解析の意味が何となくつかめてきたと思います。
オプション/OUT=1を使用した場合、実はaoutは複素数のwaveです。テーブルで確かめるには、
メニュからWindow -> New Table ...とするか、もしくはコマンドから
edit aout
テーブルを見ると、aoutは2049の複素数からなるwaveであり、テストデータのsine波にあたるところは、point 160のところで、虚数 5.0となっているのが確認できます。
次にオプションの/OUT=2を使用してみます。これは実数部分のみを結果として残します。
FFT /OUT=2 /DEST=aout a
aout = aout/2048
aoutのグラフを見ると、sine波の成分がありません。ちなみに、この場合、aoutは実数のみのwaveでサイズは2049です。
更にオプションの/OUT=3を使用してみます。これはmagnitudeを結果として残します。magnitudeは、実数部をR、虚数部をIとした場合、|R|²+ |I|²の平方根の値です。
FFT /OUT=3 /DEST=aout a
aout = aout/2048
aoutのグラフを見ると、/OUT=1の場合と似ています。しかし上の場合と異なり、aoutは実数のみのwaveです。
オプション/OUT=4は、magnitude powerを出力します。基本的にはmagnitudeの平方ですが、scalingを2048²で行わなくてはなりません。
FFT /OUT=4 /DEST=aout a
aout = aout/(2048*2048)
それぞれのピークの値が、（ほぼ）テストデータの振幅の2乗の値になっているのがわかると思います。

以上が、IgorProを用いた、フーリエ解析、パワースペクトルム解析の基本でした。
では、周波数解析の意味をもう少し理解するために、テストデータにノイズを加えます。
a = 3*cos(3*dw*x) + 5*sin(5*dw*x) + 10 * cos(10*dw*x) + gnoise(100)
gnoise(v)は、標準偏差sdがvの正規分布となる乱数を発生させる関数です。
パワースペクトルムを計算すると、周波数10のところにピークは見えますが、他の成分はわかりません。
しかしパワースペクトルムを10回加算平均すると、他の成分もなんとか見えてきます。
加算平均の部分を含めたスクリプトをのせておきます。

#pragma rtGlobals=1 // Use modern global access method.

function powerSpec(s, n, f)
  wave s
  variable n
  variable f

  variable nsq1 = 4/(n*n)
  variable dt = 1.0/f
  variable dw = 2.0 * pi

  make/o/d/n=(n) a SetScale/P x 0, dt, "", a
  a = 3.0*cos(3.0*dw*x) + 5.0*sin(5.0*dw*x) + 10.0*cos(10.0*dw*x) + 10.0*cos(20.0*dw*x+5)
  a = a + enoise(200)

  FFT /Dest=pws /OUT=4 a
  pws = nsq1 * pws

  s = s + pws
end

function averageSpectrum()
  variable n = 4096 // number of data points
  variable f = 128 // sampling rate

  variable iter = 10  // number of averaging
  variable ix
  variable n2 = n/2

  make/o/d/n=(n2+1) sw        // output
  SetScale/P x 0, f/n, "", sw // scaling in frequency domain
  sw = 0
  for(ix=0;ix< iter;ix+=1)    // loop
    powerSpec(sw, n, f)
  endfor

  sw = sw/iter               // averaging

end
Igor Proは計算速度が速いので、単純な計算の場合計算速度が問題となることはありません。しかし長時間記録の周波数解析を行う場合には、速度が問題となってくる可能性があります。そのような場合は、C言語等を用いてプログラムを作成する必要があるかもしれません。フーリエ変換に関しては優れたライブラリがあるので、プログラム言語で周波数解析行うことも、格段に難しいことではありません。この件に関しては、いずれ検証を行いたいと思います。

普通のデスクトップ機（Core 2 Duo 2.8GHz, Windows 7 64bit）でおおよその時間を計測しました。上記の4096データポイントの計算を3600回繰り返したところ、約10秒かかりました。1秒ごとのデータを解析するとして、これは1時間のデータの相当するので、かなりの規模のデータまでIgor Proで解析できるのではないかと思います。