SCILABのコンボリューション

scilabでコンボリューションを実施するときは、c=convol(b,a)を使います。

しかし信号処理の場合、周期Tをもって無限に繰り返し循環する「継続時間が無限大」の関数とみなして取り扱いたいことがあります。具体的には波の終端に始点がつながっているとします。

f:id:S_E_Hyphen:20180128122318j:plain

このようなとき、

[c,e]=convol(b,a); [d,e]=convol(b,a,e);

とすると、2回目のconvolの結果は尻尾と頭がつながっているとしてオーバーラップ加算を実施してくれます。

図中の黒線がオリジナルの例題波です。これにconvol を用いて [1/3 1/3 1/3] という3要素分の移動平均を実施した例です。

青線ではオーバーラップ加算を考慮していません。始点の青丸は1と2の平均ですから、ほとんどゼロです。一方赤線はオーバーラップ加算を考慮しているため、1の前に16の要素があります。ですから、その平均値は12.33となっています。

2から15に関しては、前後に最低1個以上要素がありますので、オーバーラップ加算を考慮しても、しなくても同じ値となります。

ハニングやハミングウィンドウによるスペクトルの平滑化を行う場合、このオーバーラップ加算は役に立つものと思います。

#!/bin/bash
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 10p,GothicBBB-Medium-UniJIS-UTF8-H,black

scilab -nw << _SCILAB_
 a=[38 -33 -19 -10 1 -8 -20 10 -1 4 11 -1 -7 -2 5 32 ];
 b=[1/3 1/3 1/3];
 rot=(size(b,'*')-1)/2;
 [c,e0]=convol(b,a);
 [d,e1]=convol(b,a,e0);
 c=[c(rot+1:16) c(1:rot)];
 d=[d(rot+1:16) d(1:rot)];

 fd=mopen('smoothed.txt','w');
 for i=1:16,
  mfprintf(fd,'%d %f %f %f\n',i,a(i),c(i),d(i));
 end,
 mclose(fd);

_SCILAB_

cat smoothed.txt |\
 awk '{print $1,$2}' |\
 gmt psxy -JX20c/4c -R0/17/-40/40 -BnSeW+t"b=3⊿t" -Bxa1+l'pt' \
 -Bya20g100 -W1p,black -K > convol_test.ps

cat smoothed.txt |\
 awk '{print $1,$3}' |\
 gmt psxy -J -R -B -W1p,blue -O -K >> convol_test.ps

cat smoothed.txt |\
 awk '{print $1,$4}' |\
 gmt psxy -J -R -B -W2p,red -O >> convol_test.ps

 なおコンボリューションの性質上、38*1/3+(-33)*1/3+(-19)*1/3 の結果が第一要素に入ってしまいますが、これを第二要素として取り扱いたいという要求が有ります。このためc=[c(rot+1:16) c(1:rot)];という、やや巧妙な(小狡い?)ことを行っています。変数rotはbの要素数を考慮して1が入るようになっています。