Mathematicanの冒険「フーリエ解析入門」

数学に戻る

by H.Yanase

　メニュー
１．フーリエ級数の基礎
２．弦の振動への応用
３．離散フーリエ変換
４．音声解析実験

本講座はMathematicaを利用してフーリエ解析を学び、実際に音声解析をしようとするものである。
Mathematicaを利用する際の注意も記した。
対象は高校卒業程度であるがこだわらない、時として大学程度の数式も登場するが数学や物理に
興味があるもの、手元にMathematicaがあるものは挑戦してほしい。

Mathematicaについて不慣れな人はまず前講座の「Mathematicaの基礎」を一読することを勧める。
まず、このファイルを別な名前で保存しておこう。そうすれば自由に変えたり、失敗もできる。
Mathematicaがない場合でも以下のサイトからMathReaderを無料でダウンロードしてくれば
Mathematicaと同じようにアニメーションや３次元の回転などが可能である。

【はじめに】
パソコンでは音楽や画像などを圧縮したり、再生したりすることができる。これら多くの
情報を処理を解析するのにフーリエ変換の技術が使われる。その基礎をここで見ていこう。

　１．フーリエ級数の基礎

はじめに次のような三角関数の和を描かせてみよう。
xの関数の定義は次のようにする。
ここではｆ１とｆ２の２つの関数を作る。
まずｆ１につてグラフを描かせてみると

三角関数を合成すると実に多様な関数ができる。その基本は次の関係にある。

このように三角関数には周期2πで積分するとサインやコサイン同士はπ、
ことなったもの同士では０になる性質がある。
これを利用するとサインやコサインがごちゃごちゃ混ざった関数からサインやコサイン
がどの程度の強さで混ざっているかを知ることができる。

たとえば次のようなサインが ,コサインが の強さで混ざった関数を異なった
周期で足しあわせたものを定義しよう。

例えば３つ足しあげると次のようなグラフになる

さて、これに次のようにコサイン、サインをかけて積分してみよう。
平均して正しい結果を得るために全体を しておく。

みごとにコサイン、サインの前の係数が出てきた。
足しあわせる数を増やしても、また適当に係数を変えても結果は正しく出てくる。
自分でいろいろ試してみよう。

このようにサインとコサインをかけて積分すればその係数がわかるから、そしたら同じように
係数をつけてサイン、コサインの足しあわせをつくれば元の関数が再現できることになる！

えっ！「元の関数がサインやコサインでなければダメだろう？」それはいい質問だね。
ところが先の足しあわせの関数はｎをいくらでも大きくできるし、周期も変えられる。
するとほとんどの関数はサインやコサインでできてしまうんだよ。

次にこの様子を見ていこう。
例えば次のような凸凹した関数はさすがにサインやコサインでは表現できないだろうと感じるね？

ところが次のようなサインの足しあわせの関数を定義しよう。

もっとｎを大きくしてグラフにしよう。ｎを大きくしても苦にならないのはPCの魅力である。

おみごと！
凸凹関数ができてしまった！。
サインやコサインの係数を積分して求め、次にこの係数でサインやコサインを足しあわせると
どんな関数でも再現できる！？

このサインやコサインを足しあわせた関数をフーリエ級数というよ。
ではフーリエ級数の基本を感じたところで次の練習問題でためして見よう。

【問題１】次の関数を三角関数の和で表現せよ。

【例解】
奇関数（-πから０と０からπの対称性）であるから先のサインを用いる。そして範囲を０からπとして２倍する。
まず、元の関数にサインをかけて係数を求めてみよう。

そして全体を して足しあわせる関数を次のようにつくる。

以下nが3、20、100の場合を図示する。

なんとか再現できたね。
このように元の関数をサインやコサインで分解して適当な係数をかけて
足しあわせたものにすることをフーリエ変換というよ。
でもよく見ると元のグラフの不連続なところ（原点）では
ギザギザがでてくるね。これをギプス現象という。

【問題２】次の関数をフーリエ級数で表せ

【例解】
今度は偶関数（y軸対称）であるからコサインを使う。
さてコサインはｘが０でも１という値をもつから次のようにa0という項
を余分につけることを忘れないように！

係数が求まったから次にコサインの足しあわせを次のように作る。

同様にnが３と10の場合を描く。

元が滑らかで、連続な関数であるためギプス現象が見られない。
つまりフーリエ変換がうまくいく関数はなめらかに連続しているという性質が
あることがわかるね。

　２．弦の振動への応用

フーリエ変換の応用範囲はさらに広い。例えばギターのような弦の振動の解析にも使える。
次のような初期状態から弦をはじくことを最初に考えよう。
例えば弦の長さがここではL=10としｘが４のところで１だけ持ち上げる。

次に前と同じようまず係数を求めよう。上の関数をフーリエサイン変換する。
グラフは２種類あるため場合分けしてその後、加える。

結果を簡単にして整理し、係数Aを決める。

あとは先と同じようにこの係数にサインをかけて足しあわせればよい。
本当は無限大まで足すのだったけどここでは後で実際に計算させるために
次のように最大ｊまで足すようにする。
もう一つ、弦の振動は時間的にも振動するそのために次のように
角振動数をωnとして Cos[ωn t]をかけておく。

このノートブックでは上の図の右タグに下矢印がついてある部分を選択してCTRL＋Yで動く

さらに時間変化も一緒に３次元のグラフにすることもできる。
自由に回転ができるようにするためにパッケージRealTime3Dをロードしておこう。(Ver4以上)

これでマウスで回転して好きな方向から眺めることができる。

外力がある場合（少し高度に）

上の例では実は次のような運動方程式をフーリエ変換を利用して解いている。

さらに次のような外力があるときをここでは考えよう。
フーリエ変換はこのような場合でも簡単に解を出すことができる。
ここでは少々数式を用いることになるので
興味ないものはとばしてよい。

まず一般的に考えよう。弦の長さをｌとする。
外力ｆが働くとき運動方程式は解ｕを次のよう変数分離しておく

外力ｆが次のようにフーリエ変換できたとすると係数Fnは次のように与えられる。

初期条件を次のように与えると

運動方程式からTnについて次が求まる。

この解は次のようになる

従って解は

ここで具体的に先の外力ｆを用いる。
上の一般解に従って係数Fnと解ｕを次のようにおく。
ただし、ｎについては最大ｊまで足すようにしておく。

では適当に数値を代入して図を描かせてみよう。

このノートブックでは上の図の右タグに下矢印がついてある部分を選択してCTRL＋Yで動く

今度は振動が一定ではないのがわかるね上のωnやωの値を変えるとどうなるか試してみるといい。
さらに上の式のAn、Bnはｔ＝０のはじめの状態での弦の様子を決める。
例えばAn=1としてみよう。

最後のにこの様子を３次元にして比較してみよう。

　３．離散フーリエ変換

パソコンなどで扱うデータはデジタル化されている。
こうしたデータを効率よくフーリエ変換するのがＦＦＴと呼ばれる高速フーリエ変換だ。
Mathematicaには任意のデータをFFT処理するためのコマンド、Fourier[データ]が用意されて
いる。ここではこの関数の性質をみていくことにする。
興味ある人は　Fourie　のところにカーソルを持っていってｆ１キーを押してヘルプを見てみよう

なにやらむずかしいことがいっぱい書いてあるようだけどくじけることはない、使いながら理解していこう。

まず離散化したデータをTableコマンドで作成し、ListPlotを利用して作図する。
次の例ではｘが1から200までのデータ列を作成している。
点をつなげない例とPlotJoined->Trueを利用し、つなげた例を以下に示す。

ここで実際に作成したｔ１の中身を覗いてみよう。１０個みるのは次のようにする。

では次のようにこのデータをフーリエ変換し、作図させよう。

データー数200の中点100を中心に左右対称なグラフができた。
グラフのピークを詳しくするために次のようにレンジを変えてみよう。

9のあたりにピークが出ているのがわかるね。
では刻みを２倍にしてデーター数だけ半分にしたらどうなるだろうか。
Mathematicaのマスターのこつは「思いついたらすぐ実行」「困ったときはすぐヘルプ」である。
ではさっそく試行錯誤といこう！

次の例では１から２００まで２ずつ増やしてデータを作っている。

予想通り先と半分の50のところで対称になっている。
でもピークの位置はかわらないね。ただ強さ（ピークの高さ）は小さくなっている。
これだけでデーターの数がピークの高さに関係していることがわかる。

では次の謎解きはピークの位置は何で決まるのかである。
そのためには先の関数f2を使う。
同じようにやってみよう。

今度は１１，２１，６１の３カ所にピークがある。さて、最初は９だった。
これははじめの関数ｆ１，ｆ２をよく見ると出てくる。
データ数は２００、ｆ１の振動数は０．０４かけると８、１を加えると９となる。
ｆ２の場合もデータ数が２００振動数は０．０５，０．１，０．３だから
かけて１を加えると１１，２１，６１でぴったりである。

そう、フーリエ変換のあとのグラフのピークは振動数に対応してる。
こうした振動数の表示をスペクトル表示というよ。
逆にいろいろな情報からその振動数ごとの強さがわかってしまうことになる。
これはすごい！なぜって？
例えば録音された声がだれのものかとか、ある色から赤がどれだけ混ざっているのかとかがわかることに
つながるのだ。

しかし、どうもデータ数をかけたり１をたしたりは面倒である。このあたりをうまくできるように
関数を作ってしまおう。
MathematicaのFourie関数は最初のデータを１ではなく０から数えるので１ずれていたのだ、
これらを解消させるため次のような関数をつくろう。

　SplotとSRplotの２種類をここでは作った。SRplotはｘ軸の範囲を指定するときに利用する。
　上のプログラムでは最大値の半分までをｘの値を１ずらして作図するようにしてある。
　また、データの最大値と数をあらかじめ求めて目盛りが合うように調節した。

この関数はMathematcaのTable関数でデータを作成したときに便利である。次のように使用する。
Splot[リスト名、刻み幅]
SRplot[リスト名、刻み幅、ｘ軸のはじめの値、ｘ軸の終わりの値]
という書式になっている。刻み幅は通常は１でよいが0.1でｘを増やしたいときにはここに１０を入れればよい。
さっそく先のt2のデータを使って試してみよう。

ピークの位置があったし、同じ形が２つも出てこないのですっきりした。
振動数の値もあっている。
でも、FFTではデータの数がデータの間隔（秒数）を決めていることを知っておこう。
これは次のような問題をやるとき重要になる。

【問題3】理解ができたら次の２つの関数をフーリエ変換してパワースペクトルを表示せよ。
    　そして、ピークとなる振動数を[Hz]単位で求めてみよ。

--略解--

最初のは先の結果がそのまま使える。

従って振動数は0.005Hzと0.01Hzだね。

２番目は振動数がとても大きい。従って同じように刻み幅を１にすると問題が生じる。
失敗した人もいるだろう。
繰り返しの周期が出てこないのだ。ＦＦＴを正確にするためにはなるべくたくさんの
周期をデータにしなくてはいけない。データの１間隔の位相差が2π以上になったら意味
をもたない。

そこで次のように刻み幅を細かくしよう。間隔を0.01にして刻み幅に100を代入する。

これで答え５[Hz]と１０[Hz]が見つかったね。
このようにFFTでは重要な定理があることがわかる。
高い振動数の時には刻みを細かくしなくてはいけないということだ。
１秒間にどれだけのデータがとれるかを表すのがサンプリング周波数だ。
ここでの刻みははこの値に一致する。最近のパソコンでは44.1kHzであることが多い。
従ってこのパソコンでFFTを有効にするには最大はその半分の22.05kHz
程度までだということを知っておこう。

FFTの詳しい内容に興味あるものは後記述の参考文献を読んでみるといい。

　４．音声解析実験

では次にデータ解析を実際にしてみよう。
　　このファイルと共にサンプルとして「うなり.wav」を添付したのでこれを利用してもよい。
　　自分のPCにマイクがあって録音が簡単にできるならサウンドレコーダーを開いて自分の
　　声を録音して名前を付けて保存しておこう。

１．録音した音声などの解析

　　サウンドファイルの作成
    マイク・スピーカーのあるウィンドズマシンではサウンドレコーダーで簡単に
    音を取り込むことができる。
    ここではすでに作成した「うなり」というWaveファイルを使おう。
    このファイルをダブルクリックして実際に自分の耳でこの音の様子を聞いておこう。

　　サウンドファイルの取り込み方
    パスを表すのは￥ではなくて半角で/といれる。
    下の例はDドライブにある「うなり」というWAVEファイルを読み込んでいる。
    ファイル名やファイルの場所は異なるので自分の利用する内容に変えよう。
    同じ要領で画像ファイルも取り込むことができる。

　　パッケージの読み込み
　　　Mathematicaは特別な用途のための多くのパッケージを用意している。ここではサウンドに関する次の
　　　パッケージを読み込む

　　次のようにするとファイルの中身の情報を教えてくれる。

　　では次のようにしてサウンドファイルを変数sndに取り込もう。変数名
　　はなんでもよい。

　　「Show」コマンドでこのwaveファイルを図にしてみよう。

　減衰していくうなりの様子が見てとれる、これと先ほど聞いた音と比べてみるといい。
　次のようにするとＭａｔｈｅｍａｔｉｃａからも直接音を鳴らせる

　またＭａｔｈｅｍａｔｉｃａで作った関数は次のようにして鳴らすことができる。
220Hzの振動を１秒間鳴らすには

とする。音楽の得意な人はどれくらいの振動数の差を聞き分けられるか挑戦してみるといい。　
次に技術的だがサウンドファイルからデータを例えば44100個抜き出す関数を次のように定義する。
これで以前にやったリストが44100個できたことになる。

中身を確認するために；をはずしてこのうちのいくつかみてみよう、後にこの間隔が重要になる。

　PlotRangeオプションを変更してフーリエ変換してみよう

　前にみたように真ん中で対称になってるね。これだとあまりに範囲が大きくてよくわからない。
　そこで次のようにレンジを変えよう。

　ピークの一つが３６０付近に出ているね。
　そこでさらにレンジをしぼってみよう。

　３６４のあたりがこの音の振動数が一番大きいことがこれでわかった。ではこれは３６４Ｈｚといって
　いいのだろうか？。
　先にこの音声ファイルの情報をみたときに
       Sampling rate: �`44100
    があった。実はこの数字にデーター数を合わせたから目盛りの値は正しく振動数を示してくれたのだ。
    このサンプリング周波数はとても重要でMathematicaを使わなくてもWindowsでファイルのアイコンを
    右クリックしてプロパティの詳細をみることでもわかるよ。
    
    でも、上の図はMathematicaのFourie関数をそのまま利用したからずれが残っている。
    正しくは次のように先のSRplotを利用し、刻み幅にはサンプリング周波数を入れればよい。

これで正しい振動数がわかった。？
そう、物理を少し知っている君なら不思議に思うだろう。
もともとこの音声ファイルは「うなり」を録音したものだ「うなり」とは
振動数が非常に近い２つの音からできるんだったね。
ところがフーリエ解析したみても振動数のピークは１つしか見えない。
これはおかしい！！。
原因も実は依然の問題演習で経験している。このリストのデータ数はちょうど
44100個、サンプリング周波数も44100Hzだからこれだと１周期分のデータしかない！
そこで多少時間がかかるがこのファイルのデータを次のようにすべて代入する。

今度はMathematicaのFourie関数はつかえないから自作したSRplotを次のように利用する。

みごと２つのピークが見て取れるこの振動数の差がうなりの振動数になって聞こえる。
上の図から差は約２Hzだからおよそ、１秒間に２回うなりが聞こえることになる。
実際と比べてみるとよい。
これらはドップラー効果としてスピードガンや飛行機、潜水艦の速さを知ることにも応用されているよ。

【問題４】同じように自分で作成したファイルのサンプリング周波数を求めてフーリエ解析し、
　　　　　振動数成分を求めよ。
　　　　　　　他の解析ソフトがあれば同じようにＦＦＴ解析した結果と比べてみるとよい。

　少し高度に・・たたみこみ

では次にフ-リエ変換+Mathematica で実験応用例を考えよう.

フーリエ変換には上記のようなたたみこみを満たす関数を利用するとフーリエ変換の積にになるという
性質がある。
これは雑音の入ったデータから目的の情報を取り出す「スムージング処理」に応用できる。
例えば波動でおなじみな次のようなベッセル関数を元にする。

これにランダム関数を使って人工的に雑音を入れる。

たたみこみで利用する次の関数を用意する。

離散フーリエ変換のためのたたみこみの核をずらしておく

ランダム関数Random[]は０から１までの数字をランダムに出すここでは-0.5から+0.5までとした。
前の係数0.1をかえると雑音のレベルが変わる。

　参考文献

    &circlef; Mathematicaの基礎    Mathematicaで何ができるか,主要なコマンドの解説
        Mathmaticaの基礎    http://www.ne.jp/asahi/buturiwa/ai/educationM.htm
        Mathematicaハンドブック　秀和システム    木村　広著　
               
    &circlef; フ-リエ変換などの数学的基礎
        フーリエ解析　　　    岩波書店　　　    福田礼次郎著
        偏微分方程式    　    岩波書店    　　　及川正行著
        
    &circlef;　FFTやプログラム関係
        Cによる科学技術計算    CQ出版社        小池真一著
  
  ◎　Mathematica公式サイト:    http://www.wolfram.com
  
  ２００１．２        H.Yanase　http://www.ne.jp/asahi/buturiwa/ai/educationM.htm