RukeとLuNaYuの日記
I know the truth.
I know whole.
And I...know you.
平凡な大学生活の日記です。時折まじめな長文を書く病気になります。興味がなければ読み飛ばしてください。
201707 << 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 >> 201709
スポンサーサイト (--/--/--(--) --:--:--)
上記の広告は1ヶ月以上更新のないブログに表示されています。
新しい記事を書く事で広告が消せます。
スクリプトゾク (2005/12/16(金) 00:51:23)
いくらか改良した。パワースペクトルはどの程度詳細が失われているかが明確ではない平滑化を止めて、ヒストグラムにした。
#ところで平滑化をする時に、(1,1,1,,,1)のようなベクトルとパワースペクトルとの畳み込みをとる事になった。そのために昨日バージョンではfftconvを利用していたのだが、これにfftとついていることからもわかるように、これは、フーリエ変換をする前の時系列データにある関数を乗じている操作そのものだ。で、これって窓関数って奴そのものじゃないか、と思い、octaveにはbartlettとかhamming,hanningとか色々ついてくるので、それをかけてやれば勝手にスペクトルが滑らかになるのかな、と試してみたのだがどうもうまく行かない。これらの窓関数って離散的なスペクトルを観測する事に主眼があるようで、この目的には適さないようだ。。。って理解で良いのかな?

せっかくだから出力結果。これが二分で生成される。わはは。自分が偉くなった気分だ。偉いのはoctaveだけど。


#Note that power spectrum S(f) is scaled so that
#the integral of S(f)df from 0Hz to +infinity Hz 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]';
freq=[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));
rowS=(deltaT/(2*pi*N))*fftdata.*conj(fftdata);
cor=shift([tau,real(ifft(rowS))*2*pi/deltaT],N/2);

#getting together several neighbor points of power spectrum
rowS=[freq(1:N/2),2*pi*2*real(rowS(1:N/2))];
smooth=floor(deltaT*N/(tau_f*40));
newN=floor(N/2/smooth);
rowS=rowS(1:newN*smooth,:);
S=zeros(newN,2);
for i=1:newN
S(i,:)=[sum(rowS((i-1)*smooth+1:i*smooth,1))/smooth,sum(rowS((i-1)*smooth+1:i*smooth,2))/smooth];
endfor

#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.3\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.6\n\
set format y \"10^{%%L}\"\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 xrange[*:*]\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.70,graph 0.84\n\
set bar 0\n\
plot \"histgram.subplot\" using ($1*1e6):2 notitle with histeps lw 0.7, \"histgram.subplot\" using ($1*1e6):2:(sqrt($2)) notitle with yerrorbars lt 1 lw 0.7 ps 0, bestgauusian(x/1e6) title \"best fit by gaussian\" with lines lt 1 lw 1.2\n\
\n\
set xrange[0:*]\n\
set yrange[*:*]\n\
set nolabel\n\
set size 0.7,0.7\n\
set origin 0,0.8\n\
set xlabel \"frequency f[Hz]\"\n\
set ylabel \"Power S(f)[V^2s]\"\n\
set nologscale x\n\
set logscale y\n\
set title ""\n\
plot \"powerspectrum.subplot\" using 1:2 notitle with histeps\n\
\n\
set xrange[*:*]\n\
set yrange[*:*]\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ヶ月以上更新のないブログに表示されています。新しい記事を書くことで広告を消せます。