◎フーリエ計算ソフト
等間隔の数値列から、その数値列の正弦波成分の周波数、
振幅、位相を算出するためのソフト群です。

○内容物
fftcyc.exe:高速フーリエ変換で、データに含まれる周波数の候補を列挙します。
dftcyc.exe:fftcyc.exeで算出した周波数の候補が実際にあるかを確認し、
           その振幅と位相を計算します。
subtract.txt:データ補正ファイルのサンプル

es2000.txt:サンプルデータ
  （JPL-Horizons,太陽位置AD2000-2500分4日間隔,平均黄道座標)
 *) 1- 7行目:2000年1月1日からの通日(時刻は0h-UTC)
 0) 9-18行目:太陽黄経−平均太陽黄経(秒)
 1)22-29行目:太陽黄経−平均太陽黄経−(359.9905度/年の成分)−(719.981度/年の成分)(秒)
 2)33-39行目:太陽黄緯(秒)
 3)43-51行目:太陽地球距離(mAU)
 4)56-62行目:太陽地球距離−平均距離−(359.9905度/年の成分)(mAU)
setting0.ini:es2000.txt用設定ファイル[0)]
setting1.ini:es2000.txt用設定ファイル[1)]
setting2.ini:es2000.txt用設定ファイル[2)]
setting3.ini:es2000.txt用設定ファイル[3)]
setting4.ini:es2000.txt用設定ファイル[4)]

===============================================

◎fftcyc.exe(第一段階)
通常の高速フーリエ変換(FFT)は、データ長(2の累乗に限る)の
全体を１周期とする周波数及びその整数倍の周波数しか
計算出来ません。

そこで、データ全体を、開始点をずらしながら数回FFTを行います。
たとえば200,000サンプルのデータにおいて、開始点を10,000点ずつ
ずらしながら、2^17=131,072サンプルのFFTを7回行います。
各FFTでは131,072サンプルを１周期とする周波数及びその整数倍の
周波数のみでフーリエ計算しますが、７回の計算における、
同じ周波数での振幅と位相を比較して変化を見ます。
この比較により、FFTでの周波数分解能より２〜３桁高く
周波数を求める事が出来ます（運が悪くなければ）。

○使用方法
1)まずはデータが入ったファイルに関するデータを入力します。
コメント欄として無視する行の先頭文字、
計算に用いるデータ各行の開始文字位置と文字数、ファイル名です。
「Subtract File」はとりあえず空白にしておいてください。

2)上記を入力後にCheck Dataを押すと、ファイルを１回読み込んで、
データ数と最大値、最小値が表示されます。
必ずしも実行する必要はありませんが、以降で設定する数値を
決定する上での参考値となります。

3)１回のFFTで用いるデータ幅（の２の累乗の指数）を指定します。
当然データ数全体よりも小さい必要があります。
大きい方が個々のFFTの計算精度は良くなりますが、
データ数全体に近すぎるとFFTの回数が取れなくなります。

4)２回目以降のFFT計算をする際の、先頭位置のシフト数を
設定します。
 int{(全データ数−[3)で定めたデータ幅])÷[シフト数]}
だけFFTを行います。
シフト数は、小さければ「極端に大きな振幅の成分がある場合の
誤動作」が防げますし、FFTの回数が増える事で振幅・位相変化の
ばらつきが分かりやすくなりますが、
時間がかかる上に、位相変化が小さくなって出力ファイルの
有効桁数が減ってしまうので、ある程度の大きさは必要でしょう。
このため、全データ数と3)で定めたデータ幅が近過ぎると
FFT回数を確保出来なくなります。

5)角速度を表示する際の単位として、「１周期」に相当する
データ数を入力します。
たとえば１日１回のサンプリングで１年を１周期単位と
したい場合は、この数値を365（または365.25）とします。
この時、１年周期が角速度360deg/cycleとなります。

6)データ先頭から周期の原点に相当する点までの周期単位数を
入力します。
たとえば、周期単位を365.25日、データ先頭が1980年1月1日、
2000年1月1日を原点とする場合は20です。
fftcyc.exe１回目の計算ではそれほど重要な項目ではありませんが、
dftcyc.exeでの位相計算、算出した成分をfftcyc.exeで
差し引く場合に必要となります。

7)FFTではデータ幅[3)で定めた数]の半数と同じ数の周波数について、
フーリエ変換をまとめて行います。131,072サンプルでFFTを行えば、
65,536個の周波数についてフーリエ変換を行います。さらに
数回のFFTを行うため、全部出力するととんでもない量になります。
そこで閾値を指定し、それよりも大きな振幅の成分だけを
取り出します。複数回FFTを行い、指定した閾値を１回でも
超えた周波数だけをファイルに出力します。2)で表示された
最大値、最小値、その差などを参考に決めてください。

8)出力するファイル名を指定します。

9)以上の入力項目は、設定ファイルとして保存できます。
また、後述するdftcyc.exeから読み取り、
必要項目をそのまま使用する事ができます。

10)「Calculate」を押すと計算してファイルを出力します。
FFT自体はそこそこ早いのですが、FFTの回数次第では
意外と時間がかかります。
また、FFT計算のために巨大な配列を用いるので、
同時に起動しているソフトの動作に影響を与える場合があります。

○出力データについて
出力データは次のような数値が、閾値を超えた角速度の数だけ
並びます。

Number of data 45656 
Length of FFT 32768 
Step for FFT-block 3652 
--------------------------------
Avg-Swing(St-D)/D-Swing,Avg-Phase(St-D)/D-Phase:     -0.0014cycle----:     39.9932cycle----:     79.9877cycle----:    119.9822cycle----
>>>>   298(   298.9503deg)-> 299.267170518419 deg
   0.1358(   0.0023), 124.6561deg(  14.1136deg):   0.1324, 105.3560dg:   0.1368, 118.8698dg:   0.1386, 131.0289dg:   0.1355, 143.3695dg
   0.0010(   0.0031),  12.6712deg(   0.6004deg):   0.0000,   0.0000dg:   0.0044,  13.5138dg:   0.0019,  12.1591dg:  -0.0032,  12.3406dg
------------------------
>>>>   314(   315.0014deg)-> 315.557157933504 deg
   0.1172(   0.0036), 262.5140deg(  24.8259deg):   0.1232, 229.2900dg:   0.1165, 251.3999dg:   0.1139, 273.3911dg:   0.1152, 295.9750dg
  -0.0027(   0.0033),  22.2283deg(   0.2561deg):   0.0000,   0.0000dg:  -0.0067,  22.1099dg:  -0.0026,  21.9911dg:   0.0013,  22.5840dg
>>>>   315(   316.0046deg)-> 315.559010845009 deg
   0.1484(   0.0023), 112.9935deg( 141.0746deg):   0.1445,  49.8200dg:   0.1493,  31.6868dg:   0.1505,  14.1059dg:   0.1494, 356.3611dg
   0.0017(   0.0025), -17.8196deg(   0.2316deg):   0.0000,   0.0000dg:   0.0049, -18.1332dg:   0.0012, -17.5809dg:  -0.0011, -17.7448dg
------------------------
>>>>  4816(  4831.3586deg)-> 4832.01960102154 deg
   0.2445(   0.0001), 161.5849deg(  29.5561deg):   0.2444, 121.9321dg:   0.2446, 148.3630dg:   0.2444, 174.8083dg:   0.2444, 201.2363dg
  -0.0000(   0.0001),  26.4347deg(   0.0076deg):   0.0000,   0.0000dg:   0.0001,  26.4308dg:  -0.0001,  26.4453dg:  -0.0000,  26.4280dg
>>>>  4817(  4832.3618deg)-> 4832.01996209175 deg
   0.4727(   0.0001), 281.4432deg(  15.2872deg):   0.4728, 301.9534dg:   0.4726, 288.2808dg:   0.4727, 274.6041dg:   0.4726, 260.9346dg
  -0.0001(   0.0001), -13.6729deg(   0.0030deg):   0.0000,   0.0000dg:  -0.0002, -13.6725dg:   0.0000, -13.6767dg:  -0.0000, -13.6695dg
>>>>  4818(  4833.3650deg)-> 4832.02002352936 deg
   0.1202(   0.0000), 221.2386deg(  60.1421deg):   0.1202, 301.9365dg:   0.1201, 248.1269dg:   0.1202, 194.3322dg:   0.1202, 140.5589dg
   0.0000(   0.0001), -53.7925deg(   0.0149deg):   0.0000,   0.0000dg:  -0.0001, -53.8096dg:   0.0001, -53.7947dg:  -0.0000, -53.7733dg
------------------------

「>>>>」で始まる行は、角速度に関するデータです。
最初の数はFFT計算での角速度の通し番号、
その次の()内がFFT計算に用いた角速度です
(指定した周期単位を１周期＝360度とした場合の数値)。
その次の「->」の後の角速度が
「実際に含まれると思われる角速度の候補」です。
上記の例の上２〜３番目、４〜６番目において、
それぞれ小数点以下２桁目までほぼ一致しているので、
この程度の精度はあると思われます。

が、必ずこの成分があると保証されている訳でもありません。
ありもしない成分が出てくる可能性としては

a)比較的近い２つの周波数が同じ程度の振幅を持っている場合、
その中間の周波数の成分があるように見える事があります。
２つの周波数が近いほど、またはFFTのデータ幅が短いほど
出やすくなります。３つ以上の周波数が近くに集まっている場合は
もっと複雑な現象が起こります。

b)極端に大きな振幅がある場合、やや離れた角速度においても
その影響が出ます。「振幅が大きく、離れた角速度成分」があると、
２つのFFT計算での位相変化が±180度を超えてしまいます。
しかし±180度を超える位相角変化は検出できないので、結果として
「計算間違い」が表示されます（一種のエイリアス現象です）。

逆に、実際にはある周波数成分が現れないケースとして、

c)極端に大きな振幅の成分がある場合に、それに近い周波数で
小さな振幅の成分が隠れてしまう事があります。

b)とc)に関しては、dftcyc.exeで大きな振幅の周波数成分の
振幅と位相を確定し、データから差し引いた後で、
再度操作を行うと、「計算間違い」の部分が消え、
「隠れていた小さな成分」が現れます。
a)は、２つの周波数成分がFFT分解能以上離れていると
検出出来る場合がありますが、近過ぎると識別出来なくなります。
近過ぎる場合は、dftcyc.exeで細かく観察する必要があります。


角速度データの下２行は、フーリエ変換の結果（振幅、位相）
が表示されています。
１行目左から、
振幅の平均値（標準偏差）、位相の平均値（標準偏差）：１回目の振幅、位相：２回目の振幅、位相：…
２行目左から、
振幅変化の平均値（標準偏差）、位相変化の平均値（標準偏差）：0,0:１回目から２回目までの振幅変化、位相変化：…

１行目は振幅以外あまり重要ではありません。振幅も目安程度です。
振幅と位相の正確な値はdftcyc.exeで求めます。
２行目は、振幅変化の大きさや標準偏差が大きい場合、
「変な事が起きている」と疑われます。位相変化の標準偏差が
大きい場合もちょっとあやしいです。
「変な事」が何なのかは、この数値だけからは分からないので、
dftcyc.exeでその周波数を打ち込んで表示してみる必要があります。

=========================================================

◎dftcyc.exe(第一段階)
指定した角速度に対する離散フーリエ変換を、一定データ幅ごとに
計算し、画像表示またはファイル出力します。
また全データに対する離散フーリエ変換の結果も表示します。
ここでいう離散フーリエ変換は、元々の定義に基づく
素直なフーリエ計算です（遅いともいう）。

○使い方
1)「~~~~~~~~」より上は、fftcyc.exeと同様の設定項目です。
fftcyc.exeで保存した設定ファイルを読み込めば完了します。

2)フーリエ計算のための角速度を入力します。
fftcyc.exeの計算結果で出てきた候補などを入れます。

3)１回のフーリエ変換で計算するデータ幅をcycle単位で入力します。
Disp(画像表示用)とFile(ファイル出力用)がありますが、
ファイル出力用の周期で画像表示での色が変わるように
してあります
(したがって画像表示のみでも0以外の値を入れる必要があります)。
短いとそれだけ点数が増えて曲線らしい曲線になりますが、
データ次第では暴れまわります。
長いとそれだけ「平均化」されるので、落ち着いた動きに
なりますが、点数が減ります。
また、データの周期性との兼ね合いで一種のエイリアス現象が
起こる事もあります。数値を変えて何度か試した方がよいでしょう。

4)画像表示の場合は振幅に合わせて倍率を入力します。
画像領域は縦横とも-200〜+200ですので、振幅1程度の場合は
100倍程度が良いでしょう。

5)「Draw」を押せば画像表示を始めて、最後に全データによる
フーリエ計算の結果（振幅と位相）が左側の表示枠に表示されます。

6)出力ファイル名を指定して「Calc.」を押せば、計算結果を
ファイルに出力します。
また、全データによるフーリエ計算の結果を左側の表示枠に
表示し、後述する「Pole data」を黄色の丸で図示します。


○画像の解釈
a)原点以外のほぼ１点に留まっている、１点の周囲を小さく
動き回る程度なら、
「その周波数成分が確かにある。計算結果の振幅と位相を
元にデータから引けば、その周波数成分はほぼ取り除ける」
b)原点の周りを小さく動き回っているなら
「その周波数成分はない」。

c)原点以外の１点の周りで円を描く場合は、
「その周波数成分はある。加えて、その近くの周波数の成分がある。
『近くの周波数』とは、おおよそ
『フーリエ計算に用いている角速度＋円描画の１周期単位当たりの回転角度』」。
d)原点を中心に円を描く場合は
「左周りならばフーリエ計算に用いている周波数は低過ぎ、
右周りならば高過ぎ。修正量は、円描画の１周期単位当たりの回転角度」。

e)直線状に単振動する場合は、
「その周波数成分自体はない。その周波数を中心に挟んで
２つの周波数成分がある。その２つの角速度は
『フーリエ計算に用いている角速度±(単振動の角速度)』」
f)直線、または原点を通るゆるやかな円弧状でほぼ一定速に動く場合、
f1)振幅が１次関数で変化している。
f2)実はe)のケースに相当するが、単振動の周期が長過ぎて
  直線運動のように見える。
f3)（元のデータから近似式を差し引いた後のデータを処理している場合）
  近似式の角速度が実際とは微妙にずれている。そのずれは
  『(近似式の位相と直交する方向への１周期単位当たりの移動距離)÷(近似式の振幅)』(rad/cycle)
g)いわゆる「リサージュ模様」を描く場合は
「２つ（またはそれ以上）の周波数が、計算に用いている周波数の
  そばにある」

d),f)などで振幅や位相の変化を正確に知りたい場合は
ファイル出力を実行して、その内容を確認します。


○ファイル出力
ファイル出力の例です。

Angular Speed: 359.990502 
cycle-----)Swing------,Phase-------:d-Swing----,d-Phase-----:d-Vec(Lng)-,d-Vec(Ang)--
   50.0068)    0.17189,   9.2889deg:-----.-----,----.----deg:-----.-----,----.----deg
  150.0041)    0.49746, 295.0495deg:    0.32558, 279.4850deg:    0.48018, 274.8971deg
  250.0014)    0.08836, 179.2322deg:   -0.40911,-109.5371deg:    0.54182, 123.4911deg
  349.9986)    0.04557,  85.0143deg:   -0.04279, -87.9371deg:    0.10236,  25.5931deg
cycle-----)Swing------,Phase-------:d-Swing----,d-Phase-----:d-Vec(Lng)-,d-Vec(Ang)--
Average---)    0.07398, 175.4160deg:   -0.04211,  27.3369deg:    0.37479, 141.3271deg
          )---Standard Deviation---:    0.29994, 178.5135deg:    0.19427, 102.5564deg
Pole data-)    0.22634, 287.0868deg

左端のcycleは、周期の原点を０として、フーリエ計算の中心点を
周期単位で表示したものです。
Average行での振幅と位相は、全データでフーリエ計算した値です。

中央２列が振幅変化と位相変化です。実際のフーリエ計算は
計算期間に端数を付ける事が出来ませんが、振幅と位相の変化は
指定した計算期間に合わせて補正してあります。
そのため単純な差とは微妙に異なってきます。
上記の例だと、原点付近を横切る形で直線状に変化しているため、
振幅変化と位相変化の平均値・標準偏差が変な値になってます。

右側２列が、振幅・位相の複素平面表示における変化ベクトルを
極座標で表示したものです。中央２列が原点中心に
振幅・位相の変化を測ったものであるのに対して、
こちらは各点からの変化をみたものです。
振幅そのものが変化するケースなどでこちらを利用します。

入力した各速度以外の近い各速度成分がある場合に、
画像描画で見ると弧を描く場合があります。
この弧の中心点を求めた結果が「Pole data」です。
本来は最小自乗法を用いて計算するべきものですが、
手抜き計算なのでかなり誤差があります。
参考程度に使ってください。
具体的に説明すると、
ファイル出力で計算した点(上の例では４点)から重ならないように
２点のペアを２組選び、各ペアの二等分垂線の交点を求めます。
この計算を可能な点の組合せで行い、重み付け平均を出力します。
(上の例では[1,2と2,3],[1,2と2,4],[1,2と3,4],
[1,3と2,4],[1,3と3,4],[2,3]と[3,4]の６通り)


=======================================================
◎計算で得られた成分を差し引いて操作する

fftcyc.exe と dftcyc.exe である程度大きい振幅の成分について
角速度、位相、振幅が分かりました。
小さな成分が大きな成分に隠れている場合もありますし、
大きな成分の振幅が小さく変化している場合もあるので、
算出したデータを元にその成分を差し引き、その上でさらに
詳しく調べます。

まず、差し引くデータを書いたデータ補正ファイルを作成します。
テキストファイルに次のようにデータを並べます。

"deg=3"
359.9905,267.52574deg,6892.3386,-0.1737,0,0
719.9810,265.0823deg,+72.006,0,0,0
#

まず"deg=n"で振幅多項式の次数を決めます。
以降、１行につきひとつの角速度成分です。先頭から
「角速度ω(deg/cyc),位相α(deg),振幅(定数)A,振幅(1次係数)B,振幅(2次係数)C,振幅(3次係数)D…(以下指定した次数まで)」
です。この場合に差し引く値は、周期原点からT(cycle)の点において
(A+B*T+C*T^2+D*T^3+…)*cos(ω*T+α)になります。
実際には、dftcyc.exeから角速度と計算結果をコピー＆ペーストして、
振幅の1次2次係数を計算して求めて埋めるか、分からなければ
とりあえず「,0,0,0,…」で埋めて各行３＋ｎ個にします。
個数を間違えるとエラーで止まります。

最大40行、#または(end of file)で終了です。
#を使えば以降コメントに出来ます。


このようにして作成したデータ補正ファイルを、
fftcyc.exe,dftcyc.exeにおいて「Subtract File」として
指定した上で、それぞれの操作を行うと、
元のデータから各成分を差し引いた値を用いて計算を行います。
fftcyc.exeの「Check」を押せば、最大値最小値は補正後の値が
表示されますし、「Calculate」では各成分を差し引いた後の
データで計算を行い出力します。
したがって差し引かれた角速度の成分は打ち切り誤差程度まで
小さくなります。
（振幅や角速度などが変化する場合は、それに応じた結果が
出力されます。）

以降、他の角速度の成分での計算を行い、振幅と位相が
計算されれば、さらにデータ補正ファイルに書き加え、
最大値最小値をさらに小さくし、さらに小さな成分で
計算をする事が出来ます。

ただしあまり付け加え過ぎると遅くなるので、
適当な所で補正を織り込んで計算し直して出力し、
そのデータファイルで操作した方がよいと思われます。

=======================================================
◎Reference
A.V.Aho, J.E. Hopcroft and J.D. Ullman:The Design and Analysis of Computer Algorithms
        (和訳:アルゴリズムの設計と解析II,野崎昭弘,野下浩平訳,サイエンス社)
=======================================================
H20/07/20:とりあえずリリース
H20/07/22:dftcyc.exeのファイル出力に変化ベクトル(極座標表示)追加
H20/07/26:データ補正機能追加、それに合わせて計算結果出力形式を変更、バグ取り
H21/04/19:dftcyc.exeの色変更,dft幅の端数処理,極座標線表示,
          差し引きファイルの次数設定,
          差し引きファイルの空白行処理,"Pole data"計算
=======================================================
Apr.19th,2009

ftcenter@mth.biglobe.ne.jp
http://www2s.biglobe.ne.jp/~ftcenter/
