スポンサーサイト

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

Scilabで連続ウェーブレット変換-9 (実データを解析)

こんにちは。だいぶ間が空いてしまいましたが、今まで作ったプログラムを使って、実際の地震動データをウェーブレット解析してみます。
プログラムの修正点は、下記です。
 ・CSVファイル(地震動のファイル)を読み込んでいます
 ・ウェーブレット解析して表示する機能を関数化しています(設定値などはハードコーディングしている部分があります)
 ・時系列波形もあわせて表示するようにしています

//行列の画像表示
//Matplot_with_axes(data,rect)
// data : 画像化する行列
// rect : [xmin,ymin,xmax,ymax] X軸、Y軸の最小値、最大値を指定(省略可)
function Matplot_with_axes(data,rect)
if argn(2)<2 then rect=[1,1,size(data,2),size(data,1)],end
rectx=(rect(3)-rect(1))/(size(data,2)-1)/2;
recty=(rect(4)-rect(2))/(size(data,1)-1)/2;
rect=[rect(1)-rectx,rect(2)-recty,rect(3)+rectx,rect(4)+recty];
Matplot1(data($:-1:1,:),rect);

ca = gca();
ca.data_bounds = matrix(rect,2,2)';
ca.tight_limits = "on";
ca.axes_visible = ["on" "on" "off"];
ca.axes_reverse = ["off" "on" "off"];
// ca.x_location="top";
endfunction

//データのウェーブレット解析
//cwt_matplot(time,data,scales,wname)
// time : 時系列リスト
// data : 解析するデータ
// scales : スケール係数
function cwt_matplot(time,data,scales,wname)
//周波数の有効数字
sig_figs=3;

//連続ウェーブレット解析
coef=cwt(data,scales,wname);

//ウェーブレット係数を画像に変換
cmap_n=512;
cmap=jetcolormap(cmap_n);
n2=int(abs(coef)/max(abs(coef))*(cmap_n-1)+1);
f=gcf();f.color_map=cmap;
Matplot_with_axes(n2,[t(1),1,t($),size(scales,2)]);

xlabel("Time[s]");
ylabel("Frequency[Hz]");

//Y軸を変更する
scales_ticks=[1:2:size(scales,2)]'; //切りが良い数値のみ抜出し

a=gca(); //AXISを取得
yticks=a.y_ticks; //Y軸を取得
yticks.locations=scales_ticks; //位置を指定
freqs_=scal2frq(scales(scales_ticks)',wname,1/fs); //スケールを周波数に変換
//有効数字にまるめる
tmp_sig=10^(sig_figs-1);tmp=10^floor(log10(freqs_));freqs=(round(freqs_./tmp*tmp_sig)/tmp_sig).*tmp;
yticks.labels=string(freqs); //表示する内容を指定
a.y_ticks=yticks;
endfunction


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

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

//データ読みこみ
fp=mopen("L3118A41.csv");
header=mgetl(fp,7);
data=mfscanf(-1,fp,"%f,%f,%f");
mclose(fp);
x=data(:,1);

//時間軸生成
t=(0:size(x(:),1)-1)'/100;


// スケーリング列
scales = 2^[2:1/8:9];

//表示
scf(0);clf;
coef=cwt_matplot(t,x,scales,wname);

//元グラフの時系列表示
scf(1);clf;
plot(t,x);
ca=gca();
ca.zoom_box=[t(1),-max(abs(x)),t($),max(abs(x))];



m.jpg
n.jpg
研究ならばこの後、様々な考察を加えてデータを見ていくことになりますが、私は地震動についてはさっぱりわかりませんので、ここで終わります。

注意点として、このデータは36000点ありますが、解析に数十秒かかります。
データが大きすぎると、時間がかかりすぎるのと、表示するときもとても重いので、現実的ではありません。
Scilabはメモリ管理がいまいちなので、メモリ不足のエラーがでることもあります。(これはScilab6で修正されるらしいですが)
巨大なデータは、適当に分割するか、ダウンサンプルして解析してください。
スポンサーサイト

コメントの投稿

非公開コメント

coefの中身

有用な記事をありがとうございます。一点質問があるのですが、このコードを使って実データをウェーブレット変換すると、変換後のcoefのマトリクスの要素の値を見ることができません。coefの内容をファイルに書き出そうともしてみましたが、coefが定義されていないというエラーメッセージが返されます。何か解決策はありますでしょうか?
全記事表示リンク

全ての記事を表示する

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