スポンサーサイト

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

Scilabでデータをプロット(Scilabチュートリアル1)

今日は入門的なチュートリアルを書いて見ます。

Scilabでテキストデータを読み込んでグラフをプロットしてみましょう。

まずはテキストデータを用意してください。
ここでは、気象庁が配布している地震の観測波形データを利用してみます。
新潟県中越沖地震の小千谷市のデータをダウンロードしました。(H7165321.CSV)

Scilabを起動してください。

まず、カレントディレクトリをデータがあるフォルダに変更します。
ツールバーの「現在のディレクトリを変更」をクリックして変更できます。
20100515 01


CSVデータを読み込みます。下記のコマンドを入力してください。コピー&ペーストも可能です。

fp=mopen('H7165321.csv'); //ファイルオープン
txt=mgetl(fp,7); //ヘッダ部分読み込み
data=mfscanf(-1,fp,'%f,%f,%f'); //データ部分読み込み
mclose(fp); //ファイルクローズ



プロットします。

plot(data);


20100515 02
プロットできました。

このままですと、横軸はサンプル数になっていますが、これを時間に変えてみます。
サンプリング周波数は100Hzですので、時間系列を作ります。

fs=100; //サンプリングレート
time=[0:size(data,1)-1]/fs; //時間系列を作成


プロットします。先ほどのグラフはいったん閉じてから下記のコマンドを実行してください。

plot(time,data)


ついでに、タイトル、凡例、XY軸ラベルをつけてみましょう。

title('新潟県中越沖地震(小千谷市)'); //タイトル
legend('NS','EW','UD'); //凡例
xlabel('$Time[s]$'); //X軸ラベル
ylabel('$acceleration[cm/s^2]$'); //Y軸ラベル


20100515 04
$で囲むことでTeXの形式を利用することができますので、数式を綺麗に表示できます。

ちなみに1chだけ表示したい場合は下記のようにします。

plot(time,data(:,1));



一部分だけ表示したい場合、簡易的にはツールバーのズームエリアで拡大することができます。
いくつもグラフを作りたい場合など、決まった範囲のみを表示させたい場合は下記のようにします。

a=gca();
a.data_bounds(1:2,1)=[18;60];


20100515 05

今回はこの位にします。
地震データがなかなか興味深いのでこのデータを使って他にもチュートリアルを書こうと思いますので、
乞うご期待!


スポンサーサイト

Scilabで移動平均

Scilabで移動平均処理をしてみましょう。

移動平均をする専用の関数はありませんが、フィルター用の関数を用いれば簡単に実現できます。
移動平均は、方形のデジタルフィルターと等価です。

ここでは、Convol関数という畳み込み積分を実行する関数を利用します。
Forを用いるよりも高速に処理することができます。

下記は移動平均を行う関数です。

ret = moving_average(data,len)
 data: 元のデータ
 len: 移動平均長さ
 ret: 移動平均後のデータ

function [ret] = moving_average(data,len)
data=data(:);len_d=length(data);
tmp=[ones(len-1,1)*data(1);data;ones(len-1,1)*data($)];

h = ones(1,len)/len;
ret = convol(h,tmp);

side=(length(ret)-length(data))/2;
ret = ret(floor(side)+1:$-ceil(side));

endfunction



上記プログラムの特徴は
 ・Convol関数で高速な処理
 ・平均値が変化しない
 ・(lenが奇数の場合)ピーク位置が変化しない
 ・データ長が変化しない
です。
データ長が変化しないように、両端にデータを付け加えています。

テストプログラムは下記になります。

exec('moving_average.sci');
a=ones(11,1);a(6)=2;
b=moving_average(a,3); //3点移動平均
scf();
plot(a);
plot(b,'r');

scf();
t=linspace(0,1,1000);
sig=sin(2*%pi*1*t);
sig=sig+0.2*rand(1,length(t));
plot(sig)
plot(moving_average(sig,32),'r');
p=get("hdl");
p.children.thickness=5


test_moving_average.pngピーク位置が変化していないことがわかります。
test_moving_average2.pngノイズ除去に有効です。


全記事表示リンク

全ての記事を表示する

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