Octave で FFT してみた

きのうの続き。久しぶりに計算してみたくなったので、きのう思い出した logistic map のパワースペクトルを描いてみた。fftw とか使えばいいんだがすっかり忘れてるので、どうせなら新しいものをってことで Octave-forge を Windows にインストールしてみた。MATLAB はむかし使ったことあるけど Octave は初めて。ダウンロードを始めてから1時間、マニュアルを見ながら、こういうスクリプトを書いた:

function retval = f(a,x)
	retval = a*x*(1-x);
endfunction

function retval = iter(n,a,x0)
	retval = [1:n];
	x = x0;
	for i = 1:10000
		x = f(a,x);
	endfor
	for i = 1:n
		x = f(a,x);
		retval(i) = x;
	endfor
endfunction

コンソールからはこういう風に使う

octave:1> source "logistic.m"
octave:2> plot(abs(fft(iter(2048, 3.62, 0.1)))(2:1024));

そうすると、Octave から呼び出された gnuplotパワースペクトルのグラフを描いてくれる。すげえ簡単。

ひさしぶりに描いてみて分かったが、カオスだったらどこでも対称ってわけでなくて a<3.65 くらい。その先にもところどころあるっぽい。それと完全に対称ってわけでもなくて対称っぽいという感じ。