久しぶりに大崎順彦(1994)「新・地震動のスペクトル解析入門」(ISBN:9784306032705)を読み返す機会があり、エルセントロ波の振幅スペクトルをSCILABとGMT ver.5 で描画させてみました。
右下の片対数グラフを描く際の設定ですが、 -JX9cl/5c は横軸対数目盛で9cm、縦軸線形目盛で5cmを示しています。また目盛の設定 -Bxa2f3 のうち a2は1と2と5に数値を表示、f3は全て(1,2,3・・・9)×10nに「ヒゲ」を表示することを指示しています。対数軸の場合、a,f,gの後に指定できる値は1か2か3しかないのだそうです。
#!/bin/bash temp=`mktemp ./XXXX.tmp` gmt gmtset FONT_ANNOT_PRIMARY 8p,GothicBBB-Medium-UniJIS-UTF8-H,black gmt gmtset FONT_LABEL 10p,GothicBBB-Medium-UniJIS-UTF8-H,black gmt gmtset FONT_TITLE 16p,GothicBBB-Medium-UniJIS-UTF8-H,black cat エルセントロ.txt |\ sharp_cancel |\ awk 'NR<=800{print $1,-$2}' > $temp gmt psxy $temp -JX20c/5c -R0/8/-400/400 -BnSeW+t"El Centro NS" -Bxa1+l'秒' \ -Bya100g1000+l'gal' -K -X6c -Y12c > test.ps scilab -nw << _SCILAB_ fd=mopen('$temp','r'); [n,t,a0]=mfscanf(-1,fd,'%f %f'); mclose(fd); ns0=size(a0,'r'); // ns0=800 fs=100; // サンプリング周波数 // 後続のゼロ a=[a0;zeros(1024-ns0,1)]; ns=size(a,'r'); // ns=1024(2の10乗) f=[0:fs/ns:fs]; // 周波数軸の設定 x=fft(a,-1); // FFTの実行 fd=mopen('$temp','w'); for k=1:ns, mfprintf(fd,'%f %f\n',f(k),abs(x(k))/fs); // 振幅スペクトルはFFTの結果をサンプリング周波数で割る end, mclose(fd); _SCILAB_ gmt psxy $temp -JX9c/5c -R0/10/0/200 -BnSeW+t"フーリェスペクトル(周波数表示)" \ -Bxa1+l'ヘルツ' -Bya50+l'gal*秒' -O -K -X0c -Y-8c >> test.ps cat $temp |\ awk '0.2<=$1&&$1<=20{print 1/$1,$2}' |\ gmt psxy -JX9cl/5c -R0.02/5/0/200 -BnSeW+t"フーリェスペクトル(周期表示)" -Bxa2f3+l'秒' \ -Bya50+l'gal*秒' -O -X11c -Y0c >> test.ps rm $temp convert -density 200 test.ps エルセントロ.jpg