これまでもノイズ解析など、主に周波数分解を目的としてFFTを使った解析を行うことはあったが、
もう少しまともな解析が必要になったので、それに際して自分用のメモ。
Contents
FFT の振幅成分
後述する窓関数を回避するため、あるいは測定精度を確保するため、
入力/サンプリング周波数などの測定条件を変えながら測定を行い、個々のデータに対して FFT を適用することがある。
この時、FFT 後の振幅の絶対値について、データ間で比較を行う際には注意が必要。
サンプリング周波数, サンプリング数、実効測定時間などが変わり、
単位スライスあたりの周波数が変わることにより、FFT 後の振幅の絶対値は変わる。
これらのデータを比較したい場合は、
- 入力側も同様の測定・FFT処理を行い、これを用いて規格化を行う
- 単位スライスあたりの周波数幅を用いて規格化を行う
FFT の位相成分
これまで用いることがなかったが、FFT には位相成分も存在する。
これは sin/cos に分解した際に、各周波数の波をどれだけ時間方向にずらして重ねるかという量である。
入力側, 出力側の FFT位相の差分が、周波数空間上で大きな傾きを持っている場合は、
この周波数領域の波形を通すと、出力側で波形が歪んでしまうことになる。
オフセットの除去
例えば、純粋な単色の sin/cos 関数であっても、const. なオフセットを加えると、
この const. なオフセットを低周波数成分で記述しようとするので、
FFT で周波数分解した際には、低周波数領域にゴミのような構造が現れる。
FFT に掛ける前に、定数でフィットを行い(あるいは平均値を計算して)、
このオフセットを取り去ってやると、低周波数領域の余分な成分は除去できる。
窓関数
ある有限サンプル数のデータを FFT に掛ける場合、
暗黙的に、データの左端と右端が繋がった無限長のデータを仮定している。
したがって、両端が滑らかに接続しない場合、FFT はうまく行かない。
例えば単色の sin/cos の場合、
データ区間内にちょうど整数個の波が入っている場合は、正しく単色の強度スペクトルが得られるが、
非整数個の場合は、設定周波数周りにブロードなスペクトルになってしまう。
この問題を回避するために、窓関数というものを用いる。
データ領域の両端で0(に近い値)になるようなウェイト(窓関数)を掛けることで、
滑らかに連続的に接続できるようにしてFFT を行う手法である。
ただし、窓関数により測定データの強度を一部抑えてしまっているため、
データ間の比較を行う際は、振幅の強度の値には注意が必要である。
補正 or 回避方法は以下。
- 入力データの関数系が既知の場合は、窓関数との積をデータ区間に渡って行い、
全体のパワーのうち、窓関数適用後に残るパワーの割合を計算し、これを用いて補正する。 - 入力側も同様の測定・FFT処理(含: 窓関数)を行い、これを用いて規格化を行う
この際、測定データ区間には多数の周期(sin/cos なら波の数)が含まれるようにする。
多数の周期が含まれていると、入力/出力の相対位相が多少ズレていても補正可能。 - そもそも窓関数を使わないようにする。
区間内に整数個の周期が含まれるようにして、測定を行う。
(私は面倒なので、この方法で窓関数を回避する)
実装
ROOT での実装
直近で用いた方法
$ROOTSYS/tutorial/fft/FFT.C を参考にして python で記述。
1 2 3 4 5 6 7 8 |
th=TH1D(...) for k in range(grWfm.GetN()): th.SetBinContent(k,grWfm.GetY()[k]-tOffset) TVirtualFFT.SetTransform(0) hM=th.FFT(hM,"MAG") hP=th.FFT(hP,"PH") #hP.Scale(TMath.RadToDeg()) |
以前用いた方法
何かを参考に書いたC++。
こちらも $ROOTSYS/tutorial/fft/FFT.C ベース?
上記との実装の違いは、ノイズ除去をした後に再度逆変換をしたかったからか?
1 2 3 4 5 6 7 8 9 10 11 12 |
double tArray[sampleLength]; TH1D* hFFT=new TH1D("hFFT","",sampleLength,0.,sampleLength); TVirtualFFT* vFFT[2]; vFFT[0]=TVirtualFFT::FFT(1,sampleLength,"R2HC P K"); vFFT[1]=TVirtualFFT::FFT(1,sampleLength,"HC2R P K"); for ( int n=0;n<sampleLength;n++ ){ hFFT->Fill(n,data[n]); } vFFT[0]->SetPoints(hFFT->GetArray()); vFFT[0]->Transform(); vFFT[0]->GetPoints(tArray); TH1* tHist=(TH1*)TH1::TransformHisto(vFFT[0],NULL,"MAG"); ... |
TH1::FFT の1つめ, TH1::TranformHisto の2つめの引数と、返り値の関係 or 必要性がよくわからない。。。
窓関数の実装
TBF