Octave入門

Octave（正確にはGNU Octave）は行列・ベクトルの数値演算を行なうことのできる，コマンドラインインタフェースを備えたフリーの汎用演算ソフトです．本家webサイトは http://www.octave.org/ です．

Octave自体に組み込まれている関数やコマンドは基本的なものに限られていますが，多くの関数やコマンドを世界中の人が開発し，公開するプロジェクトがあります（ http://octave.sourceforge.net/ ）．その成果はOctave-forgeというパッケージとして公開されています（関数一覧は http://octave.sourceforge.net/index/index.html ）．研究室のマシンにインストールされているのもOctave-forgeパッケージです．

基本操作

起動と終了

起動はWindowsPCならスタートメニューのGNU Octave 2.1.50 を選びます．MacOS XやlinuxなどPC-unix系のシステムなら，ターミナルを開いてoctaveとタイプすればいいです．起動するとターミナルが開き，コマンドを受け付ける状態になります．終了はexitとタイプします．

ヘルプ

Octave上で関数やコマンドのヘルプを表示するには，

>>help 関数名

とします（ここで>>はプロンプトの例）．見にくい場合は，前述のwebサイト上のmanualページを参照してもいいです．

また，

>>help -i キーワード

とすると，キーワード（関数名など何でも良い）に関連した説明が表示されます．

カレント・ディレクトリ

Octaveを起動して開くウィンドウにコマンドを打ち込んでいけば利用はできますが，結果をファイルとして残したり，実験データなとを読み込む場合にはディレクトリを意識する必要があります．今現在自分が作業しているディレクトリをカレント・ディレクトリと呼びます．unix系のOSならすでに各自のホーム・ディレクトリで作業し始めているので，問題は少ないですが，Windowsの場合は注意が必要です．

カレント・ディレクトリを確認するコマンドは，pwdで，Windows上でOctaveを起動した直後は，

>> pwd
/octave_files

となっていると思います．ただし，Cドライブ直下のoctave_filesディレクトリというわけではなく，GNU Octaveをインストールした場所（普通はC:/Program Files/GNU Octave 2.2.50/）の下にあります．カレント・ディレクトリをこのままで作業すると，各ユーザのファイルが混在しますし，PCのローカルなHDD上なのでログインするPCが変わると作業が続けられません．そこで，ネットワーク上の各自のホーム・ディレクトリで作業する必要があります．カレント・ディレクトリを移動するにはcdコマンドを利用しますが，cygwim上でoctaveを利用している場合，ドライブをまたいで移動するときは頭に/cygdrive/ドライブ名/をつける必要があります．研究室のWindowsPCにログインするとネットワークドライブはH:になるので，その下のDocumentsフォルダをカレント・ディレクトリとするには，

>> cd /cygdrive/h/Documents

とすることになります．

基本数値・算術コマンド

数値，ベクトル・行列の表現の仕方

Octave上でスカラー定数は，

実数：1.6e-31
複素数：1.2+0.8*i (1.2+0.8*jと入力しても良い)

のように入力し，表示もされます．

ベクトル・行列は，

行ベクトル：[1 2 3]
列ベクトル：[1;2;3]
行列：[1 2 3;4 5 6;7 8 9]

のように入力します．要素は実数や複素数，後述の変数（数式）でも構いません．また，要素の区切りはスペース以外にカンマ（, ）でも大丈夫です．行列はベクトルをセミコロン（; ）で改行しながら入力するというイメージです．

変数

Octave上では，実数，複素数，行列，ベクトルの区別なく変数に格納することができます．変数名は英文字で始まる英数字の文字列で，大文字小文字の区別をします．

>> b = 1.2e-3
>> A = [1 2*i;3  1 ]
>> w = -2+3*i

現在宣言されている変数名一覧はコマンドwho -variablesでリストアウトされます．

>> who -variables
*** local user variables:

A  b  w

定義した変数を削除するには，clearコマンドを使います．引数無しだと，全ての変数をクリアし，

>> clear b w

のように，特定の変数だけをクリアすることもできます．

コマンドライン上での入力と表示

基本的に数値やコマンド・関数をタイプしてenterを押せば，入力されて実行されます．標準状態では結果がすぐ表示されますが，表示させないようにすることもできます．一番簡単な方法は行の最後にセミコロン（; ）をつけることです．

>> a =[1 2 3]
a =
  1  2  3
>> a =[1 2 3];

四則演算，超越関数

四則演算は実数，複素数，行列，ベクトルおよびそれらを内容とする変数に対して実行可能です．ただし，行列，ベクトルの除算は以下のような意味となります．

例：\texttt{A/B}は$A\,B^{-1}$の意味

sinなどの関数も実数，複素数，行列，ベクトルおよびそれらを内容とする変数に対して実行可能ですが，行列，ベクトルに対しては各要素ごとの関数値となります．

>> A = [pi/2  pi;0   pi/6]
>> sin(A)
ans =
  1.00000  0.00000
  0.00000  0.50000

行列はベクトルの和や差はもちろん各要素ごとの和や差ですが，積や除は行列の積の定義に基づいた計算になります．しかし，場合によっては各要素ごとの積や除をとりたい場合があります．その時は，ピリオド（. ）を演算子の前につけます．

>> A = [1 2;3 4]
>> A*A
ans =
   7  10
  15  22
>> A.*A
ans =
   1   4
   9  16

行列・ベクトルの操作

要素，行，列，小行列の取り出し

行列およびベクトル内の特定の要素を取り出すには，変数名とカッコを用いて次のようにします．

>> a = [1 2 3;4 5 6;7 8 9]
>> a(2,3)
ans = 6

行列の特定の行や列を取り出すには要素の取り出しの際に列や行番号の代りにコロン(: )を使います．

>> a = [1 2 3;4 5 6;7 8 9]
>> a(2,:)
ans =
  4     5     6

小行列を取り出すにはコロンを用いて行もしくは列の範囲を指定します．

>> a = [1 2 3;4 5 6;7 8 9]
>> a(2:3,1:2)
ans =
  4     5
  7     8

行列の結合

複数の行列やベクトルをまとめて一つの行列にすることができます．ただし，対応する行数や列数は合わせておく必要があります．

>> a = [1;2;3;4];
>> b = [5;6;7;8];
>> c = [a b]
c =
  1  5
  2  6
  3  7
  4  8
 

転置行列

行列やベクトルの転置はシングルコーテーション（' ）を使います．

>> a = [1 2;3 4];
>> a'
ans =
  1     3
  2     4

特定の数列からなるベクトルの作成

Octaveの関数では，シミュレーションの時間軸やボード線図の周波数軸などをベクトルで与えることになります．そのときは，以下の関数を用いると便利です．

• 等差数列：linspace(始値,終値,点数)

たとえば，サンプリングタイム0.01で0から0.99までの点数100のベクトルを作るときは以下のようにします．

>> t = linspace(0,0.99,100);

• 等比数列：logspace(log(始値),log(終値),点数)

たとえば，0.01=10^{-2}から1000=10^3までを対数分割して100点のベクトルを作るとき．

>> omega = logspace(-2,3,100);

数値表示のフォーマット

演算結果や変数の内容を表示する桁数を指定できます．標準では5桁の数値で表示されます．これでは桁数がすくなすぎるとか，指数表示をしてもらいたい時はコマンドformatを使います．

>> 10/3
ans = 3.3333
>> format long
>> 10/3
ans = 3.33333333333333
>> format long e
>> 10/3
ans =  3.33333333333333e+00

デフォルトの5桁の表示に戻すには，引数無しでformatを実行します．なお，これらはあくまで表示上の変更で，内部計算は倍精度で行われています．

グラフ

Octaveではグラフ機能はgnuplotを呼び出して利用しているので，細かなコントロールにはgnuplotのコマンドを理解する必要があります．しかし，octave-forgeのパッケージにはMATLABとの互換性を高めるため，グラフ関連のコマンドが多く追加されていて（例えばgridなど），だいぶ使いやすくなっています．

plot関数の基本

1次元のベクトルの組を２次元グラフに表現できます．ベクトルXをx軸，Yをy軸としてグラフ化するには，plot(X, Y)とします．

>> X = linspace(0,10,100);
>> Y = sin(X);
>> plot(X, Y)

線の種類や色を変えたいときは，変数の後に続けてダブルコーテーションかシングルコーテーションでくくったフォーマット文字列を書けばいいです．

>> plot(X, Y, "b+")

この例は，色をblueとし，線ではなくて+記号でプロットします．その他の色や線類はhelpで調べて下さい．

なお，plot関数の引き数としてx軸, y軸, 線種を繰り返し用いることで複数のグラフを重ね書きすることもできます．

>> plot(X1, Y1, "r-", X2, Y2, "g-")

その他の装飾

plot関数でグラフを書いた後，xlabel, ylabel, titile関数で軸ラベルとタイトルをつけることができます．

>> xlabel("time (s)")
>> ylabel("sin(x)")
>> title("This is a sample graph")

また，グリッド（格子）をつけるには，

>> grid on

を実行すればいいです．

プロットの説明（凡例）が右上に自動的にでますが，これはlegend関数で変更できます．不要な場合は，

>> legend("off")

で消えます．

なお，環境によってはlegend関数がうまく機能しない場合があります．その場合は，plot関数でプロットする際に，色や線種を指定するのと同時に説明文をセミコロンで囲むことで指定できます．

>> plot(X, Y, "r-;sin;", X, Z, "g-;cos;")

また，グラフの表示範囲は自動的に設定されますが，複数のグラフを比較するときなど，一律に範囲を指定してプロットすべきです．その場合は，axis関数を使います．

>> axis([xmin xmax ymin ymax])

とすることで，x軸y軸それぞれの範囲を指定できます．なお，axis関数はグラフをプロットした後に実行しても効果がないので，後から実行した場合はもう一度プロットし直す必要があります．また，自動調整に戻す場合は，

>> axis("auto")

とします．

対数グラフ

plot関数のかわりにsemilogx関数を使うと片対数（x軸だけ対数）のグラフになります．両対数はloglog関数です．線種の指定などはplot関数と同じです

文書への張り込み

グラフをTeX文書などに取り込む場合は，EPSファイルとしていったん保存します．print関数を使うことでファイルに出力できます．

>> print("test.eps","-deps2")

一つ目の引数はファイル名で，二つ目はファイルの形式を指定するオプションです．-dに引き続いて，eps2はEPSのレベル2と言う意味です．epsc2とすればカラーになります．

スクリプトファイル

スクリプトファイル（mファイル）とは

これまでのようにコマンドを1行1行打ち込んでいっても作業はできますが，ある程度まとまった処理や繰り返し試行する処理などは，面倒ですし，間違いも見つかりにくいです．そこで，その一連の作業をmファイルというテキストファイルとして保存しておくと，実行や修正が非常に楽になります．

mファイルの作成と実行

mファイルはただのテキストファイルですから，メモ帳など使い慣れたエディタで作成・編集できます．保存する場所は自分のホームディレクトリ以下の適当なところとします．

ファイル名は変数名と同じく数字で始まらない英数字とし，拡張子は.mとします．実行はファイル名から.mを抜いた部分を入力すればいいです．カレント・ディレクトリにその名前のファイルがあれば実行されます．

path

自作のmファイルをあるディレクトリにまとめておき，カレント・ディレクトリがどこであっても，そのmファイルを実行できるようにすることができます．汎用性の高い関数やツールを作り，将来にわたり使いたい場合に便利です．

path関数によって，現在のコマンド検索pathの表示と変更ができます．

>> path

とすることで現在のpath設定が表示されます．（ゼミ開講現在，Windows環境ではなぜか表示されません．）これらに，自身のディレクトリを追加するには以下のようにします．

>> path(path, "/cygdrive/h/Documents/m_files/")

この例では，Documentsディレクトリの下のm_filesディレクトリにmファイルを入れておけば，カレントディレクトリがどこであっても実行できるようになります．なお，間違って

>> path("/cygdrive/h/Documents/m_files/")

とすると，pathがm_filesだけになり，自作のコマンド以外は全く実行できなくなりますので，注意です．（一度Octaveを終了して再度起動すれば元に戻ります）

制御系解析・設計コマンド類

システムの表現

Octave上では，システムは伝達関数表現（分子分母それぞれの係数を降べきの順に並べた2つのベクトルの組），状態空間表現（4つの行列の組），システム行列表現（1つの行列で表現）のいずれかで表現できます．しかし，後述する解析や設計の関数はシステム行列表現を引数としてとるので，その他の表現で与えられた物をシステム行列表現に変換する必要があります．変換を行う関数はtf()です（古いバージョンのoctaveの場合はtf2sys）．

以下の例は，&math(\frac{1}{s^2 + 2 s + 3});の伝達関数表現をシステム行列表現に変換する例です．

>> sys = tf([1], [1 2 3]);

周波数応答の計算

代表例としてボード線図を書いてみます．関数はbodeを利用します．

>> sys = tf([1], [1 2 3]);
>> bode(sys);

でゲイン線図と位相線図が描かれます．しかし，実際には複数のプロットを重ねたり，ゲイン線図だけが欲しかったりするので，bode関数は周波数応答の計算だけに利用して，グラフはあらためてplotまたはsemilog関数でプロットする方がいいです．

以下の例は，bodeの結果を変数gain, phaseに入れ，あらためてsemilog関数でプロットしています．なお，横軸である周波数軸を指定するために，logspace関数で変数omegaを作って利用しています．

>> sys = tf([1], [1 2 3]);
>> omega = logspace(-1,2,100);
>> [gain phase] = bode(sys, omega);
>> gain_db = 20*log10(gain);
>> axis([0.1 100 -120 0]);
>> semilogx(omega, gain_db);

データファイルの取り扱い

全変数の保存と復元

一旦Octaveを終了して次回続きをやりたいときなど，変数の内容をとっておきたいことがあります．スクリプトファイルとして作業手順を保存することはできますが，時間のかかる計算が必要な場合その度に計算し直すのも面倒です．

このような時，saveコマンドで変数の内容をファイルに保存できます．

>> save オプション ファイル名

とすると，カレントディレクトリに変数の内容を保持したファイルが作成されます．例えば，

>> a = [1 2 3];
>> B = [1 2 3; 4 5 6; 7 8 9];
>> save -ascii alldata.dat

とすれば，ファイルalldata.datに，現在のすべての変数の内容がアスキー形式（テキスト形式）で書き出されます．メモ帳など適当なテキストエディタで内容を確認できるはずです．なお，バイナリ形式で保存する場合は-binaryというオプションをつけます．テキスト形式にくらべて精度が高く，ファイルサイズも小さくなります．

保存したファイルの読み込みにはloadコマンドを使います．読み込むとすべての変数が復活します．試しに，いったんclearコマンドで変数を全て消してから，loadして見てください．ファイルがアスキー形式かバイナリ形式かは自動的に判別されます．

>> clear
>> load alldata.dat

なお，読み込もうとした変数がすでに存在する場合はエラーとなります．強制的に上書きする場合は，

>> load -force ファイル名

のように-forceオプションをつけます．

指定した変数だけの保存と復元

すべての変数を保存する必要がない場合や，一つのファイルに保存してある複数の変数のうちのいくつかだけを復元したい場合もあります．そのような場合，saveやloadコマンドでファイル名の後に変数名を羅列します．

>> save オプション ファイル名 変数名 変数名 ...

>> load オプション ファイル名 変数名 変数名 ...

できます．なお，先と同様にloadコマンドではファイルの形式は自動的に判別されるので，オプションは通常不要です．

タブまたはスペース区切りのテキストデータの読み込み

実験結果をOctaveで処理したりグラフ化する場合，実験用のプログラム（通常はC言語などで自作するでしょう）側はOctaveに渡したいデータをテキスト形式のファイルとして保存するのが簡単です．

saveコマンドに- asciiオプションをつけて保存したファイルは，一般に複数の変数について，変数名，タイプ，行サイズ，列サイズ，数列が列挙されています．このフォーマットを再現するようにプログラムを自作すれば，loadコマンドで先と同様に読み込めます．しかし，通常はタブまたはスペース区切りの数列だけのテキストファイルで保存しておき，Octave側で必要に応じて各変数に切り分けることが多いです．

例えば，次のような内容のテキストファイルresult.txtがカレントディレクトリにあるとします．

1 2 3
2 4 6
3 6 9
4 8 12
5 10 15

このファイルを

>> load result.txt

と読み込むことができます．この時，変数名はファイル名から拡張子をとったものとして読み込まれます．すなわち，この例の場合resultという名前の変数ができます．

先のsaveコマンドで保存したファイルの読み込みと異なるのは，このようにファイル名に基づいて変数名が決まることと，一度にひとつの変数しか読み込めないことです．ただし，後者は複数の変数の内容を結合して一つの行列としておき，それを一つのファイルでやり取りすることで対応できます．

例えば，先の例では，1列目が時間tのデータ，2列目が入力uのデータ，3列目は出力yのデータというように3つの変数を一つの行列にまとめたものであるという具合です．この場合，loadした後で，次のように列に切り分ければ，3つの変数として分離できます．

>> t = result(:,1);
>> u = result(:,2);
>> y = result(:,3);

他のアプリケーション向けのテキストデータ書き出し

Octaveでのシミュレーションや解析の結果を他のアプリケーション（例えばExcelなどのスプレッドシートやIgorやNGraphといったグラフ化ソフトなど）に渡したい時はテキストファイルとして書き出します．ただし，save -ascii 変数名 変数名 ...で保存できるOctave固有のフォーマットのファイルをそのまま読み込めるとは限らないので，以下のような工夫が必要です．

まず，書き出したい変数を結合して一つの変数にしておきます．このときサイズ（行数）が揃っていないとエラーになりますので，揃えるか，一緒にすることをあきらめて別のファイルにするか決めます．シミュレーション結果などは，時刻ベクトル，入力ベクトル，出力ベクトルなどが大抵同じ行数でしょうから，まとめて1つの行列にできるでしょう．

例えば，時刻t，入力u，出力yの3つのベクトルを変数outfileにまとめてからテキストファイルoutfile.txtに保存するには以下のようにします．ただし，3つのベクトルは列ベクトルであるとします．

>> outfile = [t u y];
>> save -ascii outfile.txt outfile

保存されたファイルの最初の5行は，日付や変数名，サイズなどの情報で先頭文字が#になっています．受け取る側のアプリケーションでその行を無視するように設定できれば，そのまま使用できますが，そうでない場合はテキストエディタで削除してから読み込ませる必要あります．

最初から先頭の不要な行をつけないで保存するには，saveコマンドを利用せず，C言語と同様にfopen()によるファイルのオープン，fprintf()による書式付き書き出し，fclose()によるファイルのクローズの手順を踏みます．

>> outfile = [t u y];
>> fid = fopen("outfile.txt", "w");
>> fprintf(fid, "%g %g %g\n", outfile');
>> fclose(fid);

このようにfprintf()の引数として行列やベクトルがとれて，いっぺんに複数行書きだしてくれるところがC言語より便利な点です（C言語ならfor()ループが必要でしょう）．ただし，fprintf()の引数の行列にアポストロフィー（' ）を付ける必要があります．すなわち，行列を転置して引数とする必要があります．また，書式指定内には%gが行列の列数と同じ数だけ必要です．

なお，Octave-forgeパッケージに限りますが，dlmwrite()という関数を用いると簡単に数列だけのファイルとして保存できます．

>> dlmwrite("outfile.txt", outfile, " ");

とすればいいです．3つ目の引数は区切りにする文字の指定で，ここではスペースにしています（省略するとカンマになります）．

Octaveスクリプトのいろいろ

キーボード入力，画面への出力

スクリプト実行中にキーボードから数値や文字列を入力させるにはinput命令を使います．同時にメッセージとして表示させる文字列を指定することができます．

画面への出力はC言語と同じくprintf命令を使います．

a = input("Please input a: ");
printf("a = %g\n",a);

条件判定

条件判定にはC言語と同じくif命令が利用できます．ただし，C言語のように括弧で条件式は囲みません．また，中括弧（{ }）は使用せず，endifでif命令の範囲を指定します．さらに，C言語ではelseとifを分けて書きますが，Octaveではelseifと続けて書きます．

a = input("Please input a: ");
if a > 0
	printf("a > 0\n");
elseif a > -5
	printf("-5 < a <= 0\n");
else
	printf("a < = -5\n");
endif 

繰り返し

while命令

while命令は条件が成り立っている間，命令群を繰り返します．やはり中括弧は用いずに，命令群の終わりはendwhileで指定します．なお，中断は\break命令です．

while 条件
　命令群
endwhile

sum = 0;
while  sum < 100
　a = input("Please input a: ");
　sum = sum + a;
　printf("sum = %g\n", sum);
endwhile

for命令

for命令はC言語のforを限定的にした感じで，インデックスとする変数の初期値，増分，終値を指定して命令群を繰り返します．命令群の終わりはendfor，中断はbreakなのはwhileと同様です．なお，増分を省略すると1ととられます．

for インデックス変数 = 初期値:増分:終値
　命令群
endfor

sum = 0;
for  i = 1:2:100
　sum = sum + i;
　if sum >100
　　break
　endif
endfor
printf("On i = %d, sum = %g\n", i, sum);

初期値:増分:終値という書き方は，ベクトルの作成にも使えますので覚えておくと便利です．

ユーザ定義関数

スクリプト内でユーザ独自の関数を定義したり，別ファイルで定義したりできます．

スクリプト内で定義する場合は以下のようにします．

function [戻値1, 戻値2, ...] = 関数名(引数1, 引数2, ...)
% 関数の説明
　命令群
endfunction

関数の終わりを示すためにendfunction命令が必要です．また，function命令直後のコメント文はhelpで表示されるので，簡単な説明を入れておくと良いでしょう．

別ファイルで関数を定義する場合は，ファイル名を関数名に拡張子.mをつけます．この場合endfunctionはあっても無くても構いません（MATLABでは書かないので書かない方がいいかも）．さらに，関数名は既存の関数とかぶらないように，オリジナルなものになるよう注意が必要です．

なお，C言語の関数と同様に，関数内で変数はローカル変数となるので，変数名がかぶっても内容は別物となることに注意が必要です．

以下の例は，二つのベクトルの各要素ごとに差をとり，それを各々2乗したものの和を計算する関数の例です．（標準関数でsumsqというものがあるので作るまでもないのですが...）

function err = sqerror(a, b)
% usage: err = sqerror(a, b)
%  a: column vector
%  b: column vector, same siaze as a
%  err: sum of square errors between a and b
c  = a - b;
err = c' * c;
endfunction