Scilabで連続ウェーブレット解析

id:ikeuchihiroki:20070222 で書いたように、連続ウェーブレット解析ができるフリーの環境というのはなかなか良いものがない。
Scilab上で動作するTOOLBOXでは以下のものがある。

Wavelab
Fraclab
Scilab Wavelet Toolbox

Scilab Wavelet Toolbox http://scwt.sourceforge.net/ は、小規模なウェーブレットツールボックスなのだが、地味に更新が続けられている。
連続ウェーブレット変換のコードはCで書かれていて、なかなか速い。

テストプログラム

fs=512;

x=[0:2047]/fs;
y=sin(2*%pi.*x.*(x/max(x)*64));

//        increase memory size
stacksize(10000000);

//        continous wavelet transform
scal=2^[1:1/32:6];
coef=cwt(y,scal,'cmor');

//        colormap obtaining
cmap=jetcolormap(512);
//        matrix normalization
n=wcodemat(abs(coef),512);

scf(0);
clf;
f=gcf();f.color_map=cmap;
xlabel Time[s]
ylabel Frequency[Hz]

getf scal2freq.sci

a=gca();

Matplot(n);

なかなかいい感じ。
でも、このままだと縦軸の意味がわからないので、スケールを周波数に変えて表示するようにしたい。
MATLABには、スケールを周波数に変換する関数が存在したので、同じような機能をする関数を書いてみる。

function [frq]=scal2freq(A,wname,fs)

  y=zeros(512,1);
  y($/2)=1;
  
  coef=cwt(y,A,wname);
  
  frq=[];
  for k=1:length(A)
    fftcoef=abs(fft(coef(k,:)));
    [m ind]=max(fftcoef,'c');
    
    tmp=ind-1;
    if tmp > length(y)/2
      tmp=length(y)-tmp;
    end
    frq=[frq;tmp/length(y)*fs];
  end
    
    
endfunction

Aがスケール、wnameがウェーブレットの名前、fsはサンプリング周波数。
理屈は簡単で、インパルスをウェーブレット変換して、各スケールの中心周波数を見ているだけだ。たぶん理屈はあっているはずだ。CWTの関数の制約上、スケールは長さ２以上のベクトルを入れないといけないので注意。

次に、縦軸横軸を調整する関数。

function coefs=cwt_plot(data,scal,wname,fs)
//plot settings
y_ticks_distance=16;

//        increase memory size
stacksize(10000000);

coef=cwt(data,scal,'cmor');

//        colormap obtaining
cmap=jetcolormap(512);
//        matrix normalization
n=wcodemat(abs(coef),512);

scf;
clf;
f=gcf();f.color_map=cmap;
xlabel Time[s]
ylabel Frequency[Hz]

getf scal2freq.sci

a=gca();

ay=a.y_ticks;
ay.locations=[0:y_ticks_distance:length(scal)]';
frq=scal2freq(scal($:-y_ticks_distance:1),'cmor',fs);
ay.labels=string(frq);
a.y_ticks=ay;


Matplot(n);

ax=a.x_ticks;
ax.labels=string(ax.locations/fs);
a.x_ticks=ax;
endfunction

dataは解析対象信号,scalはスケール,wnameはウェーブレット関数の名前,fsはサンプリング周波数。
ウェーブレット関数の名前は、cwt.cを見るとわかる。今のところ以下の関数が使える。

sinus,poisson,mexh,morl,DOG,
gaus1,gaus2,gaus3,gaus4,gaus5,gaus6,gaus7,gaus8,
cmor,shan,fbsp,cauchy,
cgau1,cgau2,cgau3,cgau4,cgau5,cgau6,cgau7,cgau8

使用例：

fs=512;

x=[0:2047]/fs;
y=sin(2*%pi.*x.*(x/max(x)*64));

cwt_plot(y,2^[1:1/32:8],'cmor',fs);

縦軸、横軸がサンプリング周波数に合わせて適正な表示になった。

意外にScilab Wavelet Toolboxが良い感じだ。