スポンサーサイト

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

Scilab5.5.1は再インストールが必要


Windows users, reinstall Scilab 5.5.1
10月15日までにScilab5.5.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で修正されるらしいですが)
巨大なデータは、適当に分割するか、ダウンサンプルして解析してください。

Scilab5.5.1リリース

Scilab5.5.1がリリースされました。



Scilab 5.5.1 Release / News / Community / Home - Scilab


バグフィックスが中心のようです。

Scilabで連続ウェーブレット変換-8 (縦軸を周波数,横軸を時間に)

こんにちは。
今回は、周波数の軸をもう少し見やすくして、さらに横軸を時間軸に変更します。
前回の表示ですと、周波数の軸が、無駄に桁数が長く表示されていました。そこで、有効数字で丸めてみます。
また、横軸がSample数のままでしたので、時間軸に変更します。

//関数定義
//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




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

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

//周波数の有効数字
sig_figs=3;

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

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

// スケーリング列
scales = 2^[1:1/8: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;
Matplot_with_axes(n2,[t(1),1,t($),size(scales,2)]);

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

//Y軸を変更する
//scales_ticks=find(fix(scales)==scales)'; //切りが良い数値のみ抜出し
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;


l2.jpg



見栄えの問題だけですが、見やすくなりました。

以前に Scilabで行列の内容を画像化する で紹介した、Matplot1を使って軸の最大・最小値を指定できる自作関数を使っています。
横軸を時間で表示できるようになりました。
縦軸は、有効数字3桁でまるめる処理を行っています。まるめたら見やすくなったので、表示を密にしてみました。

Scilabで連続ウェーブレット変換-7 (縦軸を周波数にする)

こんにちは。
今回は、前回に引き続いて、さらに縦軸を修正して見ます。

縦軸の表示をスケールから周波数に変えて見ましょう。

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

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

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

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

// スケーリング列
scales = 2^[1:1/8: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");

//Y軸を変更する
scales_ticks=find(fix(scales)==scales)'; //切りが良い数値のみ抜出し

a=gca(); //AXISを取得
yticks=a.y_ticks; //Y軸を取得
yticks.locations=scales_ticks; //位置を指定
yticks.labels=string(scal2frq(scales(scales_ticks)',wname,1/fs)); //表示する内容を指定
a.y_ticks=yticks;


k.jpg

プログラムは前回から1行変わっているだけです。scal2frqを足しました。
ウェーブレットの周波数が中途半端なので、半端な表示になりますが、わかりやすくなったのはないでしょうか。
全記事表示リンク

全ての記事を表示する

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