第31回の講座で、パターン構造を彷彿させるイメージデータを二次元フーリエ変換によって作成する方法を説明した。今回は原点に帰って一般的な１次元フーリエ変換を行ってみよう。テストデータを使ってフーリエ変換を行い、さらにそのスペクトルの大きさと位相の評価を行うので、まずはここからデータをダウンロードしてほしい。ZIPファイルを解凍したら、そのファイルを Load General Text 命令によって IGOR にロードしよう。第１カラムは xx 、第２カラムは yy とでも名づけておけばよい。

　yy を xx に対してプロットすると以下のようなグラフが得られるのだが、フーリエ変換した後の新しい x 軸を設定するのが面倒なので、setscale コマンドによって yy の X スケーリングを変更しよう。

スケーリングについて確認したい場合は、第13回の講座を参照すること。ちなみに、スケーリングを使用せずにそのままフーリエ変換した場合、オリジナル wave のＸ軸は自動的に最大値0.5にスケールされる。この事と、FFT の後にデータの長さが半分になることをわかっていれば、スケーリングを行わずにFFT 後の x wave も自分で作成できるだろう。

setscale/I x xx[0], xx[numpnts(xx)-1], "", yy

numpnts 関数によって xx のデータ点数が得られるため、yy をプロットすると、データは xx の最小値（最初のＸ値）と最大値（最後のＸ値）の範囲に分配される。簡単に言うと、yy vs xx のプロットは今後不要で、yy のみで同等のグラフが描けるということになる。（Y軸 からyy をX軸から "_calculated" を選んでプロット）

　IGOR Proのバージョン５以降では、FFTに様々なオプションが追加され、便利になった。ただ、旧バージョンとの互換性を保つ事と、プログラムの仕組みを理解する意味で一つずつ丁寧に理解していこう。

　バージョン５で可能になったオプションを使わない限り、通常FFTコマンドを実行するとその結果がオリジナルデータに上書きされる。そこであらかじめ FFT の結果を格納するための wave を作成しておこう。

duplicate/O yy yy_FFT

それでは複製した wave に対して FFT の実行だ。IGOR Pro は FFT プログラムを搭載しているので、わざわざ自分でプログラミングする必要はない。これまでに様々な FFT プログラムが公開されているが、IGOR は、Numerical Recipes in C に掲載されているアルゴリズムを採用している。FFT のプログラム原理が気になる人は原著を参照のこと。

FFT yy_FFT

さて、フーリエ変換されたデータは複素数ウェーブだ。スペクトルの大きさを求めるため、データ格納用の wave を作成する。データの長さに注意して（FFTするとデータ長が半分になる）、yy ではなく、yy_FFT をコピーしよう。

duplicate/O yy_FFT yy_mag

振幅は、実数部(real)の二乗と虚数部(imag)の二乗の和を計算し、その平方根をとることで求められる。つまり、(A^2+B^2)^0.5だ。今回は複素数データの共役をかけて平方根をとった。((A+iB)*(A-iB)）^0.5）という計算だが、同じ事だ。ちなみにIGOR において、共役関数は、conj()で与えられる。

yy_mag=(yy_FFT*conj(yy_FFT))^0.5

得られた wave の虚数データはもはや不要なので、yy_mag の虚数部を捨て、実数 wave に変換しよう。Redimension 命令は、wave の長さ、次元、計算精度など、様々なデータ変換操作が行えるが、/Rフラッグを使えば実数 wave に変換する機能をもつ。

Redimension/R yy_mag

yy_mag をプロットして下の赤丸ようなスペクトルが得られたら成功だ。

ちなみに、yy_mag が実数 wave なら、

yy_mag=real(r2polar(yy_FFT))

と入力しても同様の結果が得られるはずだ。r2polar 命令は座標変換のための命令で、この関数で変換された wave の実数部はスペクトルの大きさを与える。上の図で黒い実線がその結果だ。

　さて、次は位相だ。データ格納用の wave として、yy_FFT の複製を用意し、この wave を実数に変換しよう。

duplicate/O yy_FFT yy_phase

Redimension/R yy_phase

複素平面上でのデータは（cos（位相角）, sin（位相角））で表されるので、位相角は（虚数部／実数部）のアークタンジェントで与えられる。IGOR には atan 関数があるが、この関数を使うと結果は、-pi/2 から +pi/2 の範囲で変換されてしまう。そこで、-pi から +pi の範囲の値を正常に返してくれるatan2 関数を使おう。atan 関数と少々書式が異なるので注意が必要だが、以下のコマンドにより位相に変換できる。

yy_phase= atan2(imag((yy_FFT)),real((yy_FFT)))

ちなみに以下のように、r2polar 関数で変換したデータの虚数部をとっても同様の結果が得られるはずだ。

yy_phase=imag(r2polar((yy_FFT)))

　最後に unwrap について説明しよう。位相データは、-pi から +pi の範囲に制限されるため、連続的に増大するデータを作成するためには、２πのギャップを補正する必要がある。下図はこれを説明するための例だ。

上の図の一番下はある周期のsinカーブであり、真ん中の図がサインカーブ上の赤丸の点の位相だ。真ん中の図で、0, (2/5)pi, (4/5)pi, の次は、(6/5)pi となるはずが、突然負の値(-4/5)pi となっている。こういった位相の性質を考慮して、所定の範囲を飛び越えたデータを補正し、スムーズな位相データに変換するのがunwrap 命令だ。unwrap すると、上段の図のように正しく位相が評価できる。

unwrap 2*pi, yy_phase

yy_phase をプロットしたら以下のようなグラフが得られただろうか？

当然だが、unwrap する前のデータが思っていたようなスムーズなデータと違ったからと言って、下の図のように「絶対値をとって積分する」などという誤った計算をしてはならない。abs(-4pi/5)は、6pi/5とは異なるからだ。

最後に余談だが、FFT直後のデータ(yy_FFT)をプロットすると以下のよう実数部と虚数部の両方が表示される。

ここでプロットされたデータ上をダブルクリックすると、「Modify Graph Appearance」パネルが出現する。図中の赤枠で囲んだ部分を変更する事で、実数部のみ、虚数部のみ、振幅、位相など様々なプロットを得ることができる。

log 表示と同様、データが変わるわけではないので注意が必要だが、データをプレビューしたい場合には便利だ。