2018-01-17

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

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

f:id:S_E_Hyphen:20180117160530j:plain

右下の片対数グラフを描く際の設定ですが、 -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

2018-01-16

お揚げさんの油抜き

すぐに忘れてしまうので備忘録的に…。

「市販のうすあげ　2枚入り」を、それぞれ半分に切って沸騰したお湯で5分間煮ます。

f:id:S_E_Hyphen:20180116104201j:plain

煮上がったらザルにあげて軽く水洗いします。パリパリだったおあげさんが、空気が抜けてしっとりした状態となります。

再びフライパンに戻して「かえし醤油　25cc(カンロレードル1杯分）」と水200ccを加えて再び加熱します。

ブリの照り焼きでも書いたように、調味料は、食材をそのまま食べるのに適量かどうかが目安です。一方、加える水の量はいい加減でよいのかもしれません。煮汁がきっちり無くなるまで茹でることが大切です。早く火を止めすぎると味が入りませんし、遅すぎると焦げつきます。今回は写真を撮っていて、うっかり焦げ臭くなってしまい慌てました。

f:id:S_E_Hyphen:20180116103352j:plain

きつねうどんのように、主に麺類のトッピングに使っています。刻んでご飯に載せると木ノ葉丼ですね。もちろん、すし飯を包むとおいなりさんですが不器用なので未だ作ったことはありません。

2018-01-12

位相の異なる正弦波の足し合わせ

図の上下とも、赤線は周波数0.5Hzの正弦波、青線は0.6Hzの正弦波です。どちらの正弦波も振幅は1.0としています。2つの正弦波を足し合わせたものが、両図とも太い黒線で示したものになります。この太黒線にFFTを実施して振幅スペクトルを描くと両者とも全く同じになります。

上下の違いは、青線で示した0.6Hzの正弦波の位相だけです。しかし、位相の異なる正弦波を足し合わせると、その最大値は全く違う値となります。

f:id:S_E_Hyphen:20180113075635j:plain

太黒線だけが与えられていて、仮にローパスフィルタで0.6Hzを取り除いたならば、上図の場合は最大振幅が小さくなりますが、下図の場合は除去後の方が大きな値となることでしょう。

#!/bin/bash
temp=`mktemp ./XXXX.tmp`
temp2=`mktemp ./XXXX.tmp`
otfile=test
gmt gmtset FONT_TITLE 16p,GothicBBB-Medium-UniJIS-UTF8-H

scilab -nw << _SCILAB_
fs=100; //サンプリング周波数
ns=8192; //サンプル数
t=[0:1/fs:(ns-1)/fs];

fp_a=0.5; // 正弦波の周波数
a=sin(2*%pi*fp_a*t);
fp_b=0.6;
phi_b=180*%pi/180;
b=sin(2*%pi*fp_b*t+phi_b);
c=a+b;
max(c)
fd=mopen('$temp','w');
for k=1:ns,
mfprintf(fd,"%f %f %f\n",t(k),a(k),b(k));
end
mclose(fd);

x=fft(c,-1);
fs=1/t(2);
f=[0:fs/ns:100];
fd=mopen('${temp2}','w');
for k=1:ns,
mfprintf(fd,"%f %f\n",f(k),abs(x(k)));
end
mclose(fd);

_SCILAB_

cat $temp |\
 awk '{print $1,$2}' |\
 gmt psxy -JX15c/5c -R0/1/-2.2/2.2 -BnSeW \
   -Bxa0.1+l"sec" -Bya0.5 -W0.5p,red -K > $otfile.ps
cat $temp |\
 awk '{print $1,$3}' |\
 gmt psxy -J -R -B -W0.5p,blue -O -K >> $otfile.ps
cat $temp |\
 awk '{print $1,$2+$3}' |\
 gmt psxy -J -R -B -W1p,black -O -K >> $otfile.ps
cat << + |\
gmt pslegend -R -J -Dx12c/4.8c+w2.5c/1.2c+jTL -F+p+gwhite -O -K >> $otfile.ps
N 1
S 0.2c s 0.3c 255/0/0 0.0p 0.5c 0.5Hz
G 0.02c
S 0.2c s 0.3c 0/0/255 0.0p 0.5c 0.6Hz
+
printf "0.8 -1.9 max=%.2f\n" 0.59 |\
 gmt pstext -J -R -F+jBL+a0.0+f12p,Helvetica,0/0/0 -O -K >> $otfile.ps 

cat ${temp2} |\
 gmt psxy -JX5c/5c -R0/1/0/2000 -BnSeW+t"振幅スペクトル" \
   -Bxa0.2+l"Hz" -Bya500 -W1.0p,black -X18c -O -K >> $otfile.ps

scilab -nw << _SCILAB_
fs=100; //サンプリング周波数
ns=8192; //サンプル数
t=[0:1/fs:(ns-1)/fs];

fp_a=0.5; // 正弦波の周波数
a=sin(2*%pi*fp_a*t);
fp_b=0.6;
phi_b=0*%pi/180;
b=sin(2*%pi*fp_b*t+phi_b);
c=a+b;
max(c)
fd=mopen('$temp','w');
for k=1:ns,
mfprintf(fd,"%f %f %f\n",t(k),a(k),b(k));
end
mclose(fd);

x=fft(c,-1);
fs=1/t(2);
f=[0:fs/ns:100];
fd=mopen('${temp2}','w');
for k=1:ns,
mfprintf(fd,"%f %f\n",f(k),abs(x(k)));
end
mclose(fd);

_SCILAB_

cat $temp |\
 awk '{print $1,$2}' |\
 gmt psxy -JX15c/5c -R0/1/-2.2/2.2 -BnSeW \
   -Bxa0.1+l"sec" -Bya0.5 -W0.5p,red -X-18c -Y8c -O -K >> $otfile.ps
cat $temp |\
 awk '{print $1,$3}' |\
 gmt psxy -J -R -B -W0.5p,blue -O -K >> $otfile.ps
cat $temp |\
 awk '{print $1,$2+$3}' |\
 gmt psxy -J -R -B -W1p,black -O -K >> $otfile.ps
cat << + |\
gmt pslegend -R -J -Dx12c/4.8c+w2.5c/1.2c+jTL -F+p+gwhite -O -K >> $otfile.ps
N 1
S 0.2c s 0.3c 255/0/0 0.0p 0.5c 0.5Hz
G 0.02c
S 0.2c s 0.3c 0/0/255 0.0p 0.5c 0.6Hz
+
printf "0.8 -1.9 max=%.2f\n" 1.98 |\
 gmt pstext -J -R -F+jBL+a0.0+f12p,Helvetica,0/0/0 -O -K >> $otfile.ps 

cat ${temp2} |\
 gmt psxy -JX5c/5c -R0/1/0/2000 -BnSeW+t"振幅スペクトル" \
   -Bxa0.2+l"Hz" -Bya500 -W1.0p,black -X18c -O >> $otfile.ps

rm $temp ${temp2}

2018-01-09

小豆シチュー

小豆粥を作るのに買った小豆を使ってみたくてシチューに入れてみました。

まず茹でこぼしまでを行った小豆を準備しておきます。

ホワイトシチューのパッケージに書かれてある指示通りの分量のじゃがいも、人参、玉ねぎ、鶏肉を炒めます。指示通りの分量の水を加えて、煮込み始める段階で小豆も加えます。あとは指示通りの時間煮込んだ後、ルーを割り入れます。

もっと強い発色を期待していましたが、意外と普通の色でした。

でも、植物性繊維は増えてヘルシーですよ。

f:id:S_E_Hyphen:20180109090858j:plain

2018-01-02

GMT ver.5で表の描画

GMT ver.5 で貸出記録の印刷や GMTでニュース見出しをグラフ化するで個別には実現していましたが、今回一般化してみました。

#!/bin/bash
cntl=`getparstr $# "$*" "cntl"` #入力する制御用ファイル
psfile=`getparstr $# "$*" "psfile"` # 出力するPSファイル

# カラムの数は制御用ファイルの１行目で決定します
num_column=`cat $cntl |\
 awk 'NR==1{print NF}'`
# 行の数は制御用ファイルの２行目の第１変数で決定します
num_row=`cat $cntl |\
 awk 'NR==2{print $1}'`
# 行の幅は制御用ファイルの２行目の第２変数で決定します
int_row=`cat $cntl |\
 awk 'NR==2{print $2}'`
hight=`bc -l << _BC_
 $num_row * $int_row
_BC_`
temp=`mktemp ./XXXX.tmp`
cat /dev/stdin > $temp # 標準入力を一時ファイルに持っておきます

rm $psfile
x=2.5
for (( c=0; c<$num_column; c++ ))
do # カラムのループ

 c1=`expr $c + 1`
 if [ $c1 -eq 1 ]
 then
  flg1="-K"
  flg2="-O -K"
 elif [ $c1 -eq $num_column ]
 then
  flg1="-O -K"
  flg2="-O"
 else
  flg1="-O -K"
  flg2="-O -K"
 fi

 width=`cat $cntl |\
        awk -v c="$c1" 'NR==1{print $c}'`
 gmt psbasemap -P -JX${width}c/${hight}c \
  -R0/${width}/-${num_row}/0 \
  -Bwesn \
  -Byg1 \
  -X${x}c \
  ${flg1} >> $psfile
 cat $temp |\
  awk -v c="$c1" '{printf "0.2 -%f %s\n",NR-0.1,$c}' |\
  gmt pstext -J -R -F+jBL+a0.0+f8p,GothicBBB-Medium-UniJIS-UTF8-H,0/0/0 ${flg2} >> $psfile

 x=$width

 
done # カラムのループ終端
rm $temp

cntlで指定するファイルは制御用のものでカラムの幅や行数・幅を指定しています。

例えば

1.0 10.0 5.0
20 0.4

だと、第一カラムが1.0cm、第二カラムが10.0cm、第三カラムが5.0cmの3列からなります。行数は20行で各行の高さが0.4cmであることを意味しています。

表形式に成型されているテキストファイルを標準入力から流し込んで使用します。

f:id:S_E_Hyphen:20180102111353j:plain

本日、早朝から報道されているニュースの見出しを表形式にしてみました。

2017-12-31

最大値の表示

二つのテーブルから、それぞれ最大値を検索します。それら最大値のうち、大きな方の値に着色して表示します。

テーブル original の最大値は変数max_oに、テーブルlowpass40の最大値は変数max_fにそれぞれ代入されています。

最後はawkで大小判定して、着色しています。

#!/bin/bash
called="xxxxx"
tab="T_YYYY"
for node in 8001 8002 8003 8004 8005 8006 8007 8008 8009 8010\
            8011 8012 8024 8036 8037 8038 8039 8040 8041 8042\
            8043 8044 8045 8046 8047 8057 8067 8068 8069 8070\
            8071 8072 8073 8074 8075 8076 8077 8078 8079 8080\
            8081 8082 8083 8084 8085 8086
do
max_o=`
mysql -N filter_design << _SQL_
 select max(abs(acc)),time from ${tab}
  where called="$called" and node=$node
  and proc="original" and abs(acc)=
  (select max(abs(acc)) from ${tab} where called="$called" and node=$node
  and proc="original");
_SQL_`

max_f=`
mysql -N filter_design << _SQL_
 select max(abs(acc)),time from ${tab}
  where called="$called" and node=$node
  and proc="lowpass40" and abs(acc)=
  (select max(abs(acc)) from ${tab} where called="$called" and node=$node
  and proc="lowpass40");
_SQL_`

echo $node $max_o $max_f |\
 awk '$2>$4{printf "%d \033[0;31m%.2f\033[0;39m (%.2fsec) %.2f (%.2fsec)\n",$1,$2,$3,$4,$5;next}
           {printf "%d %.2f (%.2fsec) \033[0;31m%.2f\033[0;39m (%.2fsec)\n",$1,$2,$3,$4,$5}'  
done

f:id:S_E_Hyphen:20171231095943p:plain