ルンゲクッタ法の説明と実装例[C言語]

3日かけてルンゲクッタ法を理解する所までにきました．

とは言うものの，この記事では1階常微分方程式についての4次ルンゲクッタについて説明しています．

前回，前々回の記事である

・オイラー法と修正オイラー法の説明と実装例
・ニュートン・コーツの公式とシンプソンの公式の導出と実装

を読むと今回の記事で説明するルンゲクッタ法の導出部分が理解しやすいかと思います．

1.(4次)ルンゲクッタ法の簡単な説明

詳しい説明はしたところであまり興味がない人が多い気がするので，詳しい説明は最後にすっ飛ばして，実装するにおいて最低限の部分だけ示しておきます．

次の初期値がわかっている1階常微分方程式を解く事を考えます．
$$\begin{eqnarray*} f(x_0) &=& y_0 \\ \frac{df(x)}{dx} &=& g(x,f(x)) \end{eqnarray*}$$ この時，4次ルンゲクッタ法では次のようにして計算します． $$\int_{x_k}^{x_{k+1}} g(x,f(x)) dx \approx \frac{1}{6}(k_1+2k_2+2k_3+k_4 )\Delta x$$ ただしk1,...,k4は $$\begin{eqnarray*} k_1 &=& g(x_k,f(x_k)) \\ k_2 &=& g(x_{k+\frac{1}{2}},f(x_k) + k_1 \frac{\Delta x}{2}) \\ k_3 &=& g(x_{k+\frac{1}{2}},f(x_k) + k_2 \frac{\Delta x}{2}) \\ k_4 &=& g(x_{k+1},f(x_k)+k_3 \Delta x) \end{eqnarray*}$$ として計算しています．
この辺の導出は一番最後でしているので気になる人は見てってくださいな．

2.実装例

何を実装しようか悩んで適当にWikipediaを漁っていたらロジスティック方程式というモデルがあったので，これのシミュレーションをしてみたいと思います．

Wikipediaによれば，生物の個体数がどのように変化していくかを説明するモデルの一つであるとのこと．

求める常微分方程式はこれ． $$\frac{dN(t)}{dt} = r \frac{K-N(t)}{K} N(t)$$ rは一匹の個体の増加率，N(t)は時間tにおける生物の個体数，Kは個体数のしきい値で，最終的に安定するNの値．

例えば水槽なら最大でも20匹ぐらいしか入らない大きさならK=20みたいな値になると思う．

ちなみに厳密解は普通の変数分離法で多少計算が面倒だけど求まる．初期値をN_0とすると次のようになる． $$N(t) = \frac{N_0 K e^{rt}}{K - N_0 + N_0 e^{rt}}$$ 厳密解はさておき，この常微分方程式にルンゲクッタ法を適用する．

まず見通しを良くするために，1で示した式のg(x,f(x))に当たる部分g(t,N(t))を定義しておく． $$\begin{eqnarray*} \frac{dN(t)}{dt} &=& r \frac{K-N(t)}{K} N(t) \\ &=& g(t,N(t)) \end{eqnarray*}$$ そうすると$\Delta t$を微小な値として，$t_k = k \Delta t$とすると，漸化式はこのようになる． $$\int_{t_k}^{t_{k+1}} g(t,N(t)) dt \approx \frac{1}{6}(k_1+2k_2+2k_3+k_4 )\Delta t$$ ただしk1,...,k4は次の式で計算する． $$\begin{eqnarray*} k_1 &=& g(t_k,N(t_k)) \\ k_2 &=& g(t_{k+\frac{1}{2}},N(t_k) + k_1 \frac{\Delta t}{2}) \\ k_3 &=& g(t_{k+\frac{1}{2}},N(t_k) + k_2 \frac{\Delta t}{2}) \\ k_4 &=& g(t_{k+1},N(t_k)+k_3 \Delta t) \end{eqnarray*}$$ それではC言語による実装．関数g,を定義することで，数式とほぼ変わらないプログラムが書けます．実装の係数は適当なのでアレですが．

	/* 4次ルンゲクッタ法によるロジスティック方程式の解の計算 */
	#include <stdio.h>
	#define K 1000.0 // 個体数のしきい値
	#define r 0.9 // 1個体の増加率
	#define N_LOOP 10000 // ループ回数
	// 積分する関数g(t,N)
	double g(double t,double N_t){
	return r*(K-N_t)*N_t/K;
	}
	int main(void){
	int k; // ループカウンタ
	double dt = 1e-3; // 刻み幅Δt
	double N = 1.0; // 個体数t（初期値を代入しておく)
	double k1,k2,k3,k4; // 計算に使う諸変数
	double a,m,b; // t_k,t_k+0.5,t_k+1を格納する
	for( k = 0 ; k < N_LOOP ; k++ ){
	// t_k,t_k+0.5,t_k+1を計算
	a = k*dt;
	m = (k+0.5)*dt;
	b = (k+1.0)*dt;
	// k1,...,k4の値を計算
	k1 = g(a,N);
	k2 = g(m,N+k1*dt*0.5);
	k3 = g(m,N+k2*dt*0.5);
	k4 = g(b,N+k3*dt);
	// 積分して代入
	N = N + (dt/6.0)*(k1+2.0*k2+2.0*k3+k4);
	// 値を表示
	printf("%lf,%e\n",b,N);
	}
	return 0;
	}

view raw RK4_KogisticEquation.c hosted with ❤ by GitHub

シミュレーション結果．

このシミュレーションだとほとんど違いがわかりませんね．誤差をプロットするとこうなります．

オーダーが1e-5といったところでしょうか．誤差が振動しているのがわかります．シミュレーションの精度としては十分なのでしょう．

4次のルンゲクッタ法では4次の項まで級数展開が一致しているとかで，誤差は$O(\Delta t^5)$に比例して小さくなるとのことです．

流石にルンゲクッタ法でシミュレーションすれば精度どうなの？とか言われても大丈夫そう……？まだその辺はよく知らないけれど．

3.ルンゲクッタ法の説明

手元の本には書かれていないのではっきりしたことはわからないのですが軽く調べたところによると，修正オイラー法は2次ルンゲクッタ法とも呼ぶらしいです．

なので2次ルンゲクッタについては前回の記事を参照してください．ここでは4次ルンゲクッタについて説明します．

今回はとりあえず1階微分方程式についての4次ルンゲクッタについて説明するので，次の常微分方程式について考えます．
$$f(x_0) = y_0 \\ \frac{df(x)}{dx} = g(x,f(x))$$ ここで　$x_k = x_0 + k \Delta x$　と置くと　$f(x_{k+1})$　と　$f(x_{k})$　の漸化式を得ます． $$\begin{eqnarray*} \int_{x_k}^{x_{k+1}} \frac{df(x)}{dx} dx &=& f(x_{k+1})-f(x_k) \\ &=& \int_{x_k}^{x_{k+1}} g(x,f(x)) dx \\ \\ ∴f(x_{k+1}) &=& f(x_k) + \int_{x_k}^{x_{k+1}} g(x,f(x)) dx \end{eqnarray*}$$

この $\int_{x_k}^{x_{k+1}} g(x,f(x)) dx$ の部分をどう計算するかによってオイラー法，修正オイラー法，ルンゲクッタ法が分類できます．

オイラー法では矩形積分，修正オイラー法では台形積分，ルンゲクッタ法ではシンプソンの公式を使います．

シンプソンの公式は次のようになります． $$\int_{x_k}^{x_{k+1}}h(x) dx \approx \frac{1}{6}(h(x_k)+4h(x_k+\frac{\Delta x}{2}) + h(x_{k+1}) )\Delta x$$ これを使って積分値を計算すれば良いわけです．$x_{k+\frac{1}{2}} = x_k+\frac{1}{2}\Delta x \\ $ として $$\begin{eqnarray*} g_k &=& g(x_k,f(x_k)) \\ g_{k+\frac{1}{2}} &=& g(x_{k+\frac{1}{2}},f(x_{k+\frac{1}{2}})) \\ g_{k+1} &=& g(x_{k+1},f(x_{k+1})) \\ \end{eqnarray*}$$ という値を考えると， $\int_{x_k}^{x_{k+1}} g(x,f(x)) dx$ は次のように近似できます． $$\int_{x_k}^{x_{k+1}} g(x,f(x)) dx \approx \frac{1}{6}(g_k + 4g_{k+\frac{1}{2}} + g_{k+1} )\Delta x$$
しかしこの式には未知の値$f(x_{k+\frac{1}{2}})$や$f(x_{k+1})$を使っています．これらの値を矩形積分で求めて計算することもできますが，良い精度が得られません．

ルンゲクッタ法ではこれらの値を二段階の計算で計算します．
$$\begin{eqnarray*} k_1 &=& g(x_k,f(x_k)) \\ k_2 &=& g(x_{k+\frac{1}{2}},f(x_k) + k_1 \frac{\Delta x}{2}) \\ k_3 &=& g(x_{k+\frac{1}{2}},f(x_k) + k_2 \frac{\Delta x}{2}) \\ k_4 &=& g(x_{k+1},f(x_k)+k_3 \Delta x) \end{eqnarray*}$$ $g_{k+\frac{1}{2}}$に当たる部分をk2,k3という形で計算しています．これらを先ほどのシンプソンの公式に代入することで以下のような近似を行うものを4次のルンゲクッタ法と呼ぶそうです． $$\int_{x_k}^{x_{k+1}} g(x,f(x)) dx \approx \frac{1}{6}(k_1+2k_2+2k_3+k_4 )\Delta x$$

今日はここまで．

ニュートン・コーツの公式とシンプソンの公式の導出と実装[C言語]

ルンゲクッタ法はどうやって導かれたのか謎だったのですが，『情報処理入門コース　数値計算』によれば，シンプソンの公式から来ているとのこと．
なのでシンプソンの公式がどういう経緯でできたのかメモしていきます．

※数式がうまく表示されない時はPCのブラウザから見ることを推奨します

1.ニュートン・コーツの積分公式

まずラグランジュ補間の公式が与えられているものとして使います．f(x)が近似したい式，p(x)が近似式．x_0,...,x_nが関数f(x)上のわかっている点． $$p(x)=\sum_{i=0}^{n}f(x_i)N(i,x) \\ N(i,x) = \prod_{j=0,j \neq i}^{n} \frac{x-x_j}{x_i-x_j}$$ ラグランジュ補間の出処はよくわかんないけど関数N(x)がクロネッカーのデルタ関数のようになっていることから意味をお察しください．とりあえずラグランジュ補間については前書いたのでそっちも見てくださいな．

次に，このラグランジュ補間で求めた近似式を用いて関数f(x)を区間[a,b]で積分する事を考えます． $$a = x_0,x_1,x_2,...,x_n = b$$ という様に区間をn分割し，上記のラグランジュ補間を適用すると次のようになります． $$\begin{eqnarray*} S &=& \int_{a}^{b}f(x)dx \\ & \approx & \int_{a}^{b}p(x)dx \\ &=& \int_{a}^{b} \sum_{i=0}^{n}f(x_i)N(i,x) dx \\ \end{eqnarray*}$$ ここで少しわかりづらいのですが，x_iは定数なのでf(x_i)はxの関数では無いことから， $$\begin{eqnarray*} S & \approx & \int_{a}^{b} \sum_{i=0}^{n}f(x_i)N(i,x) dx \\ &=& \sum_{i=0}^{n}f(x_i) \int_{a}^{b}N(i,x) dx \end{eqnarray*}$$ という式が得られます．ここで関数N(i,x)はx_0,...,x_nが[a,b]をn等分した値，すなわち $$x_k = a + k\Delta x \\ \Delta x = \frac{b-a}{n}$$ と表せることを用いて，次の変数変換をすることを考えます． $$t = \frac{x-a}{\Delta x}\\ dt = \Delta x \, dx$$ ここで，N(i,x)は上記のtを代入して次のように変形したものをM(i,t)とする． $$\begin{eqnarray*} N(i,x) &=& \prod_{j=0,j \neq i}^{n} \frac{x-x_j}{x_i-x_j} \\ &=& \prod_{j=0,j \neq i}^{n} \frac{(a+ t \Delta x) - (a + j \Delta x)}{(a+i \Delta x) - (a + j \Delta x ) } \\ &=& \prod_{j=0,j \neq i}^{n} \frac{t-j}{i-j} \\ &=& M(i,t) \end{eqnarray*}$$ このN(i,x)を変形した関数M(i,t)を用いると先ほどの積分はこうなる． $$\begin{eqnarray*} S & \approx & \sum_{i=0}^{n}f(x_i) \int_{a}^{b}N(i,x) dx \\ &=& \sum_{i=0}^{n}f(x_i) \Delta x \int_{0}^{n}M(i,t) dt \end{eqnarray*}$$ M(i,t)は変数xや，区間[a,b]などに依存しないから， $$C_i = \int_{0}^{n}M(i,t) dt$$ と置いて，先に手計算で計算しておけば使い回しが効いて，最終的に次のようになります． $$S \approx \sum_{i=0}^{n}f(x_i) C_i \Delta x$$ この公式をニュートン・コーツの積分公式と呼ぶとのこと．

2.シンプソンの公式

うえで求めたニュートン・コーツの積分公式にn=2を代入したものが世間一般にシンプソンの公式と呼ばれているものとのこと．C_iの値は上式に突っ込んで $$\begin{eqnarray*} C_0 &=& \int_0^2 \frac{t-1}{-1} \frac{t-2}{-2} dt &=& \frac{1}{3} \\ C_1 &=& \int_0^2 \frac{t}{1} \frac{t-2}{-1} dt &=& \frac{4}{3} \\ C_2 &=& \int_0^2 \frac{t}{2} \frac{t-1}{1} dt &=& \frac{1}{3} \\ \end{eqnarray*}$$ となる．この公式のすごいところはこの係数をあらかじめ求めておけば面倒な計算を端折ることができるところだと思います．
この値をニュートン・コーツの積分公式に代入すると， $$\begin{eqnarray*} \int_{a}^{b} f(x) dx &=& \sum_{i=0}^{2}f(x_i) C_i \Delta x\\ &=& (\frac{f(x_0)}{3} +\frac{4f(x_1)}{3} +\frac{f(x_2)}{3}) \Delta x\\ &=& (\frac{f(x_0)}{3} +\frac{4f(x_1)}{3} +\frac{f(x_2)}{3}) {\frac{b-a}{2}}\\ &=& \frac{1}{6}(b-a)(f(x_0)+4f(x_1) +f(x_2)) \end{eqnarray*}$$ となります．

シンプソンの公式を図で表すとこんな感じ．青色の線が二次式で近似した曲線．その曲線の積分を求めてる．

ちなみにn=3の場合などはニュートン・コーツの公式のWikipediaに載ってたりするのでそちらを参照すると良いかと思います．

3.プログラム

抽象的な話はこの辺にして，実際にプログラムを書いてみます．試しに$f(x) = \frac{4}{x^2+1}$を[0,1]区間で積分してみます．

ちなみに$\int_{0}^{1}f(x)dx$の厳密解はπです．

まず $\Delta x = \frac{1}{n}$ とし，$x_k = k \Delta x$ と置くと， $$\begin{eqnarray*} \int_{0}^{1}f(x)dx = \sum_{k=0}^{n-1} \int_{x_k}^{x_{k+1}}f(x) dx \end{eqnarray*}$$ という微小区間 $[x_k,x_{k+1}]$ の積分となります．この微小な区間に対してシンプソンの公式を適用すると次のようになります． $$\int_{x_k}^{x_{k+1}}f(x) dx \approx \frac{1}{6}(f(x_k)+4f(x_k+\frac{\Delta x}{2}) + f(x_{k+1}) )\Delta x$$ 実際はf(x_k)の足し算がかぶっているのでもう少しスマートにできるのですが，今回は愚直にそのままこれを実装します．

	/* シンプソンの公式による1/(x^2+1)の計算 */
	#include <stdio.h>
	// 被積分関数
	double f(double x){
	return 4.0/(x*x+1.0);
	}
	int main(void){
	int k; // ループカウンタ
	int N = 1000; // 分割数
	double dx = 1.0/N; // 刻み幅Δx
	double S = 0.0; // 和を格納する変数
	double a; // x_kの値を格納する
	double b; // x_{k+1}の値を格納する
	double m; // a,bの中間の値を格納する
	for( k = 0 ; k < N ; k++ ){
	// 近似を行う3点を計算
	a = k * dx;
	m = a + dx / 2.0;
	b = (k+1) * dx;
	// f_k + g(x_k,f(x_k))をfに代入
	S += (1.0/6.0)*(f(a)+4.0*f(m)+f(b))*dx;
	// 値の表示
	printf("%lf,%e\n",b,S);
	}
	return 0;
	}

view raw SimpsonsRule.c hosted with ❤ by GitHub

実行結果

結果が3.141593なので円周率が概ねちゃんと求まっているので大丈夫そうです．

オイラー法と修正オイラー法の説明と実装例[C言語]

学科的に流体シミュレーションとか多分そのうちやると思うのですが，ルンゲクッタ法を理解していなかったのでとりあえず常微分方程式の数値計算の基本であるオイラー法と修正オイラー法について記事を書いてみます．

1.オイラー法

オイラー法は精度は悪いですが，非常に簡単で計算が高速です．研究などでは恐らく使われることはありませんが，ゲームプログラミング等をしているとシューティングの弾の運動などは大抵オイラー法が使われています．
今回は説明のための例として，次のような常微分方程式があったとします．今回は初期値f(x_0)がわかっているものとします． $$f(x_0) = y_0 \\ \frac{df(x)}{dx} = g(x,f(x))$$ ここで　$x_k = x_0 + k \Delta x$　と置くと　$f(x_{k})$　の値は次のようにして求められることがわかります． $$\begin{eqnarray*} f(x_{k}) &=& f(x_0) + \int_{x_0}^{x_k} g(x,f(x)) dx \\ &=& f(x_0) + \sum_{i=0}^{k-1} \int_{x_{i}}^{x_{i+1}} g(x,f(x)) dx \end{eqnarray*}$$ ここで　$f(x_{k+1})$　と　$f(x_{k})$　の漸化式を考えます． $$\begin{eqnarray*} \int_{x_k}^{x_{k+1}} \frac{df(x)}{dx} dx &=& f(x_{k+1})-f(x_k) \\ &=& \int_{x_k}^{x_{k+1}} g(x,f(x)) dx \\ \\ ∴f(x_{k+1}) &=& f(x_k) + \int_{x_k}^{x_{k+1}} g(x,f(x)) dx \end{eqnarray*}$$ この $\int_{x_k}^{x_{k+1}} g(x,f(x)) dx$ の部分がキモです．この部分をどうするかで精度が大きく変わってきます．

オイラー法ではg(x,f(x))が区間$[x_k,x_{k+1}]$でほとんど変化しないものとして，次のように近似しています． $$\int_{x_k}^{x_{k+1}} g(x,f(x)) dx \approx g(x_k,f(x_k)) \Delta x$$ これを漸化式に代入したものがオイラー法です． $$f(x_{k+1}) \approx f(x_k) + g(x_k,f(x_k)) \Delta x$$

2.修正オイラー法

オイラー法は非常に簡単なのですが精度がよろしくないです．なのでもう少し精度を上げる必要があります．

修正オイラー法はその改良版の一つです．もっとも，これもあまり制度が良いというわけではないですが．

先ほど， $\int_{x_k}^{x_{k+1}} g(x,f(x)) dx$ の部分がキモと言いました．オイラー法では関数g(x,f(x))があまり変化しないものとして計算していました．イメージとしてはこんな感じです．

矩形近似とか長方形近似とか呼ばれているタイプの積分方法です．見るからに精度がよろしくなさそうなのがわかります．

そこで台形公式というものを使って修正オイラー法というものが提案されました．イメージとしてはこんな感じです．

名前の通り台形で近似するというものです．数式で表すと次のようになります．台形の面積は((高さ)×（上辺）+(下辺))÷2という小学生の時に習った公式を用いて，

$$\int_{x_k}^{x_{k+1}} g(x,f(x)) dx \approx (g(x_k,f(x_k)) + g(x_{k+1},f(x_{k+1}))) \frac{\Delta x}{2}$$ という近似式で計算します．しかしこれを漸化式に代入すると $$f(x_{k+1}) \approx f(x_k) + (g(x_k,f(x_k)) + g(x_{k+1},f(x_{k+1}))) \frac{\Delta x}{2}$$ となり，$f(x_{k+1})$ を求めるために $f(x_{k+1})$ を使っているのでこれでは使えません．なので右辺の $f(x_{k+1})$ だけは，オイラー法で使った矩形近似を用いて次のように求めておきます． $$\widehat{f(x_{k+1})} = f(x_k) + g(x_k,f(x_k)) \Delta x$$ この　$\widehat{f(x_{k+1})}$　を用いて漸化式に代入したものが修正オイラー法です． $$f(x_{k+1}) \approx f(x_k) + (g(x_k,f(x_k)) + g(x_{k+1},\widehat{f(x_{k+1})})) \frac{\Delta x}{2}$$

3.利用例

ここまで抽象的すぎたのでここで一度，次の1階微分方程式を数値的に解いてみたいと思います． $$\frac{df(x)}{dx} - 2f(x) = x$$ まぁ適当に作った式なので意味は無いです．ちなみに厳密解はCを定数として $$f(x) = -\frac{1}{4}(2x+1) + Ce^{2x}$$ と表されます．ちなみに定数Cは関数f(x)の初期値f(0) = aとした時，$C = a + \frac{1}{4}$になります．
今回は初期値 f(0) = a = 1 として数値的に求めてみます．

上記の微分方程式を変形してみます． $$\begin{eqnarray*} \frac{df(x)}{dx} &=& x + 2f(x) \\ &=& g(x,f(x)) \end{eqnarray*}$$ これでg(x,f(x))が求まりましたので，先ほど示した漸化式に代入してみれば，オイラー法では $$f(x_{k+1}) \approx f(x_k) + (x_k + 2f(x_k)) \Delta x$$ となります.

修正オイラー法では，まず　$\widehat{f(x_{k+1})}$　が $$\begin{eqnarray*} \widehat{f(x_{k+1})} &=& f(x_k) + g(x_k,f(x_k)) \Delta x \\ &=& f(x_k) + (x_k + 2f(x_k)) \Delta x \\ \end{eqnarray*}$$ となることから，次の漸化式に代入してやれば良いです． $$f(x_{k+1}) \approx f(x_k) + (g(x_k,f(x_k)) + g(x_{k+1},\widehat{f(x_{k+1})})) \frac{\Delta x}{2} \\ $$ 全部数式を代入すれば,x_kとf(x_k)とΔxだけの式にできますが，数式が長くなるので省略．

もっとも，関数をプログラム内で定義してやればそんなに気にならないと思います．

4.精度の評価

ソースコードは長くなりますので先にオイラー法，修正オイラー法の精度を示しときます．

もっとも，評価と言ってもこの例での真値との差を見るだけですが．以下数値的に求めたものと，真値との差です．
刻み幅Δx=0.001，反復回数1000回です．euler's method 2の方が修正オイラー法です．適切な英語が見つからなかったので……．

1000回繰り返して差が10^-7オーダーですね．とりあえずそれらしい値は取れているのがわかります．

これを見る限り修正オイラー法はオイラー法と比べたら優秀そうではあります．
ただ傾きが指数関数的に増えていきそうなので後日ルンゲクッタ法を実装してそれと比較したいです．

5.ソースコード

C言語による実装．オイラー法はこんな感じです．

	/* Euler法による微分方程式 f' - 2f = x の計算 */
	#include <stdio.h>
	#include <math.h>
	double g(double x,double fx){
	return x + 2.0 * fx;
	}
	int main(void){
	int k; // ループカウンタ
	int N = 10000; // 繰り返し数
	double dx = 1e-3; // 刻み幅Δx
	double f = 1.0; // 関数f(x)の値(初期値)
	double x_0 = 0.0; // 初期値の時のxの値
	double x_k; // x_kの値
	for( k = 1 ; k <= N ; k++ ){
	x_k = x_0 + k * dx; // x_kの値を求める
	f = f + g(x_k,f) * dx; // f_k + g(x_k,f(x_k))をfに代入
	// 値の表示
	printf("%lf,%e\n",x_k,f);
	}
	return 0;
	}

view raw EulersMethod.c hosted with ❤ by GitHub

修正オイラー法はこんな感じです．

	/* 修正Euler法による微分方程式 f' - 2f = x の計算 */
	#include <stdio.h>
	#include <math.h>
	double g(double x,double fx){
	return x + 2.0 * fx;
	}
	int main(void){
	int k; // ループカウンタ
	int N = 10000; // 繰り返し数
	double dx = 1e-3; // 刻み幅Δx
	double f = 1.0; // 関数f(x)の値(初期値)
	double f_bar; // f(x_k+1)を矩形積分でとりあえず求めた値
	double x_0 = 0.0; // 初期値の時のxの値
	double x_k; // x_kの値
	for( k = 1 ; k <= N ; k++ ){
	x_k = x_0 + k * dx; // x_kの値を求める
	// まず矩形積分でf_k+1を見積もって
	f_bar = f + g(x_k,f) * dx;
	// 台形公式で積分する
	f = f + (g(x_k,f) + g(x_k + dx, f_bar)) * dx /2.0;
	// 値の表示
	printf("%lf,%e\n",x_k,f);
	}
	return 0;
	}

view raw EulersMethod2.c hosted with ❤ by GitHub

今日はここまで．

自己相関関数で音声の繰り返し部分を抽出する[Python]

フーリエ変換でミクを描いたものの，フーリエ変換をちゃんと勉強したことがなかったので『メカトロニクス入門シリーズ　信号処理入門』を読んでます．

その中で「自己相関関数」というものが紹介されており，非常に面白そうだったので実装してみることにします．

1.自己相関関数とは

「あー」というような定常な音の場合，同じような音の波形が繰り返されています．

すなわち，f(t)を時間tの時の音の大きさだとすると，f(t)には次のような性質を持つことがわかります． $$f(t) = f(t + nT)$$ もっとも，実際の音声では完全にイコールにはならないと思いますが，概ねこのような性質を持っています．こういった関数を周期Tの周期関数と呼びます．例えばサイン関数は周期2πの周期関数です．

この周期Tを求める手法として上記の本では「自己相関関数」を用いています．関数fに対する自己相関関数R_ffを，今回次のように定義します．
$$R_{ff}(\tau) = \frac{(f(t),f(t+\tau))}{||f(t)||\,\,||f(t+\tau)||}\\ $$ ただし，実装するときにどうするかをはっきりしたいので，関数の内積とノルムは次のように定義します． $$(f(t),f(t+\tau)) = \int_{t_0}^{t_1} f(t)f(t+\tau)dt \\ ||f(t)|| = \sqrt{\int_{t_0}^{t_1}f(t)^2 dt }\\ ||f(t+\tau)|| =\sqrt{ \int_{t_0}^{t_1}f(t+\tau)^2 dt} \\ $$
関数f(t)とf(t+τ)の相関係数みたいなものなので，関数のノルムで割るかどうかや，関数の平均値を減じるかどうかは人それぞれだと思います．

さてこの自己相関関数R_ffを実際に計算していくとこのようになります．

ブログ用．自己相関関数の値の変化． pic.twitter.com/Qmul8V82lH

— もりとにー (@TonyMooori) 2016年3月7日

τを変化していき，τ=nT ( n = 0,1,2,... )となる点で似たパターンとなるため，R_ffが極大値を取ることがわかります．

2.抽出アルゴリズム

ソースコードを貼る前に，実際にプログラムを組んで「あー」という音を読み込ませて，その音声に対する自己相関関数を求めるとこんな感じになります．

とまぁ現実の音をやってみるとそんなうまくいかないようです．結構ギザギザしてて嫌な感じです．

恐らくフラクタル図形のように「あー」という音の中にも似た部分があるようです．それはそれで面白いのですが……．

まぁ自分の目で見ればどれが目的の周期なのか一目瞭然なのですが，できれば自動的に判定したいところです．

とりあえずFFTで高周波成分を除去するハイパスフィルタを実装するとこんな感じになります．緑色の部分がハイパスフィルタで高周波成分を除去したもの．移動平均でも良いと思います．

そんでもってscipy.signalの極大値の要素番号をsignal.argrelmaxで求めて，しきい値0.2よりも大きい相関以上かつ極大となっている要素番号を取得すると，{5,  646, 1291, 1938, 2582, 3228, 3872, 4520} みたいな等差数列になっているので差を求めてやればそれっぽい周期が求められる．

その周期で音声の波形を10個ほど取り出すとこんな感じになった．

繰り返されている部分がいい感じに取り出せました．とは言うものの周期Tをint型に型キャストしているのでズレていってしまいますがまぁこんなもんでしょう．

この繰り返されている波形データだけで音声認識とかできないか考え中なところです．

3.プログラム

コメントは多めにしたつもりではありますがうまくいかない場合はコメントを参考にしつつ頑張って下さい．

	#coding:utf-8
	import numpy as np
	import pandas as pd
	import matplotlib.pyplot as plt
	import wave
	from scipy import signal
	"""
	参考
	[1] 波形を見る - 人工知能に関する断創録 http://aidiary.hatenablog.com/entry/20110519/1305808715
	[2] 雨宮好文 監修, 佐藤幸男 著『信号処理入門』 オーム社
	"""
	def read_wave(path):
	""" waveファイルを読み込む関数 """
	wf = wave.open(path , "r" )
	fs = wf.getframerate() # サンプリング周波数
	x = wf.readframes(wf.getnframes())
	x = np.frombuffer(x, dtype= "int16") / 32768.0 # -1 - +1に正規化
	wf.close()
	# fsはサンプリング周波数
	return (x,fs)
	def HPF(x,n_dim):
	""" ハイパスフィルタ，n_dimは打ち切る周波数 """
	c = np.fft.rfft(x) # 高速フーリエ変換
	c[n_dim:] = 0.0 # 高周波成分除去
	return np.fft.irfft(c)
	def cross_correlation(v1,v2,low,high,tau):
	""" ベクトルv1,v2の相互相関関数 """
	# v1の要素番号lowからhightまで
	# v2の要素番号low+tauからhigh+tauまでを取り出す
	x1 = v1[low:high]
	x2 = v1[(low+tau):(high+tau)]
	# ノルムを1にする
	x1 = x1 / np.linalg.norm(x1)
	x2 = x2 / np.linalg.norm(x2)
	# x1,x2の内積が相関係数となる
	return x1.T.dot(x2)
	if __name__=="__main__":
	# データ読み取り
	print("reading voice.wav")
	x,fs = read_wave("./voice.wav")
	# 自己相関関数を計算
	print("calculating R_ff")
	R = np.array([ cross_correlation(x,x,0,fs,tau) for tau in range(fs) ])
	# 周期を計算してく
	print("calculating T")
	# まずハイパスフィルタで高周波成分を除去
	R_HPF = HPF(R,len(R) / 20)
	# 自己相関関数はτが小さいほうが相関が高いため，τが小さい部分を取り出す
	R_HPF[:fs/5]
	# 極大値の要素番号を取得
	local_max, = signal.argrelmax(R_HPF)
	# 極大値の中で0.2よりも大きい部分の要素番号を取得
	threshold = 0.2 # しきい値
	local_max = local_max[ R_HPF[ local_max ] > threshold ]
	# この時 local_max は(うまく行けば)
	# array([ 5, 646, 1291, 1938, 2582, 3228, 3872, 4520])
	# のような等差数列になっているので各項の差を取る
	diff = np.diff(local_max)
	# diffは(うまく行けば)こうなっている
	# array([641, 645, 647, 644, 646, 644, 648])
	# 外れ値があるかもしれないので，平均値ではなく中央値を使うと良い
	T = (int)(np.median(diff))
	# 周期Tで取り出して表示
	for i in range(10):
	plt.plot(x[(i*T):(i*T+T)])
	plt.title("T = %d" % T )
	plt.show()

view raw extract_wave.py hosted with ❤ by GitHub