gnuplotでヒストグラムのグラフを作成する

自分の研究アーカイブ用記事として残しておきます．
（間違っている可能性もあるので責任は持ちませぬ．．．）

脳波信号（安静時 vs 集中時）を時間周波数解析（時間FFT）して，時間波形から周波数波形に直した時に出てきたピークの値をヒストグラムにして表示します．
つまり，安静時の脳波と集中時の脳波の周波数のピークの推移を見たい．
的な考えで良いと思います．

時間FFTの窓幅は約1秒の矩形窓です．
時間FFT後のピークの値はすでにテキストか何かに落とせているものとします．

以前の記事でプログラムと，やりたいことは公開してます．

Duffing振動子のモデルパラメータを最小二乗法で実験的に同定してみる

gnuplotのコマンドとして以下

軸の設定

set xrange [0:240]
set xlabel "Time (s)"
set xlabel font "Times New Roman, 20"
set xtics font "Times New Roman, 20"

set yrange [0:30]
set ylabel "Frequency (Hz)"
set ylabel font "Times New Roman, 20"
set ytics font "Times New Roman, 20"

set key font "Times New Roman, 20"

set xrange [0:240]

set xlabel "Time (s)"

set xlabel font "Times New Roman, 20"

set xtics font "Times New Roman, 20"

set yrange [0:30]

set ylabel "Frequency (Hz)"

set ylabel font "Times New Roman, 20"

set ytics font "Times New Roman, 20"

set key font "Times New Roman, 20"

時間窓ごとの周波数ピークの描写

plot "relax.txt" u 1:($6/(2*pi))  w lp t "relax"lc rgb "green", "stress.txt" u 1:($6/(2*pi))  w lp t "stress" lc rgb "red"

1	plot "relax.txt" u 1:($6/(2pi)) w lp t "relax"lc rgb "green", "stress.txt" u 1:($6/(2pi)) w lp t "stress" lc rgb "red"

おいらのtxtデータでは，1列目に時間，6列目に固有角振動数(2πf)が入っているので，６列目を2πで割ってグラフに書いてます．

set term jpeg
set out "result.jpg"
replot
set term aqua

set term jpeg

set out "result.jpg"

replot

set term aqua

jpegで画像を保存する時のコマンド．
term はterminalの略．
画像の保存が終わったら元のターミナル(aqua)も戻しておきます．

今回は時系列データの60秒から180秒の間で負荷を与えたので，こちらでピーク周波数が高い結果になっています．

実はこれは第1ピークですが，第3ピークまで同様に取得してます．
次に，こいつらをヒストグラムにしていきます．

各ピークのヒストグラムの描写

ヒストグラムの幅を指定して

binwidth=0.25*0.5
bin(x,width)=width*floor(x/width)+width/2.0

1 2	binwidth=0.250.5 bin(x,width)=widthfloor(x/width)+width/2.0

relax.txtを(relax1.txtとrelax2.txt)に分けました．

set multiplot layout 3,3
set xrange [0:30]
set xlabel "frequency (Hz)"
set ylabel "numbers"
set yrange [0:30]
plot 'relax1.txt' u (bin(($6/(2*pi)),binwidth)):(1.0) smooth freq with boxes t "omega1(a-relax)" lc rgb "green"
plot 'stress.txt' u (bin(($6/(2*pi)),binwidth)):(1.0) smooth freq with boxes t "omega1(a-stress)" lc rgb "red"
plot 'relax2.txt' u (bin(($6/(2*pi)),binwidth)):(1.0) smooth freq with boxes t "omega1(a-relax)" lc rgb "green"
plot 'relax1.txt' u (bin(($9/(2*pi)),binwidth)):(1.0) smooth freq with boxes t "omega2(a-relax)" lc rgb "green"
plot 'stress.txt' u (bin(($9/(2*pi)),binwidth)):(1.0) smooth freq with boxes t "omega2(a-stress)" lc rgb "red"
plot 'relax2.txt' u (bin(($9/(2*pi)),binwidth)):(1.0) smooth freq with boxes t "omega2(a-relax)" lc rgb "green"
plot 'relax1.txt' u (bin(($12/(2*pi)),binwidth)):(1.0) smooth freq with boxes t "omega3(a-relax)" lc rgb "green"
plot 'stress.txt' u (bin(($12/(2*pi)),binwidth)):(1.0) smooth freq with boxes t "omega3(a-stress)" lc rgb "red"
plot 'relax2.txt' u (bin(($12/(2*pi)),binwidth)):(1.0) smooth freq with boxes t "omega3(a-relax)" lc rgb "green"
unset multiplot