フーリエ変換で周期を発見
波のように寄せては返したり、人口が密になったり疎になったり、経済が強くなったり弱くなったり、勢力が広がったり縮こまったり…。自然現象から人の営みに至るまで、世の中に周期的に増減を繰り返すものは多くあります。
「待てば海路の日和あり」の諺のように、天気だって良くなったり悪くなったりで、ずっと悪天候のままってことはありません。人間の経済活動も似たものがあります。景気は良くなったり、悪くなったりします。永遠に好景気はありませんし、永遠の不況もありません。大抵は数年から十年単位で好況と不況を繰り返しています。
生物界でもイエローストーン国立公園で増えすぎた鹿と新たに導入された狼の関係は有名です。現在はお互い増えれば減り、減れば増えながら均衡を保っています。鹿だけがいる環境では大繁殖して生態系を崩しますが、狼がいることで公園内の生物ピラミッドが適正に保たれています。
増減を繰り返す波の性質
このように増減を繰り返しながらも元の状態に戻ってくるようなものには、波の性質があると言います。増減の量を縦軸に時間経過を横軸にすると、こんなイメージですね。
* * *
| *
| *
+----------------------------*-----------------------------> 時間
| *
| *
| 減 * *
縦軸の増減は重さかもしれないし、電流量かもしれませんし、ある基準からの距離かもしれませんが、時系列に並べてみると同じような増減を繰り返すのが特徴です。そして、この一巡りが周期になります。
例えば、音も波の性質を持っています。空気が密の場所と疎の場所が伝播しながら空間を広がっていきます。なお、音は進行方向に対して縦の変動になるので縦波です。対して、電磁気は横波になります(彼のテスラは縦波を発見していたとの噂もあります)。
さて、知りたいのは、その周期です。一体どれくらいで一巡するのか。それが分かれば(予測できれば)何かに利用できることがあるかもしれません。
複数の要因が絡む現実の波
増減量を時系列にグラフにして波が目視できれば周期もなんとなく分かります。でも、複数の要因が絡んでいるとグラフ化してもそう簡単に周期が見えてきません。
ここにアイスクリーム屋の販売額の推移データがあるとします。販売に気温が影響しやすいことは分かりますが、湿度も関係あるかもしれませんし、近くに他社のアイスクリーム屋が出没していることが影響するかもしれません。
仮に気温や湿度にも周期があり、ライバルの屋台の出現にも周期(曜日単位で変化)があるとすると、3つの要因が絡み合うことになります。丁度、次のグラフのようにです。
** *** ***
* * ** ** **
| * * * * * * * * * * * *
| * * * * * * * *
+------------*--------**--------------*--------**----------*
| * * * * * * * *
| * ** * *
| * * * *
a点とb点を周期とする波はぼんやりと見えてきますが、複雑そうです。他に波が混じっているのかどうか明確には分かりません。実際、この例では3つの波が絡み合っています。このように、グラフによる目視で周期を見つけることは、限界があることがすぐに分かります。
でも、これらの周期を瞬時に数値的に解析することができれば、強力な武器となることは容易に想像がつきます。
フーリエ解析・フーリエ変換
200年以上前、フランスの数学者・物理学者だったフーリエ(1768～1830)は周期的に繰り返すどんな波も三角関数(sinとcos)の合成で置き変えることができることを証明しました。
鐘やピアノの音など、複数の波が複雑に絡み合っていたとしても、周期的な量の増減を1つ1つ紐解くように、全て三角関数の合成として表せるという画期的な発明です。
フーリエが生きていた頃、日本は江戸時代の後期(化政文化)でした。時は移って明治の開国。輸入された洋算に対し、日本の和算が引けをとっていなかったのは事実(微分も積分も既にあった)なのですが、そのまま和算を発展させてもフーリエの発想までは届かなかった気がします。
このフーリエ解析を通せば、単調な音から複雑な音、電圧や電流の増減、経済活動や自然現象に至るまで周期が含まれているものを全て解析することができます。
そして、知りたかった変動の周期まで数値的に炙り出すことができるわけです。まずはフーリエ解析の力の一端を示す例として、時報の音を解析してみましょう。
時報の音をフーリエ解析
デジタル化によって、最近は0分(O'clock)3秒前の「ポッ・ポッ・ポッ・ポーン」という時報を聞く機会が減りました。従来のアナログ伝送と違ってデジタル伝送には遅延が生じるため、時刻ピッタリに音を鳴らすことができなくなったためです。
時報「ポッ・ポッ・ポッ・ポーン」のメロディーは、440Hz、440Hz、440Hz、880Hzの音の並びです。ここでは、時報のサウンドファイルから「ポッ」の部分だけを切り出して、データ化することにしてみました。データはファイルの"440hz.txt"にあるものとします。内容は以下のとおりです。
0
0
-2
-303
-908
:
(中略)
:
600
339
2
1
1
フーリエ変換する際の注意としては、解析する波の最大周波数の2倍以上でサンプリングしなければならないことがあります。つまり、440Hzの2倍以上のサンプリングレート(880Hz以上)が必要です。
人間が耳で聞き取れる音はおおよそ20Hzから20000Hzです。低い周波数は聞き取りにくく、高い周波数は加齢と共に聞こえなくなっていきます。なお、CDは44.1kHzでサンプリングしてあるので、人が聞こえる周波数をカバーしていることになります。
音の高低(周波数)とは別に音の大きさ(音圧)の単位にdB(デシベル)があります。dは10の意のデシでBは電話を発明したベルに因んでいて、dBは音圧を表す常用対数の無次元単位。人間は音の大きさを絶対値ではなく、対数に比例した感覚で捉えています(ウェーバ・フェヒナの法則)。
人間が聞ける最も小さな音は0dB(雪の降る音)で10dB増える毎に10倍になります。例えば、会話の60dBは囁き声の30dBに比べて30dB大きいので10の3乗(1000倍)の大きさです。大きな音には多数の車が通り過ぎる音(70dB)、工場の騒音(90dB)、ジェットエンジン(120dB)、銃声(140dB)などがありますが、持続的な 90dB以上の音は難聴になる可能性があるので注意が必要です。
今回は4000Hzでサンプリングすることにします。この場合、半分の周波数2000Hzをナイキスト周波数といいますが、フーリエ変換ではナイキスト周波数以上の高い周波数は解析できません。なお、量子化は16ビットを使用し、-32768～+32767とします。
さっそく、440hz.txtの概要を見てみましょう。まずはデータの概要を知るための基本統計量です。
デ ー タ 412
最 小 値 -5859
最 大 値 5805
範 囲 11664
合 計 値(Σ) -1671
平 均 値(μ) -4.05583
分 散 値(σ2) 1.56845e+07
標準偏差(σ) 3960.37
分 散 値(ｓ2) 1.57227e+07
標準偏差(ｓ) 3965.18
歪度(ａ3≒0) 0.000137656
尖度(ａ4≒3) 1.57202
変動係数(ν) -977.651
サンプルデータが412個、最小が-5859で最大が5805のようです。幹葉表示にしてみます。
-5 | 000000001122222233333444444555555566666666677777777777777788888
-4 | 011112223334444445666666668888888899
-3 | 00001113333444566666778899999
-2 | 0001111233344444557777788
-1 | 00000112333444455777778899
-0 | 0000122333344466677778899
0 | 0000000011233333344466667777778
1 | 0000113334445667778889
2 | 00001112333444557777788899
3 | 0000011122333334456666677888999999
4 | 011111223344455666667788888899
5 | 00000011111222233333334444444555555555666666666677777777777777778
波の性質を持っている時系列データなので当然ですが、ガウス分布でも指数分布でもない様子が見受けられます。時系列にしてグラフ化してみると次のようになります。
^y**0*0 ** **** * * **** ** **** *** ***
| * **** *** **** **** *** **** * ** ******
| **** **** **** ****
| ** ** *** * ** ** * *** * ** *** *** *
|* *** **** **** *** *
|** ** ** *** ** ** * * * * * **
|* * ** ** * **** *
|*** ** *** ** ** * *** * * ** ****
*---*-*--**--------**-**-*-------**-*-*-------**-*-**------>
*o **** * * * **** **** **x
|* ** *** ** ** ** ** * * * * **
| * **** ***** **** ***
| ** * * ** * * ** * ** *** *** ** *
|* **** *** **** **** *
| **** * ******** ******** * ******* ** * ***
|-6*** *** ** ** *** ** * *** *** *** ****
プロットが詰まり過ぎていてよく分かりませんね。最初から30サンプリングまでの部分を拡大表示してみましょう。
^y 6000 * *
| * *
|
| * *
| * *
| * * *
| *
| * * * 30
+*-*-*-*--------------------------------------------------->
|o * * x
| * * *
| *
| *
| *
| * *
|-6000 *
キレイなカーブが見えてきました。そして、徐々に音の振幅が大きくなっているのが分かります。一定の振幅を繰り返した後、最後のほうでは徐々に振幅が小さくなっているはずです。
このサウンドデータをフーリエ変換して、周期(周波数)を見つけ出します。うまくいけば、440Hzが導かれるはずです。
離散フーリエ変換を試す
フーリエ変換を計算するには、離散フーリエ変換:DFT(Discrete Fourier Transform)を行ないますが、処理手順が複雑なのでrpnプログラムを用意しました。プログラム名はfourierです。
高速フーリエ変換:FFTもありますが、サンプル数が少ないのでDFTで十分ですし、2の冪乗個のサンプル数という制約もないので便利です。ただし、1000サンプル程度までとなります。
プログラムはn行1列の時系列データを受け取って、DFT計算後、スペクトルまで一気に出力します。
さっそく、サンプリングしたデータファイルから計算してみましょう。rpnを使ったフーリエ変換の計算はたったの1行で済みます。
計算結果はspec440.txtに格納されます。では、spec440.txtに行番号を付けて、計算結果を表示してみましょう。
0 1671
1 877.696
2 354.916
3 220.734
4 402.434
:
(中略)
:
407 974.978
408 402.434
409 220.734
410 354.916
411 877.696
スペクトルの振幅が412個計算されています。一応、基本統計量も出しておきましょう。
デ ー タ 412
最 小 値 7
最 大 値 974491
範 囲 974484
合 計 値(Σ) 5.8688e+06
平 均 値(μ) 14244.7
分 散 値(σ2) 6.25911e+09
標準偏差(σ) 79114.6
分 散 値(ｓ2) 6.27434e+09
標準偏差(ｓ) 79210.7
歪度(ａ3≒0) 9.86647
尖度(ａ4≒3) 111.741
変動係数(ν) 5.56073
さて、グラフ化したスペクトルは次のとおりです。計算結果の最初の数値は無視します。また、スペクトル表示の2つに分けた後半は無視する決まりです。
^y 1200000
|
| * *
|
|
|
|
|
| * *
|
|
|
| * *
| * ** x
****************************************
+--------100-------200-------300------4>
ピークが2つ見えますね。412÷2の206のあたりがナイキスト周波数になるので、そこから前半のピークが有効で後半のそれは無視することになります。そこで、前半の部分を拡大表示してみましょう。
>xyp -x,205 -y,120000 -s50 -n -m <tmp
^y 120000
|
| **
|
|
| *
|
| *
| * *
|
| *
| * *
| *
| **** ** x
***** * *****************************
+--------50--------100-------150------->
どうやらピークが横軸メモリの50弱のところにありそうです。今回のサウンドファイルはサンプリングレートが4kHzなので4000回/秒です。実際のサンプル数は412だったので、1サンプル分は約9.7Hzになります。
9.70874
次に最大値が974491であることが既に分かっているので、その値が何番目にあたるか調べます。
974491 45
974491 367
45番目と367番目です。スペクトルの後半は無視しますから、45番目が最大値になります。そこで、先ほどの約9.7Hzに45を掛けるとスペクトルがピークを示す周波数が分かります。
436.893
約437Hzです。理論値の440Hzと0.7%の誤差です。ほとんど一致しましたね。
メスの蚊が飛ぶときの羽音は440Hz付近(A=ラの音)という話があります。あの嫌なプーンという音はA(ラ)の音だったんですね。その音を聞いてオスの蚊がメスを探すそうです。
赤ちゃんの泣き声も440Hzとのこと。絶対音感を持つ人が実際に確認したという話もあるのでほぼ確実でしょう。赤ちゃんが要求を伝えるには泣くしかありません。きっと440Hzは、それを聞くものにとって放ってはおけない気分にさせる周波数なのでしょう。
そして、電話の受話器を上げたときのツーというダイヤルトーン:DTの音はG(ソ)です。勧告によるとダイヤルトーンの周波数は約400Hzなので、G=392Hzとはちょっと違いますが、実際にチューナーで確認するとGより少し高いくらい(GとG#の間)なので、強ち間違ってはいません。ﾌﾚｯﾄ 0--1--2--3--4--5--6--7--8--9--10-11-12 例:ギターの弦が発する周波数(0ﾌﾚｯﾄ)
1弦 E--F--F#-G--G#-A--A#-B--C--C#-D--D#-E 1弦:330Hz 2弦:246Hz 3弦:196Hz
2弦 B--C--C#-D--D#-E--F--F#-G--G#-A--A#-B 4弦:147Hz 5弦:110Hz 6弦: 82Hz
3弦 G--G#-A--A#-B--C--C#-D--D#-E--F--F#-G ベースはギターの3弦～6弦の1/2Hz
4弦 D--D#-E--F--F#-G--G#-A--A#-B--C--C#-D
5弦 A--A#-B--C--C#-D--D#-E--F--F#-G--G#-A ※電話のDTは1弦3ﾌﾚｯﾄを弾いた音(G)に
6弦 E--F--F#-G--G#-A--A#-B--C--C#-D--D#-E 等しい。440Hz(A)の音は1弦5ﾌﾚｯﾄ。
6弦ギターで出せる周波数が 82Hz～1320Hzなのに対して、88鍵ピアノは27.5Hz～4186Hzもあります。なお、音楽CDのサンプリング周波数は44.1kHzなので、上記の原則から22.05kHzまで再現可能なのですが、残念ながら人間が音階として聞き取れる範囲は20Hz～4000Hzくらいです。
次は…
フーリエ変換は自然界の周期的な変動を捉えて、それぞれシンプルな波に分解することができます。そこで、音階のA(ラ)でもある440Hzの音を実際にサンプリングしたデータを使って、フーリエ変換してみました。結果は誤差がほとんどない満足のいくものです。
ただし、時報音は三角関数の正弦波(sin)を使って合成されているので、人工的な音でもあり、フーリエ解析できれいに周期を取り出せるのは当然かもしれません。現実の世界では複数の要因が影響しあって、複雑な波を描くので簡単に分析できない可能性もあります。
そこで、太陽の黒点数と株価の時系列データを題材にして、フーリエ変換の実際への応用を考えてみましょう。
rpnプログラムを実行するには、rpn試用版かrpn標準版が必要です(バージョンの違いはこちら)。
pasteは講座サポートで公開されています。rownumはカンタン分析パッケージに同梱されています。lookupはユーティリティーパッケージに同梱されています。statinfo, stemleafはrpnマイスターパッケージに同梱されています。xypとnpdはrpnの姉妹ソフトウェアです。詳しくはプロダクトを参照ください。
/