SCILABとGMT ver.5 でフーリェ振幅スペクトルの描画

久しぶりに大崎順彦（1994）「新・地震動のスペクトル解析入門」(ISBN:9784306032705)を読み返す機会があり、エルセントロ波の振幅スペクトルをSCILABとGMT ver.5 で描画させてみました。

右下の片対数グラフを描く際の設定ですが、 -JX9cl/5c は横軸対数目盛で9cm、縦軸線形目盛で5cmを示しています。また目盛の設定 -Bxa2f3 のうち a2は1と2と5に数値を表示、f3は全て(1,2,3・・・9)×10ⁿに「ヒゲ」を表示することを指示しています。対数軸の場合、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