パソコンの画面を使って、1限の講義の補足を少々
標本化定理・エイリアシング
講義でも述べたように、ある周期で変動する信号を測定するにあたっては、その 周期の半分より短い時間間隔 で測定を行わなければならない。 これが「標本化定理 (サンプリング定理)」であった。 また測定の時間間隔が標本化定理を満たしていない場合には「エイリアシング (aliasing)」が起こり、本来の周期とは異なる 長い周期の信号として測定されてしまう のであった。 ここでは「エイリアシング」が起こる様子を、Gnuplot の画面を使って確かめてみよう。
端末から Gnuplot を起動し、以下のように入力してみよ (gnuplot> は Gnuplot のプロンプト、コピペ可)。
gnuplot> plot [0:120] sin(2*pi*x/12),"<awk 'BEGIN{for (x=0;x<=120;x=x+3) print x,sin(2*x*3.141592/12)}'" with linespointsこれにより赤線で \(y=\sin(\frac{2\pi}{12}{x})\) のグラフ、緑の点と線でこの \(y\) の値を \(x=0,3,6,9,12,\ldots\) と 3 おきに取り出した値と、それらを結んだ折れ線が描かれるはずだ。 なお、もし赤線のグラフが余りにも「カクカク」している場合には、
gnuplot> set samples 1200などとしてから上の操作を再び行うとよい。 この場合は測定の時間間隔 (= 3) が信号の周期 (= 12) よりも十分短く、標本化定理を満足しているので、緑のように「飛び飛び」の点のデータであっても元の信号 (赤) をうまく再現することができている。 では、測定の時間間隔が長くなり、標本化定理を満足しなくなったら、どのような信号として測定されてしまうことになるだろうか。
上の例にならい、測定の時間間隔を 3 からどんどん大きく していった場合 (ただし話を簡単にするため、時間間隔は整数で) に、測定される信号 (緑) がどの程度うまく元の信号 (赤) を再現できるかを確認せよ。 特に 測定の時間間隔が標本化定理を満足しない場合 (例えば 9、11、13、15 とか) に、どのような信号が「見かけ上」得られてしまうのか を確認せよ (これがエイリアシング)。
自己相関係数
これまでこの演習で用いている「松山市の月平均気温」のデータには、当然予想されるように1年 (=12ヶ月) を単位とする規則的な変動が強く表われている。 では、まだ時間変動の規則性が分からない時系列データの周期を求めるにはどうすればよいだろうか?
時間的 (あるいは空間的でも可) な変動の周期を見積る上で理想的な方法は、もちろん「フーリエ (Fourier) 変換」を用いることなのだが、 以下では、「フーリエ変換」の準備ともなる「自己相関」という概念について紹介する。
ある時系列データに対し、ある時間幅だけずらしたデータと元のデータとの関係を考えてみる。 端末から Gnuplot を起動し、以下のように入力してみよう (gnuplot> は Gnuplot のプロンプト、これもコピペ可)。
gnuplot> plot [1488:1608] "matsuyama-kion2.txt" using 1:3 with lines, "matsuyama-kion2.txt" using ($1+2):3 with linesこれにより、最近10年 (=120ヶ月) 分の毎月の平均気温が赤の折れ線グラフで、それを2ヶ月分だけ (右に) ずらしたデータが緑の折れ線グラフで表示される (目がチカチカしそうなのはご容赦を)。 なお、2つのデータの時系列上のずれを「遅れ」または「ラグ」 (lag) とよぶ。 上の例では、元の時系列データと、ラグが 2 のデータとを重ねて表示 させたことになっている。
上の例にならい、ラグを 0 から 25 程度まで変化させてグラフを描かせてみよ。 ラグの大きさの違いによって、2つのグラフの位置関係や重なり具合にどのような変化が見られるか。 また元のデータには12ヶ月を単位とする周期的な変動があるはずだが、それはグラフではどのような特徴として表われているか。
ある時系列データに対し、ある時間幅だけずらしたデータと元のデータとの間で相関をとることを 「自己相関」 (autocorrelation) という。 自己相関は、そのデータが 過去の履歴との間にどのくらいの相関をもっているか? を示すものである。 またこれによって計算される相関係数を 「自己相関係数」 と呼び、ラグの値によって定まる関数である。
余力のある人は、このプログラム を参考にして (コピペ可)、ラグを 0 から適当な大きさまでとったときの自己相関係数の値を計算 してみよ。 さらに、ラグを横軸に、自己相関係数の値を縦軸にとったグラフ (「コレログラム」ともいう) を描き、その特徴を検討 してみよ。 また元のデータには12ヶ月を単位とする周期的な変動があるはずだが、それはコレログラムではどのような特徴として表われているか。 ちなみに実際にやってみるには例えば、
$ gfortran enshu06z.f90として correlogram.dat というファイルに結果を出力し、その後 gnuplot を使って
$ ./a.out > correlogram.dat
gnuplot> plot [0:60] "correlogram.dat" with linespointsとかするとよい。
ここで挙げたように、あるデータと「少しの時間幅だけずらしたデータ」を重ね合わせ、その重なり具合によって時間変動の周期を検討することができる。 もちろん、あるデータと「少しの時間幅だけずらした 別のデータ」との重なり具合を調べることも可能である。 このような相関を「相互相関」と呼び、地球科学に限らずいろいろな分野で使われている。
バンドパスフィルター
第7回 の演習で扱った「移動平均」の計算では、速く変動する (高周波) 成分を弱め、ゆっくり変動する (低周波) 成分をとり出す「ローパスフィルター」、及びそれから導いた「ハイパスフィルター」を施したデータの特徴をみてきた。 しかし、「ローパスフィルター」と「ハイパスフィルター」に続けて通してやる と、ある特定の範囲の周波数で変動する成分のみを取り出す ことができるはずだ。 このようなフィルターを バンドパスフィルター (band-pass filter) という。 そこで今回の課題として、バンドパスフィルターをプログラミングにより作ってみよう。
バンドパスフィルターによって「ある特定の範囲の周波数の成分のみを取り出す」操作は、地球科学的観測データの処理でもよく用いられている。 例えば重力異常の解析においても、測定値を適切なバンドパスフィルターに通しておくことで、注目している解析領域の規模 (水平方向の広がりや深さ) に応じた変動の成分のみを取り出しておくことが普通である。 このような操作をすることで、地表面付近の (細かすぎる) 構造の影響や、(解析領域の広がりを超えた) 大規模な地球内部の構造の影響をあらかじめ排除できるという利点がある。
ハイパスフィルターを作る
これまでは Gnuplot のグラフの中だけで「ハイパスフィルター」の効果を見てみたが、「バンドパスフィルター」を作るための準備として「ハイパスフィルター」に通すプログラムを作ってみよう。 ここでは簡単のため、「自分自身とその前後1つずつをとった移動平均」をもとにして、「ハイパスフィルター」を作ってみることにする。 必要な手順 (の要点だけ) を以下に示す。
「ハイパスフィルターに通した結果」を入れるために、kion1retsu_highpass1 という名前の器を使うことを宣言する。 ただし kion1retsu や kion1retsu_lowpass1 という器の宣言のときと同様に、「入れられるデータの数」については必要分よりも多め にとってある (ここでは 1800/12=150年分)。
r e a l , d i m e n s i o n ( 1 8 0 0 ) : : k i o n 1 r e t s u _ h i g h p a s s 1 全てのデータを読み込んだ後に作った do ループの中で、「前後1つずつを含めた移動平均」に加えて、移動平均の値を引いて作った「ハイパスフィルター」後の値も計算 する。 計算結果の出力は、2列目に平均温度の「生」の値、3列目に「前後1つずつを含めた移動平均」、4列目に「ハイパスフィルター」後の値、となるようにしている。 なお、「ハイパスフィルター」後の値の計算には、「もとの値 kion1retsu から移動平均した値 kion1retsu_lowpass1 を引く」方法を使うのが簡単であろう。
d o b a n g o u = 1 , k o s u u i f ( b a n g o u = = 1 ) t h e n k i o n 1 r e t s u _ l o w p a s s 1 ( b a n g o u ) = & ( k i o n 1 r e t s u ( b a n g o u ) & + k i o n 1 r e t s u ( b a n g o u + 1 ) ) / 2 . 0 e l s e i f ( b a n g o u = = k o s u u ) t h e n k i o n 1 r e t s u _ l o w p a s s 1 ( b a n g o u ) = & ( k i o n 1 r e t s u ( b a n g o u - 1 ) & + k i o n 1 r e t s u ( b a n g o u ) ) / 2 . 0 e l s e k i o n 1 r e t s u _ l o w p a s s 1 ( b a n g o u ) = & ( k i o n 1 r e t s u ( b a n g o u - 1 ) & + k i o n 1 r e t s u ( b a n g o u ) & + k i o n 1 r e t s u ( b a n g o u + 1 ) ) / 3 . 0 e n d i f k i o n 1 r e t s u _ h i g h p a s s 1 ( b a n g o u ) = & k i o n 1 r e t s u ( b a n g o u ) - k i o n 1 r e t s u _ l o w p a s s 1 ( b a n g o u ) w r i t e ( * , ' ( i 4 , 3 f 6 . 2 ) ' ) b a n g o u , k i o n 1 r e t s u ( b a n g o u ) , & k i o n 1 r e t s u _ l o w p a s s 1 ( b a n g o u ) , k i o n 1 r e t s u _ h i g h p a s s 1 ( b a n g o u ) e n d d o
これらを参考にして、ハイパスフィルターの機能を追加したプログラムを作成せよ [みほん]。 そのプログラムを用いて (例の如く gfortran enshu06e.f90 などとして a.out を新たに作った後に)
$ cat matsuyama-kion2.txt | ./a.outと実行し、正しい (言い換えれば「自分の意図した通りの」) 出力が得られているかどうかを確認せよ。 なお、正しく計算できていれば、出力の先頭3行と末尾3行はそれぞれ以下のようになっているはずである。
1 4.80 6.15 -1.35 2 7.50 7.03 0.47 3 8.80 10.30 -1.50
1606 19.40 20.50 -1.10 1607 14.80 14.50 0.30 1608 9.30 12.05 -2.75正しい出力が得られていることが確認できたら、
$ cat matsuyama-kion2.txt | ./a.out > matsuyama-kion5.txtと実行することにより、出力を matsuyama-kion5.txt という名前のファイルに書き出させてみよ。 さらに、Gnuplot を使って、
gnuplot> plot [1488:1608] "matsuyama-kion3.txt" using 1:($2-$3) with lines, "matsuyama-kion5.txt" using 1:4 with linesなどとしてみることにより、2通りの方法で求めた「ハイパスフィルター」後の値が一致していることを確認せよ。
バンドパスフィルターを作る
ではいよいよ、「バンドパスフィルター」後の値を計算するプログラムを作ってみよう。 ここでは、「ハイパスフィルター」後の値を「ローパスフィルター」に通す、という方法で求めてみることにしよう (もちろん逆でもOKだが)。 必要な手順 (の要点だけ) を以下に示す。
「バンドパスフィルターに通した結果」を入れるために、kion1retsu_bandpass1 という名前の器を使うことを宣言する。 ただしこれまでに宣言してきた器と同様に、「入れられるデータの数」については必要分よりも多め にとってある (ここでは 1800/12=150年分)。
r e a l , d i m e n s i o n ( 1 8 0 0 ) : : k i o n 1 r e t s u _ b a n d p a s s 1 ハイパスフィルター後の値 kion1retsu_highpass1 の移動平均 をとる。 この操作によって、もとのデータをハイパスフィルターとローパスフィルターに続けて通す ことになるから、結果的にバンドパスフィルターに通したことになっている。 ただしこの操作は、ハイパスフィルターを計算する do ループの後 に、新しく do ループを作って移動平均を計算させる ようにするのがよい。
計算結果を出力する。 計算結果の出力は、2列目に平均温度の「生」の値、3列目に「前後1つずつを含めた移動平均」、4列目に「ハイパスフィルター」、5列目に「バンドパスフィルター」となるようにしている。 また結果の出力にあたっては、以下の例のように、プログラムの末尾にそれ専用の do ループを作ってやる ほうが明解であろう。
d o b a n g o u = 1 , k o s u u w r i t e ( * , ' ( i 4 , 4 f 6 . 2 ) ' ) b a n g o u , k i o n 1 r e t s u ( b a n g o u ) , & k i o n 1 r e t s u _ l o w p a s s 1 ( b a n g o u ) , k i o n 1 r e t s u _ h i g h p a s s 1 ( b a n g o u ) , & k i o n 1 r e t s u _ b a n d p a s s 1 ( b a n g o u ) e n d d o
これらを参考にして、バンドパスフィルターの機能を追加したプログラムを作成せよ [みほん]。 そのプログラムを用いて (例の如く gfortran enshu06f.f90 などとして a.out を新たに作った後に)
$ cat matsuyama-kion2.txt | ./a.outと実行し、正しい (言い換えれば「自分の意図した通りの」) 出力が得られているかどうかを確認せよ。 なお、正しく計算できていれば、出力の先頭3行と末尾3行はそれぞれ以下のようになっているはずである。
1 4.80 6.15 -1.35 -0.44 2 7.50 7.03 0.47 -0.79 3 8.80 10.30 -1.50 0.01
1606 19.40 20.50 -1.10 0.43 1607 14.80 14.50 0.30 -1.18 1608 9.30 12.05 -2.75 -1.22
正しい出力が得られていることが確認できたら、
$ cat matsuyama-kion2.txt | ./a.out > matsuyama-kion6.txtと実行することにより、出力を matsuyama-kion6.txt という名前のファイルに書き出させてみよ。 また、本日の演習の出席確認のため、このプログラムを 末尾のフォーム 経由で提出せよ。
さらに、Gnuplot を使って、
gnuplot> plot [1488:1608] "matsuyama-kion6.txt" using 1:4 with lines, "matsuyama-kion6.txt" using 1:5 with linesなどとしてみることにより、単なるハイパスフィルター後のデータと、バンドパスフィルター後のデータを比較してみよ。 両者のグラフの間にはどのような違いが見られるか。
余力がある人向きのおまけ: サブルーチン副プログラム
先に作ったバンドパスフィルターのプログラムの中には、「ローパスフィルターに通す」手順が2回含まれていた。 それら2箇所のローパスフィルターは、入力するデータは違っているけれど、手順そのものは全く同じ だった。 その違いのためだけに、ほとんど同じ命令を何度もプログラムに入力しないといけないのは面倒だし、どうせなら 共通で使える部分を使い回せる ようにしておきたいと思うのが普通であろう。
このような場合 (だけに限った用途ではないけど) には、副プログラム (subprogram) という機能を利用するのが有効である。 例えば Fortran には、サブルーチン (subroutine) と呼ばれる副プログラムが用意されている。 ここ に挙げたプログラム (コピペ可) は、サブルーチンの使い方の例である。 ここでは、移動平均によってローパスフィルターに通す部分を「内部サブルーチン」として独立させてある。 またこのサブルーチンの呼び出し (call) がプログラムの中で2回あるのだが、それぞれの呼び出しで異なった入力を与える ことにより、2種類のデータ (「もとのデータ」と「ハイパスフィルター後のデータ」) の両方に対して 同じ手順の処理を実行させる ことができている。 サブルーチンの詳しい説明は ここ を参照のこと。 共通部分をうまく「使い回し」たことにより、プログラムの「見た目」が劇的にすっきりしたことが分かるだろう。
サブルーチンに代表される「副プログラム」を用いることには、上に挙げた「共通部分の使い回し」の他にも数多くのメリットがある。 そのため、Fortran 90 以外のプログラミング言語で書かれたものも含め、世の中にある「実用的」なプログラムはたいがいが (何らかの形で) 副プログラムを利用して作成されている。 ということで、この先 (卒業研究などで) プログラミングを「バリバリに」使う腹積もりの人は、副プログラムの使い方に慣れておくことをおススメする。
まだまだ余裕のある人は、前回の最後で扱った「移動平均をとるデータの個数が奇数ならいくらでもOK」、あるいは「偶数でもOK」という実用的なプログラムを修正し、サブルーチンを使って書き直してみよ。
さらに余力がある人向きのおまけ: フィルターの「線形性」
ここで使用している「移動平均」に基づいたバンドパスフィルターでは、ローパスフィルターとハイパスフィルターの順序によらず同じ出力が得られる。 余力のある人は、最初にローパスフィルター、次にハイパスフィルターに通すことにより、バンドパスフィルターを作成してみよ。 また実際に両方の順番でバンドパスフィルターの計算を行い、答が変わらないことを確認せよ。 ただし実際にこの2通りでバンドパスフィルターを作るならば、例えば こんな感じ (コピペ可) のようにして、サブルーチンをうまく使うことをおススメする。
出席確認用プログラム提出フォーム
Moodle のコースページ 以下の 課題ファイル提出 により、各自で作成した (enshu06f.f90 に相当する) 「自分と前後1つずつ」でとった移動平均に基づくバンドパスフィルターの処理をする Fortran 90 プログラムの最終版を提出せよ。