スポンサーサイト

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

Scilabで連続ウェーブレット変換-5 (Matplotで表示を見やすく)

こんにちは。
今回は、表示を見やすくいじってみます。
データは前回と同じで、生成したデータを使っています。

//ウェーブレット関数
wname = 'cmor';

//サンプリング周波数
fs=100;

//時間軸生成
t = [0:1/fs:10];

//データ生成
x=sin((2*%pi.*(t)).*t);
x=x+x($:-1:1);

// スケーリング列
scales = 2^[1:1/10:6];
coef=cwt(x,scales,wname);

//ウェーブレット係数を画像に変換
cmap_n=512;
cmap=jetcolormap(cmap_n);
n2=int(abs(coef)/max(abs(coef))*(cmap_n-1)+1);

//表示
scf(0);clf;f=gcf();f.color_map=cmap;
ca=gca();
ca.axes_reverse(2)="on";
ca.tight_limits = "on";
Matplot(n2($:-1:1,:));

xlabel("Samples");
ylabel("Scale");

//scf();plot(t,x);


i.jpg

【前回からの変更点】
・ウェーブレット関数を「cmor」に変更
  cmorはComplex Morlet関数のことで、複素ウェーブレット関数です。
  複素ウェーブレット関数は、位相と振幅の情報を別々に得られるので、表示が見やすくなります。

・表示する方法をMatplotに変更
  前回までは、SWTの関数cwtplotを使っていましたが、今回はScilab標準の関数Matplotを使っています。
  そうすることで、表示をいじりやすくなります。
  
・Y軸が上下反転
  Y軸の上が高周波であるほうがわかりやすいので、上下反転するようにしています。

・カラーマップを「jet」に変更
  個人的にカラーマップはjetが好きです

・スケーリング列を2の冪数になるように生成
  スケーリング列を 2^[1:1/10:6] というように、2の冪数で生成しています。
  これは、好みの部分もありますが、ウェーブレット関数のオーバーラップが均等になり、
  冗長性が均一になる効果があります。
  実際の物理現象を理解するうえでも分かりやすい表示が得られやすいです。


7/6追記
上記のコードで表示した場合、縦軸にでている数値に意味はありません。
スケールと縦軸を一致させる例は次の記事でのせています。
スポンサーサイト

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

こんにちは。


今回は、もうちょっとわかりやすい例で紹介したいと思います。

wname = 'morl';
fs=100;

t = [0:1/fs:10];

x=sin((2*%pi.*(t)).*t);
x=x+x($:-1:1);

// Set scales
scales = [1:1:60];
coef=cwt(x,scales,wname);
cwtplot(coef,scales);

scf();plot(t,x);


h1.jpg
h2.jpg

2つの周波数が変化していく信号がMIXされていることがよくわかりますね。

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



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

Scilabで連続ウェーブレット変換-2 (ファイルの読み込みとプロット)

こんにちは。前回に引き続き、Scilabで連続ウェーブレット変換をする方法を紹介していきます。


今回は、ファイルの読み込みとプロットをしてみましょう。

以前の記事
Scilabでデータをプロット改(Scilabチュートリアル2)
を参考にします。
地震動データ(CSVファイル)をダウンロードしておいてください。

Scinoteを起動して、以下のプログラムを入力し、データファイルがある場所に保存してください。

fname="L3118A41.csv"; //ファイル名を指定
fs=100; //サンプリングレートを設定

//ファイルを読み込む
fp=mopen(fname);
comments=mgetl(fp,7); //ヘッダ部分読み込み
data=mfscanf(-1,fp,"%f,%f,%f");
mclose(fp);

time=[0:size(data,1)-1]/fs; //時間系列を作成

//グラフを作成
scf(0);//Figure0を準備
clf; //Figureをクリアする
plot(time,data(:,1)) //1列目をplotする

xlabel("時間 s"); //X軸ラベル
ylabel("加速度 gal"); //Y軸ラベル


プロットした地震動

次回はいよいよウェーブレット変換を行います。

Scilabで連続ウェーブレット変換-1 (準備)

これから、何回かにわたって、Scilabを使って、連続ウェーブレット変換で信号解析する方法を紹介していきます。

発端はYahoo!知恵袋です。Yahoo!知恵袋では時々、Scilabに関する質問が出るのでできるだけ回答しているのですが、先日こういう質問がありました。

卒研でscilabを用いてウェーブレット解析をしたいのですが、やり方がわかりません...

Toolboxを紹介したまでは良かったものの、それ以降の回答をしようとプログラムを書いていたら、回答が締め切られてしまいました。せっかく書いたので、ブログでチュートリアル風にして残しておこうと思います。


まずは、知恵袋でも書いた準備の部分です。


【Scilab Wavelet ToolboxとSIVPをインストール】
・Scilabを起動
・メニューの「アプリケーション→モジュール管理」をクリック
・「Signal Processing→Scilab Wavelet Toolbox」をクリックし、「インストール」をクリック
・あわせて画像処理用ライブラリが必要なので「Image Processing→Scilab Image and Video Processing toolbox」をクリックし、「インストール」をクリック
・インストール完了後、Scilab再起動



これで、Toolboxがインストールされました。


【デモをみてみましょう】
・メニューの「その他→Scilab デモ」をクリックし、「Toolbox swt」をクリック

多くのデモがありますが、continuous wavelet transform が連続ウェーブレット変換です。


b.jpg

c.jpg

d.jpg



今回はここまで。続きはぼちぼち書きます。




全記事表示リンク

全ての記事を表示する

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