13回目


前回の補足説明
ミルンの方法とアダムスの方法については前回省略しましたので,補足説明をします.
ルンゲクッタ法は基本になりますが,4次までの近似で,なおかつ多角形的なものであることはオイラー法と
変わりがありません.
その欠点を経験則的に解決したのがルンゲクッタ・ジル法であり,係数が細かく調整してあります.それでも,
微分方程式の曲線のカーブの変化が大きい時などは誤差が蓄積します.
ミルンの方法,アダムスの方法はこの問題点を解決する方法であり,関数中に設けたepsというパラメータが
あることに着目してみて下さい.
これは設定した誤差であり,一ステップ当たりの誤差がeps以下になるまでルンゲクッタによる解を修正すると
いうものです.
ルンゲクッタによる4ステップ分の計算をまず行い,これらの解を基に修正を1ステップごとに繰り返していくと
いうものです.ですから計算量は多くなります.epsを小さくすればするほど計算量が増える可能性が大きく
なるといえます.
ミルンの方法→comp_mlnm.c
アダムスの方法→comp_admsm.c


早いもので,13回を迎えました.
復習をかねて,関数の使い方,ポインタの使い方,グローバル変数の使い方と数値演算を学ぶために,
フィルタについて学習します.

フィルタとは何か?
皆さんも,実験や卒業研究で様々な測定を行うと思います.測定すると,必ず悩まされるのは,自らが
測定したいものをきちんと測れているのか?という問題です.
波形の生データを観察しても,何やら細かい波形が混ざっていたり,およそ測るものと関係なさそうな
うねりが延々と続いたりして,欲しい信号が取り出せていないと感じることもあるかと思います.
この,自身が測っているものとは関係なさそうなものこそが,ノイズと呼ばれるものです.
詳細は割愛しますが,以下のような測定システムの場合,ノイズはいろいろなところから発生すると考え
られます.


              図1 ノイズの経路

実験装置をきちんと作ればこのようなことは少なくなるのですが,それでも回避できないノイズは起こり
得るのです.従って,計測をする際には先に扱ったようなフーリエ変換を行うことで周波数成分の構成を
解析するのも一策ですが,対象とする信号の周波数(振動数)がある程度分かっているような場合には,
対象以外の信号を場外するような工夫をするのも大変有効なものといえます.
ここでは詳しくは述べませんが,ローパスフィルタ,ハイパスフィルタ,バンドパスフィルタ,ノッチフィルタ,
そのほかに,ローシェルフフィルタ,ハイシェルフフィルタ,ピーキングフィルタなどがあり,不要な帯域の
信号を減衰させ,必要な情報を取り出しやすくすることができます.

LPFのソースファイル→lpf.c

このソースファイルは2次バタワースIIR双一次変換フィルタというもので,デジタルフィルタとしての動作を
するところが特徴です.パラメータをいじって試してみて下さい.
詳細を話すと長くなりますので,概要だけにとどめておきます.

ランダムなノイズ
プログラムに用いる,ある関数にノイズを乗せた信号を作ってみます.今回は乱数を発生する方法を学び
ましょう.

ソースファイルを見て下さい.→random.c

ノイズは乱数列を使って発生させようというアイディアです.ランダムノイズですので偏りのない(正規では
ない)信号を乗せてみようというものです.
ランダム関数は,コンパイラに実装されている(エクセルにも実装されている)ものですが,処理系によって
差異があったりするようです.中身を知ってもらうため,ソースファイルのような関数を用意しました.

いちばんはじめのrand()がその中身であり,ある初期値に大きな数値を掛けて割り算した余りを戻すという
ものです.srand()もしくはlsrand()は初期値を変えることで乱数列を変化させるための関数です.但し,その
初期値がランダムになるような工夫が必要となります.
ポピュラーな方法としてはPC内の時刻(西暦,月,日付,時,分,秒)を取得し初期値とします.もちろん,
時刻に依存するから乱数かというと正確ではなくなりますが,今の正確な時刻はユーザーは分からない
ものとすれば,一応乱数といってよいと考えているからです.
もちろん,1秒以内に乱数列の初期化を2回以上行ってはいけないことは間違いありません.
あるサイン波にノイズを乗せるにはこのソースファイル内では,サイン波にノイズ成分を足し合わせています.
変調ノイズにしたければ,サイン波にノイズ成分を掛けてもいいと思います.
ノイズを乗せた信号を用いて,先程のフィルタに組み込めば,フィルタの効果が分かるかと思います.

プリプロセッサ
上記のプログラムにはそれとなくプリプロセッサを組み込んでいます.情報処理・演習1でも学んでいると思い
ますが,プログラムのはじめに,プリプロセッサを記述します.
#include <stdio.h>や#include<math.h>はお馴染だと思いますが,これもプリプロセッサです.
これはstdio.hというヘッダファイルをこのプログラムでは読み込みます,という意味です.
なお,#ifndef ○○のように,未だ○○を定義していなければ,#define ○○にて,○○を定義するとか,
便利な使い方ができるようになっています.例えば,RANDという定数を定義するなら,以下のようになります.

#ifndef RAND
#define RAND
#endif

グローバル変数を用いるのも方法ですが,プリプロセッサを用いて定数を文字で示してしまうこともできるので
必要に応じて使ってみて下さい.宣言は簡単です.重力加速度の例を示します.

#undef G_CONST
#define G_CONST 9.80665
#endif

一行目は,既に宣言している可能性のあるものを一旦解除し,二行目で定義し直しています.
尤も,既に宣言しているのを新たに無効にしているので本当に間違いないか注意深くプログラムを作る必要が
あります.


課題1 ミルンの方法,アダムスの方法を用いて前回と比較せよ.ソースファイルは以下のものを参考にせよ.
ミルンの方法 →comp_mlnm.c
アダムスの方法 →comp_admsm.c

課題2 ステップ波形に対し,ローパスフィルターを掛けてみよ.
その結果,波形はどのようになったか観察せよ.ファイルはlpf.cを利用すればよい.

課題3 ノイズが1/10乗った波形を作り,その波形を観察せよ.ソースファイルはrandom.cを使用してみること.
もし余力があれば,ローパスフィルタを掛けてみよ.lpf.cとrandom.cを組み合わせれば可能である.