◎フーリエ計算ソフト 等間隔の数値列から、その数値列の正弦波成分の周波数、 振幅、位相を算出するためのソフト群です。 ○内容物 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/