フーリエ変換で周期を発見 | 科学アラカルト [実践] | 逆ポーランド電卓の実践ウェブ rpn hacks!

逆ポーランド電卓の実践ウェブ rpn hacks!

逆ポーランド電卓rpnの実践ウェブ   
rpn hacks! アールピーエヌ・ハックスサイトマップ

rpn | 実戦 | 科学アラカルト | 科学に関する話題をrpnで探求!調べて実験してみよう。

HOME > 実践 > 科学アラカルト > フーリエ変換で周期を発見

hatena twitter facebook rss ソーシャルブックマーク

フーリエ変換で周期を発見

 波のように寄せては返したり、人口が密になったり疎になったり、経済が強くなったり弱くなったり、勢力が広がったり縮こまったり…。自然現象から人の営みに至るまで、世の中に周期的に増減を繰り返すものは多くあります。

「待てば海路の日和あり」の諺のように、天気だって良くなったり悪くなったりで、ずっと悪天候のままってことはありません。人間の経済活動も似たものがあります。景気は良くなったり、悪くなったりします。永遠に好景気はありませんし、永遠の不況もありません。大抵は数年から十年単位で好況と不況を繰り返しています。

生物界でもイエローストーン国立公園で増えすぎた鹿と新たに導入された狼の関係は有名です。現在はお互い増えれば減り、減れば増えながら均衡を保っています。鹿だけがいる環境では大繁殖して生態系を崩しますが、狼がいることで公園内の生物ピラミッドが適正に保たれています。

増減を繰り返す波の性質

 このように増減を繰り返しながらも元の状態に戻ってくるようなものには、波の性質があると言います。増減の量を縦軸に時間経過を横軸にすると、こんなイメージですね。

  ^ 増                                     *                  
  *    *                             *                        
  |                                              *            
  |          *                                                
  +----------------------------*-----------------------------> 時間
  |                                                    *      
  |                *                                          
  | 減                   *                                   *


縦軸の増減は重さかもしれないし、電流量かもしれませんし、ある基準からの距離かもしれませんが、時系列に並べてみると同じような増減を繰り返すのが特徴です。そして、この一巡りが周期になります。

例えば、音も波の性質を持っています。空気が密の場所と疎の場所が伝播しながら空間を広がっていきます。なお、音は進行方向に対して縦の変動になるので縦波です。対して、電磁気は横波になります(彼のテスラは縦波を発見していたとの噂もあります)。

さて、知りたいのは、その周期です。一体どれくらいで一巡するのか。それが分かれば(予測できれば)何かに利用できることがあるかもしれません。

複数の要因が絡む現実の波

 増減量を時系列にグラフにして波が目視できれば周期もなんとなく分かります。でも、複数の要因が絡んでいるとグラフ化してもそう簡単に周期が見えてきません。

ここにアイスクリーム屋の販売額の推移データがあるとします。販売に気温が影響しやすいことは分かりますが、湿度も関係あるかもしれませんし、近くに他社のアイスクリーム屋が出没していることが影響するかもしれません。

仮に気温や湿度にも周期があり、ライバルの屋台の出現にも周期(曜日単位で変化)があるとすると、3つの要因が絡み合うことになります。丁度、次のグラフのようにです。

                           (a)                      (b)
  **                       ***                      ***
  * *   **                       **                       **
  |  * *              *   *   * *  *           *   *   * *  *
  |   *   *          * *       *              * *       *
  +------------*--------**--------------*--------**----------*
  |        *  * *   *               *  * *   *
  |          *                       **   *                  *
  |         *    * *                        *


a点とb点を周期とする波はぼんやりと見えてきますが、複雑そうです。他に波が混じっているのかどうか明確には分かりません。実際、この例では3つの波が絡み合っています。このように、グラフによる目視で周期を見つけることは、限界があることがすぐに分かります。

でも、これらの周期を瞬時に数値的に解析することができれば、強力な武器となることは容易に想像がつきます。

フーリエ解析・フーリエ変換

 200年以上前、フランスの数学者・物理学者だったフーリエ(1768~1830)は周期的に繰り返すどんな波も三角関数(sinとcos)の合成で置き変えることができることを証明しました。

鐘やピアノの音など、複数の波が複雑に絡み合っていたとしても、周期的な量の増減を1つ1つ紐解くように、全て三角関数の合成として表せるという画期的な発明です。

フーリエが生きていた頃、日本は江戸時代の後期(化政文化)でした。時は移って明治の開国。輸入された洋算に対し、日本の和算が引けをとっていなかったのは事実(微分も積分も既にあった)なのですが、そのまま和算を発展させてもフーリエの発想までは届かなかった気がします。

このフーリエ解析を通せば、単調な音から複雑な音、電圧や電流の増減、経済活動や自然現象に至るまで周期が含まれているものを全て解析することができます。

そして、知りたかった変動の周期まで数値的に炙り出すことができるわけです。まずはフーリエ解析の力の一端を示す例として、時報の音を解析してみましょう。

時報の音をフーリエ解析

 デジタル化によって、最近は0分(O'clock)3秒前の「ポッ・ポッ・ポッ・ポーン」という時報を聞く機会が減りました。従来のアナログ伝送と違ってデジタル伝送には遅延が生じるため、時刻ピッタリに音を鳴らすことができなくなったためです。

時報「ポッ・ポッ・ポッ・ポーン」のメロディーは、440Hz、440Hz、440Hz、880Hzの音の並びです。ここでは、時報のサウンドファイルから「ポッ」の部分だけを切り出して、データ化することにしてみました。データはファイルの"440hz.txt"にあるものとします。内容は以下のとおりです。

  >type 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の概要を見てみましょう。まずはデータの概要を知るための基本統計量です。

  >rpn -c statinfo <440hz.txt
  デ ー タ 412
  最 小 値 -5859
  最 大 値 5805
  範    囲 11664
  合 計 値(Σ) -1671
  平 均 値(μ) -4.05583
  分 散 値(σ2) 1.56845e+07
  標準偏差(σ) 3960.37
  分 散 値(s2) 1.57227e+07
  標準偏差(s) 3965.18
  歪度(a3≒0) 0.000137656
  尖度(a4≒3) 1.57202
  変動係数(ν) -977.651


サンプルデータが412個、最小が-5859で最大が5805のようです。幹葉表示にしてみます。

  >rpn -c stemleaf <440hz.txt
  -5 | 000000001122222233333444444555555566666666677777777777777788888
  -4 | 011112223334444445666666668888888899
  -3 | 00001113333444566666778899999
  -2 | 0001111233344444557777788
  -1 | 00000112333444455777778899
  -0 | 0000122333344466677778899
   0 | 0000000011233333344466667777778
   1 | 0000113334445667778889
   2 | 00001112333444557777788899
   3 | 0000011122333334456666677888999999
   4 | 011111223344455666667788888899
   5 | 00000011111222233333334444444555555555666666666677777777777777778


波の性質を持っている時系列データなので当然ですが、ガウス分布でも指数分布でもない様子が見受けられます。時系列にしてグラフ化してみると次のようになります。

  >rpn 1 -c rownum <440hz.txt | xyp -x,412 -y-6000,6000 -m -w60
  ^y**0*0   ** **** *     * **** **      **** ***      ***
  | * **** ***     **** ****     *** **** *    ** ******
  |           ****           ****          ****         ****
  | ** **  ***    * **   ** *    ***   * **    ***    ***  *
  |*    ***           ****          ****          ***      *
  |** **    **     ***    **     **     * *    * *   * **
  |*          * **           ** *          ****          *
  |***        ** ***        ** ** *      *** * *       ** ****
  *---*-*--**--------**-**-*-------**-*-*-------**-*-**------>
  *o   ****          * * *          ****          ****     **x
  |* **     ***    **     **     **     **     * *    * * **
  | *          ****         *****          ****          ***
  |  **     * *   ** *   * **   * **    ***    ***    ** *
  |*    ****          ***           ****          ****    *
  | ****   * ********     ******** *    ******* **    * ***
  |-6*** *** **    ** *** **     * *** ***      *** ****


プロットが詰まり過ぎていてよく分かりませんね。最初から30サンプリングまでの部分を拡大表示してみましょう。

  >rpn 1 -c rownum <440hz.txt | xyp -x,30 -y-6000,6000 -m -w60
  ^y 6000                            *                 *
  |                                    *                 *
  |
  |                                *                 *
  |                * *
  |              *                       *                 *
  |                    *
  |            *                 *                 *        30
  +*-*-*-*--------------------------------------------------->
  |o       * *                                               x
  |                      *                 *                 *
  |                            *
  |                                              *
  |                        *
  |                          *               *
  |-6000                                       *


キレイなカーブが見えてきました。そして、徐々に音の振幅が大きくなっているのが分かります。一定の振幅を繰り返した後、最後のほうでは徐々に振幅が小さくなっているはずです。

このサウンドデータをフーリエ変換して、周期(周波数)を見つけ出します。うまくいけば、440Hzが導かれるはずです。

離散フーリエ変換を試す

 フーリエ変換を計算するには、離散フーリエ変換:DFT(Discrete Fourier Transform)を行ないますが、処理手順が複雑なのでrpnプログラムを用意しました。プログラム名はfourierです。

高速フーリエ変換:FFTもありますが、サンプル数が少ないのでDFTで十分ですし、2の冪乗個のサンプル数という制約もないので便利です。ただし、1000サンプル程度までとなります。

プログラムはn行1列の時系列データを受け取って、DFT計算後、スペクトルまで一気に出力します。

  時系列データ --> fourierプログラム --> スペクトル


さっそく、サンプリングしたデータファイルから計算してみましょう。rpnを使ったフーリエ変換の計算はたったの1行で済みます。

  >rpn -c fourier <440hz.txt >spec440.txt


計算結果はspec440.txtに格納されます。では、spec440.txtに行番号を付けて、計算結果を表示してみましょう。

  >rpn 0 -c rownum <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個計算されています。一応、基本統計量も出しておきましょう。

  >rpn -c statinfo <spec440.txt
  デ ー タ        412
  最 小 値        7
  最 大 値        974491
  範    囲        974484
  合 計 値(Σ)    5.8688e+06
  平 均 値(μ)    14244.7
  分 散 値(σ2)   6.25911e+09
  標準偏差(σ)    79114.6
  分 散 値(s2)   6.27434e+09
  標準偏差(s)    79210.7
  歪度(a3≒0)    9.86647
  尖度(a4≒3)    111.741
  変動係数(ν)    5.56073


さて、グラフ化したスペクトルは次のとおりです。計算結果の最初の数値は無視します。また、スペクトル表示の2つに分けた後半は無視する決まりです。

  >rpn 0 -c rownum <spec440.txt | xyp -x,411 -y,1200000 -s100 -n -m
  ^y 1200000
  |
  |   *                             *
  |
  |
  |
  |
  |
  |   *                             *
  |
  |
  |
  |   *                             *
  |   *                             **   x
  ****************************************
  +--------100-------200-------300------4>


ピークが2つ見えますね。412÷2の206のあたりがナイキスト周波数になるので、そこから前半のピークが有効で後半のそれは無視することになります。そこで、前半の部分を拡大表示してみましょう。

  >rpn 0 -c rownum <spec440.txt | rpn 1 205 -c lookup >tmp
  >xyp -x,205 -y,120000 -s50 -n -m <tmp
  ^y 120000
  |
  |       **
  |
  |
  |        *
  |
  |        *
  |      * *
  |
  |         *
  |      *  *
  |         *
  |   ****  **                           x
  *****  *   *****************************
  +--------50--------100-------150------->


どうやらピークが横軸メモリの50弱のところにありそうです。今回のサウンドファイルはサンプリングレートが4kHzなので4000回/秒です。実際のサンプル数は412だったので、1サンプル分は約9.7Hzになります。

  >rpn 4000 412 /
  9.70874


次に最大値が974491であることが既に分かっているので、その値が何番目にあたるか調べます。

  >rpn 0 -c rownum <spec440.txt | rpn x 974491 . -c lookup
  974491 45
  974491 367


45番目と367番目です。スペクトルの後半は無視しますから、45番目が最大値になります。そこで、先ほどの約9.7Hzに45を掛けるとスペクトルがピークを示す周波数が分かります。

  >rpn 9.70874 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)を使って合成されているので、人工的な音でもあり、フーリエ解析できれいに周期を取り出せるのは当然かもしれません。現実の世界では複数の要因が影響しあって、複雑な波を描くので簡単に分析できない可能性もあります。

そこで、太陽の黒点数と株価の時系列データを題材にして、フーリエ変換の実際への応用を考えてみましょう。

続き(part2)はこちらLinkIcon

警告rpnプログラムを実行するには、rpn試用版かrpn標準版が必要です(バージョンの違いはこちら)。

情報pasteは講座サポートで公開されています。rownumはカンタン分析パッケージに同梱されています。lookupはユーティリティーパッケージに同梱されています。statinfo, stemleafはrpnマイスターパッケージに同梱されています。xypとnpdはrpnの姉妹ソフトウェアです。詳しくはプロダクトを参照ください。

科学アラカルト

調べて実験してみよう

※実践コーナーのTOP

紹介 rpnの利用シーンはこちら…

講座初めての人のrpn基礎もどうぞ
講座しっかり学べるrpn入門もどうぞ
講座すぐに使えるdos入門もどうぞ
物語データを見抜くojt物語もどうぞ

実践他の分野への挑戦は実践TOP

応用rpnアプリケーションは応用TOP

part2

調べて実験してみよう

※実践コーナーのTOP

講座初めての人のrpn基礎もどうぞ
講座しっかり学べるrpn入門もどうぞ
講座すぐに使えるdos入門もどうぞ
物語データを見抜くojt物語もどうぞ

実践他の分野への挑戦は実践TOP

応用rpnアプリケーションは応用TOP

書籍紹介

記事に関連した書籍

本ウェブサイトで扱った話題に関連した書物で、スタッフが実際に読了したものを紹介。

書籍数学の書籍
数の世界は思ったよりもエキサイティング

書籍投資の書籍
失敗しない投資には広範囲で実践的な知識が必要

警告バックスラッシュはエンマークに読み替えて下さい(IE)。
バックスラッシュとエンマーク

警告文字で作られた図表や式が崩れることがあります。ブラウザによっては固定幅フォントをMSゴシックにするときれいに表示されます。それでも崩れる場合は図表や式をメモ帳にコピー後、閲覧下さい。

警告rpn試用版と標準版(2kリビジョン)はダブルクォートで囲って下さい。

rpn 1 2 + ⇒ rpn "1 2 +"
rpn 1 -c foo ⇒ rpn "1" -c "foo"

ダブルクォートは省略できることが多いのですが、慣れない間は囲んだほうが無難です。なお、本ウェブサイトの記事ではrpn標準版(98リビジョン)を使用しているため囲っていません。詳しくは技術サポートの「rpn TIPS参照ください。

注意rpnの障害情報はこちら

警告rpn試用版の場合、複雑なプログラムや処理時間のかかるプログラムの一部には動作しないものがあるかもしれません。あくまで無料提供であることを勘案・了承ください。rpn標準版は、すべてのプログラムが動作します。