総合情報処理センター 舟木慶一
2004/1/8 作成開始
2004/5/2改訂
2004/5/8 plotframe()追加
2004/5/28 配列のコピーとARスペクトルを追加
SCILABは、 フランス国立 コンピュータ科学・制御研究所(INRIA)で開発されたフリーの高機能数値演算プログラムです。Matlabライクなソフトウェアであり、Matlabソフトの変換ツールも提供されています。
(1)インストール方法
フランス国立 コンピュータ科学・制御研究所のSCILABのページより最新版をダウンロードする。
(1-1)WINDOWS
Windows版バイナリーをダウンロードし、解凍する(簡単)。
(1-2)LINUX
Linux版ソースをダウンロードし、コンパイルする。
# gtar -xzvf scilab-2.7.src.tar.gz # ./configure --with-gnu --without-tk # gmake all |
(1-3)MAC OS 10.3
準備がやや厄介です。マシンのスピードにも依存しますが、LINUXとかと異なり、コンパイルに時間もかかります。あまり勧めませんが、どうしてもMACで使いたい方は次のページを参考にインストールしてみましょう
SCILABをMAC OS 10.3にインストール(舟木作成)
MusinSystemのページ。
(2)起動方法
(2-1)UNIX(DEP,CC)にログインして利用。
LINUX、MAC OS Xからは、
(DEP)
% ssh -X 133.13.4.10 |
% ssh -X 133.13.2.10 |
WINDOWSからは、ASTEC-XなどのXエミュレーションソフトを利用して、LOGINする。
WINDOWSはアイコンをクリックするだけで起動します。
(2-2)SCILABの起動
% tcsh |
で次のような画面が立ち上がります。 (DEPからASTEC-XでWINDOWS XPに表示)
(3)実行方法
WINDOWS版もUNIX版もほぼ同じです。
(3-1)コマンドライン入力
--> chdir('ディレクトリ') |
でディレクトリ移動し、
--> exec('scilabファイル名') |
で実行できます。
(3-2)GUI入力
(3-2-1)WINDOWS版
『file』→『change directory』
でファイルの存在するフォルダに移ります。
『file』 → 『exec』 でSCILABファイルを選択して、クリックして実行します。
(注意)データ入力をする場合、そのファイルが存在するディレクトリに移動しないと、データ入力できません。
(3-2-2)Solaris版(DEP,CC)
コマンドラインで、ディレクトリを変更します。
--> chdir('ディレクトリ'); |
『File』 → 『File Operations』で次の画面が立ち上がります。
上のようにAlternativesから実行したいファイル(ここではfft1.sci)をクリックします。
をクリックすれば、選択したファイルが実行されます。
(4)HELPの利用
--> help |
または、Helpをクリックすれば、HELP画面が立ち上がります。
--> help lev |
で次のように、lev( )という関数の説明が表示されます。
(5)SCILAB言語の文法
○変数の入力
L=100;
M=100;
で変数LとMに100が代入されます。
○配列変数の入力
A(1)=100;
A(100)=100;
でA(1)とA(100)に100が代入されます。
(注意)
配列はMATLAB同様A(1)から定義されています。A(0)はありません。注意しましょう。
○テキストファイルの入出力
スカラー値の入出力
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);
1次元ベクトルの入出力
A=read('file名',1,n);
write('file名',A);
nは列
マトリクスの入出力
A=read('file名',m,n);
write('file名',A);
mは行、nは列
○バイナリファイルの入出力
data_nameを入力し,data_name.lpcというファイルからデータを入力する.
clear;
mtlb_hold('off');
data_name = input(' data name: ','s');
LPC_name = mtlb_sprintf('%s.lpc',data_name);
len = 257;
[fid,%v] = mopen(LPC_name,'r',0)
if %v<0 then fid = -1;end
LPC = mtlb_fread(fid,spc_len,'double');
status = mclose(fid);
mtlb_plot(LPC,'red');
mtlb_hold('on');
○wav形式音声ファイルの入出力
wavファイルのLOAD
x=loadwave('data/1.wav') |
でdata/1.wavというWAV形式の音声データの音声部分がx(1),x(2),x(3),...,に入力される。
wavファイルのSAVE
savewave('data/2.wav',x,[サンプリング周波数(Hz)]) |
でx(1),x(2),x(3),…,の値がdata/2.wavにWAV形式で出力される。
(注意)x=loadwave('./1.wav')で xには32768で正規化された1以下の値が入力され、savewave('./2.wav',x)で16ビットの値に変換されて出力される。その結果、1.wavとはLSBが違ってくる。
□2次元PLOTの方法
plot(x,y)
plot(y)
2つの関数をPLOTする方法
plot([sin(x),cos(x)])
plot
plot2d
plot2d1 曲線
plot2d2 折れ線
plot2d3 縦線
plot2d4 ○
(使い方)
plot2d1("onn",x,[y1(x),y2(x),y3(x)]); |
3次元PLOT plot3d()
軸のスケール補正:調査中です。
(Matlab)
(Scilab)
画面を複数設定する。
xset('window',1); |
PLOT関数でPLOT
xset('window',2); |
PLOT関数でPLOT
窓の初期化
xbasc() |
特定の窓を初期化する場合は、
xbasc(窓番号) |
で初期化される。
□plotframe()で複数の描画を1つの窓に収める方法。
plotframe(rect, tics, [A,B], ["C","D","E"],
[X0,Y0,X1,Y1]); plot2d(x,y,色の番号,"000"); |
グリッド |
スケール |
|
[%f,%f] |
なし |
あり |
[%t,%f] |
あり |
あり |
[%t,%t] |
あり |
なし |
[%f,%t] |
なし |
なし |
□グラフィック画面のファイルへのSAVE方法
(Matlab) % print -deps 111.eps により,111.epsというEPSファイルが出力されます。
(SCILAB)
グラフィックウィンドウ上で、FILE->EXPORT
でEXPORTファイルタイプを指定します。
(Windows)
(Solaris)
Export typeは、
(1)Postscript-Latexでeps、
(2)Xfigでfig、
(3)Gifでgifファイルが出力されます。なお、(3)のGIFでは窓の枠がなくなります。窓枠を付けて出力したければ、WINDOWSのALT+PrintScreenにより出力しましょう。
WINDOWSでは『OK』をクリックし、出力ファイル名を指定しましょう。
Solarisでは、FileNameのボックスに拡張子なしの名前を入力しOkをクリックします。
IF文の使い方
if (A==1) then C=D; else C=E; end; |
if (A) then .... ; else .... ; end; |
(A)は0以外の時,then以下を実行.0の時,else以下を実行する.
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になります。
(例:MAフィルタリング処理)
for n=1:len r(n)=x(start+n); for i=1:order r(n)=r(n)+b(i)*x(start+n-i); end end |
ar(i): フィルタ係数
x(n): 入力信号
res(n): フィルタ出力
order: 次数
len: 信号長
start:入力信号の開始サンプル
x(start+1)からx(start+len)をMAフィルタ(MAフィルタ係数:b(1),b(2),...,b(order))でフィルタリングし、出力をr(1),r(2),...,r(len)に出力する。
配列のコピー
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); |
(6)プログラムの例
自己相関関数の算出とLPC分析
// 自己相関関数の算出 j=1; for i=0:order r(j)=0; for n=1:len-i r(j)=r(j)+x(n)*x(n+i); end j=j+1; end [ar,sigma2,rc]=lev(r); // LPC分析 lev() |
(入力)
x(1),x(2),...,x(len): 信号
len: 分析長
order: 次数
(出力)
r(1),r(2),...,r(order+1): 自己相関関数 0次が1次になっていることに注意。
ar(1),ar(2),...,ar(order): LPC係数
sigma2: 予測誤差のエネルギー
r(0)+_sum_{i=1}^order(r(i)*a(i))
rc(1),rc(2),...,rc(order): 反射(PARCOR)係数
対数スペクトルの算出
s=20*log10(abs(fft(x,-1))); |
x(1),x(2),…x(n)(nはデータ数=FFT点数)をfft(x,-1)によりFFTし、abs()により絶対値を算出し、20*log10()でdB値に変換し、s(1),s(2),...,s(n)に対数スペクトルの値が代入される。
ただし、配列xの長さは2のべき乗になるFFT長であり、信号の長さよりも長くなり、信号長以上では0詰めする。
xがMA係数の場合は、上記で算出できる。
xがAR係数の場合は、伝達特性の逆数に対数を掛けるため、マイナスをつける。s=-20*log10(abs(fft(x,-1)));
ARフィルタの対数スペクトルの算出
a(1)=1.0; a(2:order+1)=ar(1:order); // ar(1),...,ar(order)をa(2),...,a(order+1)にコピーする。 s=-20*log10(abs(fft(a,-1)))+10*log10(sigma2); |
伝達特性を算出し、零点を求める。
H=poly(a(1:order+1),'z','coeff'); p=roots(H); |
a(1),a(2),...,a(order+1)のフィルタ係数からなるorder次のZの多項式をpoly(a(1:order+1),'z','coeff')により算出しHに代入。P=roots(H)により、Hの根を算出し、p(1),p(2),...,p(order)に代入。
極周波数から伝達特性の分母の多項式の係数を算出。
z=poly(0,'z');// 多項式の変数をzに指定。 H=1;// 伝達特性の初期値を1に設定。 for i=1:order/2 H=H*poly([A(i) B(i)],'z','root'); // 多項式(伝達関数)の算出 end c=coeff(H); // 伝達関数Hから係数を抽出しc(1),c(2),...,c(order)に代入する。 |
A(1),A(2),...,A(order/2)は極周波数。
B(1),B(2),...,B(order/2)はA(i)の共役複素数。
poly([A(i) B(i)],'z','root'); により、(z-A(i))(z-B(i))を算出。
(注意)Zは1/Zでないのでcの係数は順番が逆になる。
(7)リンク集
フランス国立 コンピュータ科学・制御研究所のSCILABのページ
広島大学大野先生によるSCILABのページ
筑波大学桜井先生によるSCILABのページ