RukeとLuNaYuの日記
I know the truth.
I know whole.
And I...know you.
平凡な大学生活の日記です。時折まじめな長文を書く病気になります。興味がなければ読み飛ばしてください。
200503 << 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 >> 200504
スポンサーサイト (--/--/--(--) --:--:--)
上記の広告は1ヶ月以上更新のないブログに表示されています。
新しい記事を書く事で広告が消せます。
スクリプト (2005/12/14(水) 05:31:59)
ブラウン運動の第一日目は、ぬるそうに見えて48000サンプリングポイントの時系列データを16種類取るので、解析と整理が糞大変だ。

んで、octaveとgnuplotのスクリプトをずんだか書いた。でっちあげ終わって見ると、(48000サンプリングポイントの時系列データのヒストグラム、パワースペクトル、自己相関関数の表示+パラメタや統計情報付きでのepsファイルへのプロット(三つの図が縦に並んだ形になる))x16が二分で済むようになった。かけた手間をかけたらこれは素晴らしい事だ。

再びgnuplotとoctaveに感謝。高級言語だが抽象言語ではないというのは、やっぱり重要だなあ。

あと、これはC言語の標準ライブラリの実装の問題なんだけど、printfの書式指定で%.3gとした時に、末尾の0が落とされるのはひどくないだろうか。何だか広く知られている謎仕様の一つのような気がするが、知らずにはまった(というか解決せずに放置)。あるいは合理的な理由があるのかな?

というわけでポリシーを曲げて書いた超ダーティーなスクリプトを載せておく。例えばこれをthermalnoise.mとしてbashでfor name in *.dat;do octave -qf noise.m dataname=$name;doneしてぼけーとしていると2分でプロットがだーっと生成されるのじゃ。
#Note that power spectrum S(w) is scaled so that
#the integral of S(w)dw from 0 to +infinity gives <V^2>.
#Self correlation function C(tau) is also scaled
#as to give <V(t)V(t+tau)>.

#If dataname=%%.dat or dataname=%% is specified
#in the command-line options, dataname.dat
#and dataname.params is processed.
#Otherwise the default name shown below is used (for debugging).
filename="group4a/10MOhm_100ms_50Hz_960s";

for i=1:__nargin__
head="dataname=";
if length(argv{i})>length(head) && strcmp(head,argv{i}(1:length(head)))==1
filename=argv{i}(length(head)+1:length(argv{i}));
endif
endfor
if strcmp(".dat",filename(length(filename)-3:length(filename)))==1
filename=filename(1:length(filename)-4);
endif

printf("processing %s.dat with %s.params...\n",filename,filename)

#time-sequence is read from .dat file
data=load(strcat(filename,".dat"));

#parameters
#R, tau_f, G1 is read from .params file
load(strcat(filename,".params"));
G0=5e4;
k_B=1.38e-23;
deltaT=data(2,1);
N=rows(data);

#normalizatoin
data=[data(:,1),data(:,2)/G0/G1];
v_max=10.0/G0/G1;

#prepare x-axis
div=40;
v=[-v_max+v_max/div:v_max/div:v_max-v_max/div]';
tau=[0:deltaT:deltaT*(N/2-1),-deltaT*(N/2):deltaT:-deltaT]';
omega=2*pi*[0:1/(N*deltaT):(N/2-1)/(N*deltaT),-(N/2)/(N*deltaT):1/(N*deltaT):-1/(N*deltaT)]';

#statistics
average=mean(data(:,2))
sigma=std(data(:,2))
sigma_sq=sigma^2
#temperature
T=tau_f*sigma_sq/k_B/R

#power and self-correlation
fftdata=fft(data(:,2)-average);
S=(deltaT/(2*pi*N))*fftdata.*conj(fftdata);
cor=shift([tau,real(ifft(S))*2*pi/deltaT],N/2);
#smoothing power spectrum
smooth=100;
S=fftconv(S,ones(smooth,1)/smooth)(floor(smooth/2)+1:N+floor(smooth/2),1);
S=[omega(1:N/2),2*real(S(1:N/2))];
#S=shift([omega,real(S)],N/2)

#histgram
histgram=[v,hist(data(:,2),v)'];
figure(1);
gplot histgram using 1:2 with histeps;

#power spectrum
figure(2);
gplot S using 1:2 with lines;

#self-correlation
cut=6;
figure(3);
gplot [-tau_f*cut:tau_f*cut] cor using 1:2 with linespoints;

#prepare for using gnuplot
save "histgram.subplot" histgram;
save "powerspectrum.subplot" S;
save "selfcorrelation.subplot" cor;

#perform some fitting
initialparams=sprintf("\n\
A=%g\n\
V_0=%g\n\
sigma_=%g\n\
",histgram(floor(rows(histgram)/2)+1,2),average,sigma);

gnuplot_commands=sprintf("\n\
A=0\n\
V_0=0\n\
sigma_=0\n\
f(x)=A*exp(-((x-V_0)**2)/(2*(sigma_**2)))\n\
fit [-%g*3:%g*3] f(x) \"histgram.subplot\" using 1:2:(sqrt($2)) via \"initialparams.fit\"\n\
update \"initialparams.fit\" \"bestparams.fit\"\n\
",sigma,sigma);

save "initialparams.fit" initialparams;
save "gnuplot_commands" gnuplot_commands;
system "gnuplot gnuplot_commands";
source "bestparams.fit";

initialparams=sprintf("\n\
sigma__=%g\n\
tau_f_=%g\n\
",sigma,tau_f);

gnuplot_commands=sprintf("\n\
sigma__=0\n\
tau_f_=0\n\
f(x)=(sigma__**2)*exp(-abs(x)/tau_f_)\n\
fit [-%g*3:%g*3] f(x) \"selfcorrelation.subplot\" via \"initialparams.fit\"\n\
update \"initialparams.fit\" \"bestparams.fit\"\n\
",tau_f,tau_f);

save "initialparams.fit" initialparams;
save "gnuplot_commands" gnuplot_commands;
system "gnuplot gnuplot_commands";
source "bestparams.fit";

#perform plotting
gnuplot_commands=sprintf("\
tau_f=%g; \n\
reset\n\
set terminal postscript eps enhanced\n\
set output \"%s.eps\"\n\
set size 0.7,2.1\n\
set multiplot\n\
set sample 5000 \n\
\n\
A=%g\n\
V_0=%g\n\
sigma_=%g\n\
bestgauusian(x)=A*exp(-(x-V_0)**2/2/sigma_**2)\n\
set size 0.7,0.7\n\
set origin 0,1.4\n\
set xlabel \"noise voltage V[{/Symbol m}V]\"\n\
set ylabel \"count\"\n\
set nologscale x\n\
set logscale y\n\
set yrange [1:*]\n\
set title \"R=%.3gM{/Symbol W}, {/Symbol t}_f=%.3gms.\"\n\
set label \"<V>=%.2g{/Symbol m}V\\n {/Symbol s}=%.3g{/Symbol m}V\\n\\nT={/Symbol t}_f{/Symbol s}^2/k_BR=%.3gK\" at graph 0.67,graph 0.84\n\
set bar 0\n\
plot \"histgram.subplot\" using ($1*1e6):2 notitle with histeps lw 1.2, \"histgram.subplot\" using ($1*1e6):2:(sqrt($2)) notitle with yerrorbars lt 1 lw 1.0 ps 0, bestgauusian(x/1e6) title \"best fit by gaussian\" with lines lt 1 lw 0.6\n\
\n\
set yrange[*:*]\n\
set nolabel\n\
set size 0.7,0.7\n\
set origin 0,0.7\n\
set xlabel \"frequency {/Symbol w}(2{/Symbol p}f)[s^{-1}]\"\n\
set ylabel \"Power S({/Symbol w})[V^2s]\"\n\
set logscale x\n\
set logscale y\n\
set title ""\n\
plot \"powerspectrum.subplot\" using 1:2 notitle with lines\n\
\n\
sigma__=%g\n\
tau_f_=%g\n\
f(x)=(sigma__**2)*exp(-abs(x)/tau_f_)\n\
set size 0.7,0.7\n\
set origin 0,0\n\
set xlabel \"Delay time {/Symbol t}[s]\"\n\
set ylabel \"Self-correlation function C({/Symbol t})[V^2]\"\n\
set nologscale x\n\
set logscale y\n\
set key right top\n\
set label \"{/Symbol s}'=%.3g{/Symbol m}V\\n {/Symbol t}_f'=%.3gms\" at graph 0.7, 0.85\n\
plot [-tau_f*6:tau_f*6] \"selfcorrelation.subplot\" using 1:2 notitle with lines lt 1, f(x) title \"{/Symbol s}'^2 exp(-{/Symbol t}/{/Symbol t}_f')\" with lines lt 2\n\
",tau_f,filename, A,V_0,sigma_,R*1e-6,tau_f*1e3,average*1e6, sigma*1e6,T,sigma__,tau_f_,sigma__*1e6,tau_f_*1e3);

save "gnuplot_commands" gnuplot_commands;
system "gnuplot gnuplot_commands";
スポンサーサイト
コメント
この記事へのコメント
コメントを投稿する

管理者にだけ表示を許可する
トラックバック
この記事のトラックバックURL
この記事へのトラックバック
(C)Copyright 2003-2007 by Ruke All rights reserved. Powered By FC2. VALID HTML? VALID CSS?
上記広告は1ヶ月以上更新のないブログに表示されています。新しい記事を書くことで広告を消せます。