Microsoft Word - SCILAB入門.doc

Transcription

1 SCILAB 入門 琉球大学総合情報処理センター舟木慶一 1. はじめに 総合情報処理センターでは数値演算ソフト Matlab を第 3 実習室の Windows パソコンに 51 ライセンス Solaris サーバに 5 ライセンス導入している ( 新システムより Windows PC は 51 ユーザまで同時使用可能 ) Matlab は行列演算だけでなく 制御 信号処理 画像処理などの多彩な Toolbox を持ち グラフィック機能が充実しているために 世界中で広く活用されている しかし 高価なツールのために 個人で購入するのは困難である また 講義で Matlab を利用しても 宿題を第 3 実習室でしかできないという生徒からの苦情に悩まされている教員方も多いのではないだろうか そこで 今回 Matlab の最も良く出来たフリークローンである SCILAB( 注 : サイラブと発音されている ) のインストール方法 使い方について説明を行う SCILAB のメリットは Matlab とほぼ同じプログラム言語を持ち 行列演算機能や 関数群ならびに Matlab からのソース変換ツールを有し Matlab に近いかなり充実した PLOT 関数が用意されていることである 2. インストール方法 フランス国立コンピュータ科学 制御研究所の SCILAB のページより最新版 ( 原稿を書いている段階では 4.1) をダウンロードする Windows の場合 ダウンロードされたファイルをクリックすると インストールが行われます Mac の場合も 10.4 以降なら ZIP ファイルをダウンロードし解凍するだけです 3.SCILAB の起動と簡単な使い方 3.1 起動方法 ディスクトップに作成されてた アイコンをダブルクリックします 次のような窓が表示されるので OK をクリックする

2 次の初期画面が表示されます --> が入力待ちカーソルです 3.2 変数と基本演算 演算方法は Matlab と同じです SCILAB の全ての変数は複素数を要素とする行列です 実数を要素として持つベクトルやスカラー値はその特殊な形式です 習うより慣れろ です 早速 使ってみましょう 1スカラー値の代入 :A に 3 を代入します -->A=3 A = 3. 2 変数の表示 : 変数に格納されている値を示します 変数 ( リターン ) か disp( 変数 ) で変数の内容が表示されます -->A A = 3. -->disp(a) 3.

3 3 複素数 : 複素単位は %i です (Matlab は i または j) -->B=3+2*%i B = i -->conj(b) i ここで ans は答えであり 変数に演算結果を代入しない場合に 直前の演算結果が格納される変数です -->abs(b) 上のように conj( 変数 ) で共役 abs( 変数 ) で絶対値が算出できます A*B の算出は次のように実現できます -->A*B i 演算結果を別の変数 C に代入したければ 次のようにします -->C=A*B C = i 4ベクトル a(a とは異なる変数 ) に行ベクトル (1 2 3) を代入します -->a=[1 2 3] a = b(b とは異なる変数 ) に列ベクトル (4,5,6) の転置を代入します -->b=[4 5 6]' b = ここで は転置 ( 厳密には複素共役転置 ) を表します 行ベクトルは a のように列ベクトルは b のように表示されます 列ベクトルは次のように ; を使っても表記されます ; は次の行への改行を表します (, はと同じで次の要素になります )

4 -->b=[4;5;6] b = では b はどうなりますか? -->b' ベクトル a,b の要素は a(i) b(j) のように表します なお Matlab 同様 配列の添え字は 1 から始まります a(0) とするとエラーが出ます -->a(1),b(2),a(3) a(1)=1,b(2)=5,a(3)=3 で合っていますね 定義されていない値を表示しようとすると次のようにエラーが出ます -->a(4)!--error 21 invalid index では a や b に先程の複素数 B を代入しましょう -->B B = i -->a(2)=b, b(1)=b' a = i 3. b = i のようにできます a,b の複素共役転置は次のようになります

5 -->a',b' i i では 複素共役転置 (Helmite 変換 ) ではなく 単なる転置を実現したい場合は 次のように 変数. で表します -->a.',b.' i i 行列 行列はベクトルの拡張ですから簡単ですね 3 3 の行列 C の例を示します -->C=[1 2 3;4 5 6;7 8 9] C = >C(2,2) 5. -->C(2,2)=B C = i >C' i

6 -->C.' i C と C. は 複素共役転置と転置です 逆行列は inv( 行列 ) で求まります -->inv(c) i i i i i i i i i -->inv(c)*c, C*inv(C) D-17i D D D D-16i D-17i D-17i D D D-17i D-16i D D-17i 8.882D D D-16i 1. 誤差により 1 や 0 になっていない要素もありますが 両方とも単位行列になります b=cx の線型方程式は次のように解くことができます -->c=c\b c = i i i 答えが正しいか確かめて見ましょう -->C*c i

7 D-16i 6. -->b b = i 正しい答えが得られていることになります なお SCILAB にも QR 分解 (qr) Cholesky 分解 (lu chol) といった関数が用意されています 興味のある方は文献 [2] を参照願います 次の表記を頭に叩き込みましょう, またはスペース行の次の要素 ; 列の次の要素 [] 変数の代入値をベクトルまたはスカラーで表現 () 関数または 変数の要素 3.3 アレイ計算 A=[A(1),A(2),A(3),,A(10)] B=[B(1),B(2),B(3),,B(10)] と2つの列ベクトルが定義されているとします 各要素の加減乗除を実行しましょう 加算 A+B で [A(1)+B(1),A(2)+B(2),A(3)+B(3),...,A(10)+B(10)] を実行します 減算 A-B で [A(1)-B(1),A(2)-B(2),A(3)-B(3),...,A(10)-B(10)] を実行します 乗算 A.*B で [A(1)*B(1),A(2)*B(2),A(3)*B(3),...,A(10)*B(10)] を実行します 除算 A./B で [A(1)/B(1),A(2)/B(2),A(3)/B(3),...,A(10)/B(10)] を実行します.* や./ のことをアレイ演算と呼びます -->A=[ ] A = >B=[ ] B = >A+B

8 >A-B >A.*B >A./B さらに べき乗は次のように実行します -->A^ >A.^B D D+10 上の2 乗はすべて同じ値のためアレイ演算にする必要はありませんが 下の場合 アレ イ演算になります A と B の内積を算出する場合 行 列で実現できますので -->A*B' 385. となります なお A*B は行ベクトル同志の演算になるため 実行できません -->A*B!--error 10 inconsistent multiplication A/B は次のようになります -->A/B 1. A'*B とすると の行列になります

9 -->A'*B 注 )A*B A+B でも A,B のサイズまたは型異なる場合エラーが出ます A,B が行ベクトルなら A*B はエラーです A,B が列ベクトルなら A*B はエラーです A が行ベクトル B が列ベクトルなら A+B はエラーです A がスカラ B がベクトルなら A+B はエラーです 3.4 入出力方法 (1) テキストファイルの入出力スカラー値の入出力 A=read('file 名 ',1,1); write('file 名 ',A); A = read(%io(1),1,1); // 標準入力 print(%io(2),a); // 標準出力 A=を表示して標準入力したい場合は 次のようにします print(%io(2),'a= '); A=read(%io(1),1,1); (2) 1 次元ベクトルの入出力 A=read('file 名 ',1,n); write('file 名 ',A); nは列の数 (3) 行列の入出力 A=read('file 名 ',m,n);

10 write('file 名 ',A); mは行の数 n は列の数 (4) バイナリファイルの入出力 load(),save() でバイナリの入出力が行われます (4-1) ファイルを直接指定する方法 a=eye(2,2);b=ones(a); // a,b に値を設定 eye(n,n) は n n の単位行列 // ones(a) は a と同じサイズのすべての要素が1の行列を生成 save('val.dat',a,b); // val.dat に a,b を出力 clear a; clear b; // a,b をクリア load('val.dat','a','b');// val.dat から a,b を入力 a,b // a,bの出力 (4-2) ファイルポインタで指定する方法 fd=mopen('foo','wb'); // foo をいうファイルをファイルポインタ fd に割当 for k=1:4 x=k^2; save(fd,x,k); // x と k を fd にバイナリで出力 end mclose(fd); // 一度 fd を close fd=mopen('foo','rb'); // fd を再度ファイルポインタに割当 for i=1:4 load(fd,'x','k'); // x と k を入力 x,k; end mclose(fd); // fd を close fd=mopen('foo','mode') で mode は read 時は rb write 時は wb です (5)wav 形式音声ファイルの入出力 wav ファイルの LOAD x=loadwave('data/1.wav')

11 で data/1.wav という WAV 形式の音声データの音声部分が x(1),x(2),x(3),..., に入力される (6)wav ファイルの SAVE savewave('data/2.wav',x,[ サンプリング周波数 (Hz)]) で x(1),x(2),x(3),, の値が data/2.wav に WAV 形式で出力される ( 注意 )x=loadwave('./1.wav') で xには で正規化された 1 以下の値が入力され savewave('./2.wav',x) で 16 ビットの値に変換されて出力される その結果 1.wav とは LSB が異なる結果になる 3.5 変数の初期化 -->clear A で A の初期化 ( 値が0になるのではなく Aという変数の存在もなくなる ) -->clear ですべての変数が初期化されます A を表示させようとしても次のようになります -->A!--error 4 undefined variable : A 3.6 ディレクトリの移動 --> chdir ディレクトリ または File Change Directory でフォルダを指定すれば current directory を変更できる

12 ただし ディスクトップなどの漢字表記フォルダに移動することはできない 1 USB メモリ上にアルファベッドのフォルダを作成し そこで実行する のが無難ではあるが 漢字フォルダで実行したければ 2 漢字表記フォルダ内の SCILAB ファイルをダブルクリックして SCILAB を起動させる SciPad というエディタも起動するがファイルをエディットしない場合 ウィンドウを消せばいい File Get Current Directory でカレントディレクトリが確認できる 漢字名のディレクトリは文字化けする

13 4.PLOT の方法 (1)2 次元 PLOT の方法 plot() x=linspace(0, 5*%pi,100); // 0 から 5π の間の 100 サンプルを x に代入 plot(x,cos(x), r+,x,sin(x), b- ); // cos(x) を赤の+ sin(x) を青の実線で PLOT r+ b- などのグラフの線に関するフォーマットを表 1にまとめる Matlab とほぼ同じで ある 表記 線 マークの種類 表記 色 マーク - 実線 b 青 : 点線 g 緑 -. 1 点鎖線 r 赤 -- 破線 c シアン + 十字 m マゼンダ * 星 y 黄 s 四角 K 黒. 点 v 三角 ( 下向き ) o 丸 ^ 三角 ( 上向き ) x 印 < 三角 ( 左向き ) d ひし形 > 三角 ( 右向き ) p 5 角形 H 6 角形 (2) 複数のグラフを PLOT する方法 plot2d() xbasc(); x=linspace(0,10,100); y1=x+sin(x); y2=cos(x); y3=exp(-x).*y2; plot2d([x' x' x'],[y1' y2' y3'],[1 2 3]); 最後の [1 2 3] が各線の色を表す 色番号は 1: 黒 2: 青 3: 緑 4: 水色 5: 赤 6: 桃色 です plot2d() には次のバリエーションがあります plot2d1() 曲線 plot2d2() 折れ線 plot2d3() 縦線

14 plot2d4() 実際に PLOT させてみましょう plot2d1() は 下記のように X 軸 Y 軸を対数軸に設定することもできます 詳しくは [2][3] を参照されたい plot2d1("onn",x',[y1',y2',y3']); plot2d1("oln",x',[y1',y2',y3']); plot2d1("oll",x',[y1',y2',y3']); (3) 画面を複数設定する 複数のグラフを PLOT する場合 次のように xset( window, 窓番号 ) で窓を設定して PLOT する 例は cos(x) と sin(x) の PLOT である x=0:0.1:4*%pi xset('window',1); plot(x,cos(x)); xgrid(); xset('window',2); plot(x,sin(x)); xgrid(); 窓 1 に cos(x) 窓 2 に sin(x) が表示される xgrid( 色番号 ) でグリッドを表示する 窓は clf; または xbasc(); で初期化される 特定の窓を初期化する場合は xbasc( 窓番号 ) で初期化される (4) 複数の PLOT を 1 つの窓に収める方法 xset() で窓 1つに付き1つのグラフを表示できたが 1 つの窓で 2 つ以上のグラフを表示した方が比較にも画像化にも便利である そのためには plotframe() という関数を用いる 設定がやや面倒だが便利な関数である

15 plotframe(rect, tics, [A,B], ["C","D","E"], [X0,Y0,X1,Y1]); plot2d(x,y, 色の番号,"000"); で y(x) を描画します plotframe() は plot2d() の前に書きます 1) rect=[x 軸の最小値, Y 軸の最小値, X 軸の最大値, Y 軸の最大値 ]; で X 軸 Y 軸の最大最小値を指定します 2) tics=[x の小さい刻みの個数,X の大きな刻みの個数,Y の小さい刻みの個数,Y の大きな刻みの個数 ]; を指定します 大きな刻み間の小さな刻みが小さい刻みの数です 3) [A,B] でグリッドのあるなしと スケール表示ありなしを指定します グリッド スケール [%f,%f] なし あり [%t,%f] あり あり [%t,%t] あり なし [%f,%t] なし なし 4) ["C","D","E"] C で図のタイトル D,E が X,Y 軸の名前を指定します 5) [X0,Y0,X1,Y1] 描画する四角の大きさを指定します X0,Y0 が四角の左上の相対座標です ウィンドウの左上の角が (0,0) 右下の角が(1,1) になります X1,Y1 は四角の X,Y の大きさを表します 4 分割したければ 左上 : [0.0,0.0,0.5,0.5] 左下 : [0.0,0.5,0.5,0.5] 右上 : [0.5,0.0,0.5,0.5] 右下 : [0.5,0.5,0.5,0.5] のように設定します plot2d() の色番号は前述したとおりである 例を示す 4 種類の窓関数 ( 方形窓 三角窓 ハニング窓 ハミング窓 ) の時間波形を左側に そのスペクトルの振幅特性を右側に表示し 4 種類の窓関数を比較する図を PLOT する L=512; N=64; wtype=['re','tr','hn','hm']; // 窓関数のタイプ 順に方形窓 三角窓 ハニング窓 ハミング窓

16 xset('window',2); xbasc(2); for n=1:4 win=window(wtype(n),n); // 窓関数の指定 x=1:n; // plot(win); rect=[1,0,n,1]; tics=[2,2,2,2]; if (n==1) then plotframe(rect,tics,[%f,%f],['rec','time','amp.'],[0.0,0.00,0.5,0.25]); // 方形窓 end; if (n==2) then plotframe(rect,tics,[%f,%f],['tr','time','amp.'],[0.0,0.25,0.5,0.25]); // 三角窓 end; if (n==3) then plotframe(rect,tics,[%f,%f],['hanning','time','amp.'],[0.0,0.50,0.5,0.25]); // ハニング窓 end; if (n==4) then plotframe(rect,tics,[%f,%f],['hamming','time','amp.'],[0.0,0.75,0.5,0.25]); // ハミング窓 end; plot2d(x,win,2,"000"); end for n=1:4 win=window(wtype(n),n); for i=n+1:l; win(i)=0; end; z=abs(fft(win,-1)); // FFT の絶対値 z=z+(z==0)*%eps; // z=0 なら %eps という小さな値を代入 Log10(0) を算出しないため x=1:l/2; y=20*log10(z(x)); // 対数スペクトル m=max(y(x)); y=y-m; rect=[1,min(y),l/2,max(y)]; // rect=[1,-130,l/2,0]; tics=[2,2,2,2]; if (n==1) then

17 plotframe(rect,tics,[%f,%f],[ 'Rec','Freq.','Amp.[dB]'],[0.5,0.00,0.5,0.25]); // 方形窓 plot2d(x,y,3,"000"); end; if (n==2) then plotframe(rect,tics,[%f,%f],[ 'Tr','Freq.','Amp.[dB]'],[0.5,0.25,0.5,0.25]); // 三角窓 plot2d(x,y,4,"000"); end; if (n==3) then plotframe(rect,tics,[%f,%f],['hanning','freq.','amp.[db]'],[0.5,0.50,0.5,0.25]); // ハニング窓 plot2d(x,y,5,"000"); end; if (n==4) then plotframe(rect,tics,[%f,%f],['hamming','freq.','amp.[db]'],[0.5,0.75,0.5,0.25]); // ハミング窓 plot2d(x,y,6,"000"); end; end なお 3 次元 PLOT である plot3d() については [2][3] を参照されたい (3) グラフィック画面のファイルへの SAVE 方法グラフィックウィンドウ上で FILE->EXPORT で EXPORT ファイルタイプを指定します 下記では BMP を指定しています

18 窓枠もファイルに入れたい場合は ALT キー +PrintScreen キーを利用しましょう 5. プログラムの実行方法 SCILAB はファイルから実行できます プログラム名.sci というファイルを作成し SCILAB のプログラムを書けば 1 行目から順に実行します 実行方法は --> exec プログラム名.sci です プログラム名の文字はアルファベットならびに数字限定です 6.SCILAB プログラミング Tips 例で示してしまいましたが IF 文 FOR 文の使い方を示します (1)IF 文の使い方 if (A==1) then C=D; else C=E; end; if (A) then... ; else... ; end;

19 (A) が 0 以外の時,then 以下を実行.0 の時,else 以下を実行する. (2)FOR 文による LOOP for i= スタート値 : インクリメント値 : エンド値処理 ; end ( 例 :(x(i) y(i),i=1,2,...,10) for i=1:10; x(i)=y(i); end; i=1:10 と i=1:1:10 は同じです インクリメント値を省略すると 1 になります T ( 例 : 内積 < a, b >= a b ) A=0; for n=1:len A=A+a(n)*b(n); end a と b とが同じ長さの行ベクトルまたは列ベクトルの場合 a b または ab で実現できます (3) 配列のコピー A(A1:A2) を B(B1:B2) にコピーする場合 ( ただし A2-A1=B2-B1) j=b1; for i=a1:a2 B(j)=A(i); j=j+1; end でもいいですが 次のように 1 行で実現できます B(B1:B2)=A(A1:A2);

20 C=A(A1:$); で A(A1) から A の配列の最後の要素までが C に代入されます 7. おわりに Matlab のフリー版と言える SCILAB について説明を行った なお 本書をコピーアンドペーストしたい方も多いと思われるので 本書の HTML 版を [4] に公開する 本稿はマニュアルとしては貧弱であり 随時追加更新する予定である なお 詳しいマニュアルは [2][3] を参照されたい 参考文献 [1] SCILAB の HP( 日本語 ) [2] 広島大学大野先生の HP: 詳細なマニュアル Matlab と SCILAB のコマンドの違いも表で掲載されている [3] 筑波大桜井先生の HP: [4] 本稿の HP 版