なんだか最近Scilabをよく使います。

• 大量のデータを読み込んで可視化したい
• グラフの見たい区間をインタラクティブに拡大したい
• 大量のデータをバッチで処理したい
• 気軽にFFT掛けたり特異値分解して解析したい
• しかもそれらを高速に処理したい

という要求を満たす良いツールが無いものかと思ったら、周りでScilabを使っている人がいたので、使ってみたらこれが便利！要はMatlabをパクったOSSなのですが、すっかり虜なのでした。

Scilab - Open source software for numerical computation

似たものだとOctaveも有名(Matlabと互換性が高い)。
https://www.gnu.org/software/octave

ひとつScilabの練習に遅延時間推定(DTE)をやってみましょう。
DTEとは、２つの信号の間の時間の遅れを推測する技術です。例えばレーダーやソナー等でも必要となる古典的技術です。

いろいろな手法があるので、興味があればサーベイ論文でも見てください。
http://www.control.isy.liu.se/research/reports/LicentiateThesis/Lic1061.pdf

ここでは簡易的に、信号 と に対して 遅れてノイズ を含む があるとして、 と から を求める問題とします。

つまり、

の を求めます。

問題をさらに簡単にするため、ノイズは無いものとして周波数領域に変換すると、 を の遅れを持つインパルス応答とすれば、

となります。これを について解けば、遅延が推定できます。


なんて簡単なのでしょう。
ただ、やってみると幾つか注意が必要。

• FFTのお約束として窓関数を掛けてデータを強引に周期的にする。しないと不要な周波数成分が発生する
• ノイズが増えると逆FFT後に求めたインパルス応答は綺麗に出ないので、各周波数ごとの位相差から遅延を求める(ゲインによって重み付けが必要と思われるが、今回は省略)

遅延を0.05s与えてそれがどの程度推定できるのか見てみます。

上から、

• x(t)
• y(t)
• X,Yのナイキスト線図
• X,Yのボード線図
• 推定したインパルス応答
• Y/Xの位相を遅延時間にしたもの

1%ほど誤差がありますが、わりと綺麗に0.05が求まってます。

この場合は問題なさそうです。(ただ、ノイズがランダムなので出方によっては誤差が10%くらいになることも)

0.064と30%近く真値とずれがあります。ただ、元々の信号の主成分が18Hzなので、そのあたりの推定値を使うようにすれば問題ないことが最後のグラフで分かります。今回の実装はゲインの低いノイズに左右されやすい周波数帯も含んでいて、それらを含めて単純な平均値を取っているのが問題です。

Scilabを使うことでアルゴリズムの試行錯誤から可視化まで、とても楽にできました。しかもこれがタダ。素晴らしい世の中になったものです。

samplingRate=1600
delayTime=0.05
NOISE=0.0

function x=originalSignal(t)
    F1=18
    F2=51
    x=sin(2*%pi*F1*t)+0.3*sin(2*%pi*F2*t)
    //x=(F1*t)-floor(F1*t)
    //x=squarewave(F1*t)*0.9
endfunction

function y=observedSignal(x)
    DELAY=floor(delayTime*samplingRate)
    y=cat(2,zeros(1,DELAY),x(0:$-DELAY-1))
    y=y + grand(1,N,'nor',0,NOISE)
endfunction

function plotSignal(t,x,tit)
    plot(t,x)
    title(tit)
    a=gca()
    a.grid=[2,2];
    a.x_label.text="time[s]"
    a.y_label.text="level"
endfunction

t=0:1/samplingRate:1
w=exp(-4*t)
N=size(t,'*')

scf(0)
clf

PLOT_ROW=6

// Make original signal
subplot(PLOT_ROW,1,1)
x=originalSignal(t)
plotSignal(t,x,"original signal")

// Make observed signal (includes delay+noise)
subplot(PLOT_ROW,1,2)
y=observedSignal(x)
plotSignal(t,y,"observed signal")

// FFT with window function
X=fft(x.*w)
Y=fft(y.*w)

M=samplingRate/2
f=samplingRate*(0:N-1)/N

// Draw nyquist
subplot(PLOT_ROW,1,3)
nyquist(f(1:M),[X(1:M);Y(1:M)],["original signal";"observed signal"],%F)
//black(f,[X(1:M);Y(1:M)],["original signal";"observed signal"])

// Draw bode
subplot(PLOT_ROW,1,4)
bode(f(2:M),[X(2:M);Y(2:M)]/*,["original signal";"observed signal"]*/)

// Calc estimated impulse response
subplot(PLOT_ROW,1,5)

Z=zeros(1,size(X,2))
Z(2:$)=Y(2:$)./X(2:$)
z=zeros(1,size(X,2))
z(2:$)=ifft(Z(2:$))

L=samplingRate*0.1
plotSignal(t(1:L),z(1:L),"impulse signal")

// Estimate the delay time
subplot(PLOT_ROW,1,6)
M=100 // Use phase data under M[Hz]
[Zdb,Zphi]=dbphi(Z(1:M))
d=-Zphi./f(1:M)/360
plot(f(1:M),d)
a=gca()
a.grid=[2,2];
a.x_label.text="freq[Hz]"
a.y_label.text="delay time[s]"
a.data_bounds=[0 0;M 0.1]

xstring(0,0,"estimated delay time="+string(mean(d(2:$)))+"[s]")