スポンサーサイト

上記の広告は1ヶ月以上更新のないブログに表示されています。
新しい記事を書く事で広告が消せます。

Scilabで連続ウェーブレット変換-3 (cwtplotで表示)

こんにちは。
今回は、「Scilab Wavelet Toolbox」に付属している関数「cwtplot」を使って、
連続ウェーブレット解析をしてみましょう。

まずは、仮想的につくったデータで行います。

以下のプログラムはcwtplotのヘルプに書かれているものです。
一番下の行は私が追加しました。

wname = 'morl';
A = 0; B = 64; P = 500;
// Compute the sampling period and the sampled function,
// and the true frequencies.
t = linspace(A,B,P);
delta = (B-A)/(P-1);
tab_OMEGA = [5,2,1];
tab_FREQ = tab_OMEGA/(2*%pi);
tab_COEFS = [5,3,2];
x = zeros(1,P);
for k = 1:3;
x = x+tab_COEFS(k)*sin(tab_OMEGA(k)*t);
end
// Set scales
scales = [1:1:60];
coef=cwt(x,scales,wname);
cwtplot(coef,scales);

scf();plot(t,x);


e.jpgf.jpg

1つめの図がウェーブレット変換の図で、2つめが元データです。
それっぽく出ているのですが、で何がわかったのか・・・といわれると難しいですよね。

このプログラムは、3つのSIN波が加算されています。
tab_OMEGAが、SIN波の周波数を指定していて、
tab_COEFSが、振幅を指定しています。

例えば、1つのSIN波で計算した場合はどうなるのでしょうか。
tab_COEFSを定義している行を下記の3パターンに書き換えて、実行してみましょう。

tab_COEFS = [5,0,0];
tab_COEFS = [0,3,0];
tab_COEFS = [0,0,2];


g1.jpgg2.jpgg3.jpg
この図をみるとわかりましたでしょうか?先ほどの結果は上記の3つの結果が合算されたものだったわけです。
縦軸は、下に行くほど、周波数が高く、上に行くほど周波数が低いことをあらわします。


【周波数を読み取るには?】
Scaleが意味するところは、Wavelet関数が、時間軸方法に、何倍に拡大されたかということです。
Scale=1ならば、Wavelet関数は定義の通りの大きさで使用され、Scale=60ならば、60倍に拡大されて使用されます。Wavelet関数が時間軸方向に大きければ、周波数は低くなります。

Scaleと周波数との関係を調べる関数「scal2frq」というのが準備されています。
「scal2frq(1,'morl',delta)」と入力すると、Scale=1の時の周波数を返してくれます。deltaは上記のプログラムで計算されたサンプリング周期です。

下記を入力し、縦軸があらわす周波数を調べてみましょう。

for sc=[5:5:60]
fq=scal2frq(sc,"morl",delta);
printf("%d , %f\n",sc,fq);
end


下記の結果が得られました。左の列がスケール、右の列が周波数です。
5 , 1.164963
10 , 0.582481
15 , 0.388321
20 , 0.291241
25 , 0.232993
30 , 0.194160
35 , 0.166423
40 , 0.145620
45 , 0.129440
50 , 0.116496
55 , 0.105906
60 , 0.097080



今日のところはこれくらいで。
正直、この関数をそのまま使って、実際に収録された信号の解析は難しいと思います。
また、細かいノウハウもありますので、色々と紹介していきます。
スポンサーサイト

コメントの投稿

非公開コメント

全記事表示リンク

全ての記事を表示する

最新記事
最新コメント
最新トラックバック
カテゴリ
検索フォーム
RSSリンクの表示
リンク
上記広告は1ヶ月以上更新のないブログに表示されています。新しい記事を書くことで広告を消せます。