FFTの結果からフーリエ係数は求められるものでしょうか？

質問

FFTの結果からフーリエ係数は求められるものでしょうか？

実測値の波形データを関数で表したいのです。
フーリエ展開でなんとかならないかと考えたのですが、
計算の速いFFTを使うと、逆変換で波形データに戻せるものの、
関数として表す事が出来ない事に、
FFT後に気付いてしまいました。
数学は専門外で、入門書と格闘しながらやっているのですが、
行き詰ってしまいました。

回答No.1

inara1 · ベストアンサー

データは加速度（振動）データでしょうか。関数で近似できれば、数値でなく式で微分も積分もできるので、積分して速度を出したり、さらに積分して位置を計算したりできますね（そういう目的でしょうか？）。

＞FFTの結果からフーリエ係数は求められるものでしょうか？
波形データの組 A(1), A(2), ... , A(N) をFFTにかけて出てきた数値の組 B(1), B(2), ... , B(N) はフーリエ係数（複素数）ですから、これを係数として元の波形を再現することはできます。ただし、係数B はN個（２の倍数）あるので、N=256 だとしても、１点計算するのに 256項のsinの和を計算しなければなりません。

もし元の波形 A(t) が、ある周波数の基本波（sin波）とその高調波からなるものであれば、以下のような有限項の Fourier級数を仮定して、Excelのソルバーを使って最小２乗近似で係数を出したほうが簡単だと思います（充分な精度ででます）。
　　　A(t) = a(0) + Σ[ k=1～N ] { a(k)*sin( 2*π*k*f*t ) }
未知パラメータは 基本波周波数 f [Hz] と Fourier係数 a(0) ～ a(N) です。a(0) は定数項ですので、振動データの場合はゼロとしてもいかもしれません。

例えば、以下のように、Ｅｘｃｅｌシートの A列の A11からA1011 に時間、その右のB列に対応する波形データがあるとします（1001点のデータがあると仮定しています）。まず、B2～B8 に未知パラメータの値を書きます（この値は適当でいいですが、基本周波数はあまりかけ離れていると収束しないので、ある程度正確な値を記入）。次に、C列に近似波形の計算式を書きます。C11に B2～B8 のパラメータを使って
　　　= $B$3 + $B$4*sin( 2*PI()*$B$2*A11 ) + $B$5*sin( 2*2*PI()*$B$2*A11 ) + $B$6*sin(3* 2*PI()*$B$2*A11 ) + $B$7*sin( 4*2*PI()*$B$2*A11 ) + $B$8*sin( 5*2*PI()*$B$2*A11 )
と書いて Enter を押したあと、このセルをコピーして C12～C1011 に貼り付ければ、C12～C1011 には近似波形のデータが作られます（A列の時間と C列の波形でグラフを描けば波形グラフとなります）。次に元の波形と近似波形の誤差の２乗を D11～D1010 に作ります（ D11 に = (B11-C11)^2 と書いて、そのセルを、D12～D1010までコピー）。最後に、C2 に = sum(D11:D1010) と書いて、ここに差の２乗の和 を作ります。その後ソルバーを使い、差の２乗和が最小となるような B2～B8の値を求めます（その方法は末尾に書きます）。

　　　　A　　　　　　　　　B　　　　　　　　　C　　　　　　　　　D
1　　項目　　　　　　　初期値　　　　差の２乗和
2　　f [Hz]　　　　　　　100　　　　　= sum(D11:D1010)
3　　a0　　　　　　　　　　1
4　　a1　　　　　　　　　　1
5　　a2　　　　　　　　　　1
6　　a3　　　　　　　　　　1
7　　a4　　　　　　　　　　1
8　　a5　　　　　　　　　　1
9
10　時間 [s]　　　　元の波形　　　　　近似波形　　　　差の２乗
11　0.00001　　　　0.021173203　　　= $B$3 + ・・ 　　= (B11-C11)^2
12　0.00002　　　　0.042339618
13　0.00003　　　　0.063492462
　　・・・
1010　0.00999　　　-0.021173203
1011　0.01　　　　　-8.2575E-16

【ソルバーの使いかた】
（１）　Excelのメニューバーの [ツール] → [アドイン] で　「ソルバーアドイン」 にチェックを入れ OK
（２）　Excelのメニューバーの [ツール] → [ソルバー] で
　　　　　　　目的セルを $C$2
　　　　　　　目標値の「最小値」のラジオボタンをクリック（差の２乗和を最小にするという意味）
　　　　　　　変化させるセルを $B$2:$B$8
　　　 とする
（３）　OKをクリックすると近似計算が始まる。収束すれば「解が見つかりました・・」と出るので、OKをクリック
（４）　B2～B8に近似パラメータが記入される（グラフを描いていれば近似グラフも変わります。元の波形も描いてあれば、ほとんど一致していることが分かるはずです）。

回答No.2

inara1

ANo.1 です。
測定データは t=0 から始まっているわけではないので、遅延時間（位相）もパラメータに加えてください。位相を加えた近似式は
　　A(t) = a(0) + Σ[ k=1～N ] [ a(k)*sin{ 2*π*k*f*( t - p ) } ]
です。C11 の近似式は
　　　=$B$4+$B$5*SIN(2*PI()*$B$2*(A11-$B$3))+$B$6*SIN(2*2*PI()*$B$2*(A11-$B$3))+$B$7*SIN(3*2*PI()*$B$2*(A11-$B$3))+$B$8*SIN(4*2*PI()*$B$2*(A11-$B$3))+$B$9*SIN(5*2*PI()*$B$2*(A11-$B$3))
としてください。Excelシートも以下のように位相を加えてください。位相の初期値は、基本波が－側から＋側に行くときのゼロクロス点での時間としてください（高調波が大きくてはっきりしないときはゼロとしてもいいと思います）。

　　　　A　　　　　　　　　B　　　　　　　　　C　　　　　　　　　D
1　　項目　　　　　　　初期値　　　　差の２乗和
2　　f [Hz]　　　　　　　100　　　　　= sum(D11:D1010)
3　　p [s]　　　　　　　0.001
4　　a0　　　　　　　　　　0
5　　a1　　　　　　　　　　1
6　　a2　　　　　　　　　　1
7　　a3　　　　　　　　　　1
8　　a4　　　　　　　　　　1
9　　a5　　　　　　　　　　1
10　時間 [s]　　　　元の波形　　　　　近似波形　　　　差の２乗
11　0.00001　　　　0.021173203　　　= $B$4 + ・・ 　　= (B11-C11)^2
12　0.00002　　　　0.042339618
13　0.00003　　　　0.063492462
　　・・・
1010　0.00999　　　-0.021173203
1011　0.01　　　　　-8.2575E-16

お礼

丁寧な解説ありがとうございます。
なるほど、最小２乗法による近似ですか。
やってみます。

＞フーリエ係数（複素数）ですから、これを係数として元の波形を再現することはできます。

絶対値を算出するまえの、複素数を係数としてsin和を計算すれば良かったんですね。絶対値で計算してました（汗
ありがとうございました。