本報告では, 総合的な科学技術計算環境の1つであるMATLABについて解説する.
MATLAB は, 行列計算を中心とした数値計算, 数式処理, ビジュアリゼーション, GUI構築などの機能を総合した科学技術計算環境であり, またプログラミング言語でもある. この種のプログラムとしてはMathematicaが有名であるが, Mathematicaが数学や物理学の研究者に特に人気が高いのに対し, MATLABは工学研究者, 特にシステム制御や信号処理の研究者に利用者が多い. システム制御や信号処理の分野では, シミュレーションレベルまでの研究ならFORTRANやCを使用せずに MATLABでプログラミングを行うことが常識化しつつあり, 教育に関しても座学にMATLABを組み合わせて 直観的に分かりやすい講義を行うことが主流になりつつある.
MATLAB の一般的な利用法や, システム制御における MATLAB の重要性については既に多くの報告が存在するので, 以下では信号処理におけるMATLABの利用について解説して行きたい. なお, 山梨大学総合情報処理センターに MATLAB が導入されたのは 前回の機種更新時 (1998年) であるが, センター研究報告にMATLABに関する報告が掲載されるのは 今回が初めてのようなので, 以下ではまずMATLABの一般的な紹介から始める.
MATLAB の最大の特徴は, 行列計算をほぼ数式そのままの形でプログラム化できることであり, これによってアルゴリズムのコーディングの能率を非常に高くすることが可能である.
また, MATLAB は基本的にはインタープリターなので, 新たに作成したプログラムをそのままコマンドとして システムに付け加えて行くことが可能である. この機能を利用して多数のアプリケーションパッケージ (MATLAB では TOOLBOX と呼ぶ) が開発されており, 特に制御理論関係の TOOLBOX は充実している.
>> A = [1 2 3 > 4 5 6 > 7 8 0] A = 1 2 3 4 5 6 7 8 0 >> A = [1,2,3;4,5,6;7,8,0] A = 1 2 3 4 5 6 7 8 0 >> b = [1 ; 2; 3] b = 1 2 3ここで, > 記号が現れている行はキーボードからの入力を表し, > 記号が現れていない行はMATLABが返してきた応答を表示している.
>> 2 * A ans = 2 4 6 8 10 12 14 16 0 >> A + A ans = 2 4 6 8 10 12 14 16 0 >> A - A ans = 0 0 0 0 0 0 0 0 0 >> A * b ans = 14 32 23一方, 数の除算に相当するものは行列の場合は連立方程式の求解であるが, これも, 以下のように計算することができる.
>> A ¥ b ans = -3.3333e-01 6.6667e-01 1.3878e-17 >> b' / A' ans = -3.3333e-01 6.6667e-01 1.3878e-17バックスラッシュ (日本語端末では ¥ 記号) は A*x=b を解く操作を表しており, スラッシュ (/) は x*A'=b' を解く操作を表している. また, シングルクォート (') は転置を示す.
行列としての乗除算ではなく, 対応する要素毎にスカラーとしての乗除算を行いたい場合には,
>> [1 2 3] .* [2 3 4] ans = 2 6 12 >> [1 2 3] ./ [2 3 4] ans = 0.50000 0.6667 0.75000のように .*,./ という記号を使う. べき乗についても, 行列としてのべき乗を表す ** と要素毎のべき乗を表す .** がある.
>> A ** 2 ans = 30 36 15 66 81 42 39 54 69 >> A .** 2 ans = 1 4 9 16 25 36 39 54 69
>> A = [1 2 3 4 5 6 7 8 0] A = 1 2 3 4 5 6 7 8 0] >> det(A) ans = 27.000 >> rank(A) ans = 3 >> trace(A) ans = 6 >> [U S V] = svd(A) U = -0.230357 -0.396071 -0.888855 -0.607284 -0.655208 0.449343 -0.760356 0.643297 -0.089596 S = 13.20146 0.00000 0.00000 0.00000 5.43876 0.00000 0.00000 0.00000 0.37605 V = -0.604629 0.273256 0.748167 -0.725676 0.198242 -0.658858 -0.328356 -0.941292 0.078432特異値分解 svd の例に見られるように, MATLAB の関数は複数の返り値を返すことができ, 返り値はスカラーだけでなくベクトルや行列であっても良い.
>> x=[0, 0.5*pi, pi, 1.5*pi, 2*pi] x = 0 1.5708 3.1416 4.7124 6.2832 >> cos(x) ans = 1.0000 0.0000 -1.0000 0.0000 1.0000この機能をうまく利用すると, for や while のようなループ制御命令をほとんど使用することなくプログラムを書くことができる.
>> A = [1, 2, 3 ; 4, 5, 6 ; 7, 8, 9]; >> A A = 1 2 3 4 5 6 7 8 9同じことを display という命令を用いて明示的に行うこともできる.
>> display(A); A = 1 2 3 4 5 6 7 8 9それ以外に, fopen,fclose, fread,fwrite, printf,fprintf,sprintf, scanf,fscanf,sscanf などの入出力関数があり, これらの使用法はほぼC言語と同じである.
if k==3 k = k + 1; else k = k - 1; end for m=1:2:10 n = m + 1; end n = 1; while n < 10 n = n + 1; end等々である.
なお, ここで 1:2:10 などとあるのは制御構造の一部ではない. 1:2:10自体が, 1から2刻みで10までの値をとるベクトルを表現しており, m=1:2:10はm=[1,3,5,7,9]と書いても同じ結果となる. したがって,
for m=[-3, 2, 10, 0.5] n = m + 1; endのようにインデクスをでたらめな順序にすることもできる. また,
>> x = [0:0.2:1] x = 0 0.2 0.4 0.6 0.8 1.0という記法によって等差数列を要素とするベクトルを簡単に入力することが可能である.
>> x = randn(1,100); >> plot(A)と入力すると下の図のように, xのグラフが描かれる.
plot の代わりに mesh を使うと, 3次元グラフが描かれる.
>> A=cos(2*pi*[0:0.05:1])' * cos(2*pi*[9:0.05:1]); >> mesh(A)
A = [1, 2, 3 ; 4, 5, 6 ; 7, 8, 9] b = [2 ; 4 ; 6] x = A ¥ bここで, MATLAB を起動して,
>> test1と打ち込んでみる. 結果は, MATLAB のコマンドラインから直接上の内容を打ち込んだ場合と全く同じになる. このように, script M-file には MATLAB のコマンドを書き並べて置くことができ, MATLAB は scriptM-file から拡張子 .m を除いた部分を新しい命令として認識する.
ところで, 上のような書き方をすると test1 コマンドの途中結果はすべて画面に表示される. 場合によっては, 途中結果の表示は必要ない場合もある. その場合には,
>> c = [2 7 3];のように, 各コマンドの終りにセミコロン ; を付ける. セミコロンを付けると実行結果を表示しなくなる.
function y = entropy(p) % y = entropy(p) は離散確率分布 p のエントロピーを計算する. y = -sum(p .* log2(p));function M-file は, 必ず function 文で始まらなければならない. ここで, function は関数を表す決まり文句, entropy が関数の名前, p が引数, y が戻り値である. 次に, 行の中に%記号があると % 以降の部分はコメントとみなされ無視される. また, function M-file の名前は関数の名前に拡張子 .m を付けたものでなければならない. 上の例では, function M-file の名前は entropy.m でなければならないことになる. 上記の部分に続けて, 通常の MATLAB の命令を書いて行くが, script M-file の場合と違って, 途中結果は保存されず, 保存されるのは返り値 y のみである.
適当なエディタで上のような内容のファイルを作成し, entropy.m という名前で保存する. その後, MATLAB を立ち上げ, 適当な離散確率分布を入力してエントロピーを計算してみる. 例えば, 統計的にある日の天気が「晴れ」30%,「曇」30%,「雨」40%ならば, その日の天気に関する不確実さは,
>> p = [0.3 0.3 0.4]; >> y = entropy(p) y = 1.5710となる.
前節のエントロピー関数を例にとって, MEX-file の作成法を示す.
>> mcc -x entropy >> ls ans = entropy.c entropy.h entropy.m entropy.mexsg entropy_mex.cmcc が,function M-file を MEX-file に変換する命令である. lsコマンドによってカレントディレクトリのファイル一覧を表示すると, entropy.mをC言語に変換したentropy.cおよびentropy_mex.c, ヘッダファイルentropy.h, およびentropy.mexsgというMEX-fileが作成されていることが分かる (mexsgの最後のsgは,計算機のアーキテクチャを区別するための記号である). この状態で entropy を呼び出すと, entropy.m ではなく entropy.mexsg の方が実行される.
また, TOOLBOX ではないが MATLAB と協調して動作するプログラムとして以下のものがある.
M-1 X(ω) = Σ x(k) exp(-jωk) k=0という積和演算が必要であるし, 信号u(k)をフィルタリングして出力y(k)を得ようとすれば,
N N y(k) = Σ b(l) x(k-l) - Σ a(l) y(k-l) l=0 l=1という積和演算が必要となる (a(l),b(l) はフィルタ係数). このように, 信号処理においては積和演算の形で書くことのできるアルゴリズムがきわめて多い.
ここで, 信号やフィルタ係数をインデックスに関して 適切な数だけまとめてベクトルで表記する ( T は転置を表す).
x = [x(0) x(1) ... x(M-1)]T e = [exp(-jω0) exp(-jω1) ... exp(-jω(M-1))]T u = [u(k) u(k-1) ... u(k-N)]T v = [y(k-1) y(k-2) ... y(k-N)]T a = [a(1) a(2) ... a(N)]T b = [b(0) b(1) ... b(N)]Tこれを用いると積和演算はベクトルの内積で表現することができる.
X(ω) = eT x y(k) = bT u - aT vさらに,
X = [X(ω0) X(ω1) ... X(ωM-1)]T y = [y(0) y(1) ... y(M-1)]T E = [exp(-jω00) exp(-jω01) ... exp(-jω0 (M-1)) exp(-jω10) exp(-jω11) ... exp(-jω1(M-1)) : exp(-jωM-10} exp(-jωM-11) ... exp(-jωM-1(M-1))] A = [a(N) ... a(1) 0 ... 0 0 a(N) ... a(1) 0 ... 0 : 0 ... 0 a(N) ... a(1)] B = [b(N) ... b(0) 0 ... 0 0 b(N) ... b(1) 0 ... 0 : 0 ... 0 b(N) ... b(1)]のように定義すると, 積和演算は行列とベクトルの積で表現されることになる.
X = E x y = B u - A vこの式を MATLAB で表現すれば, 以下のようになる.
>> X = E * x >> y = B * u - A * vこのように, 信号処理アルゴリズムをベクトルと行列で書き表すと, アルゴリズムの表記そのものが ほとんどそのまま MATLAB のプログラムになることが分かる.
そこで, フィルタ設計に関する CAD ツールが重要となる. MATLAB の Signal Processing TOOLBOX には, インタラクティブにフィルタを設計することができる Filter Desinger が付属している. フィルタの種類 (lowpass,highpass,bandpass,Hilbert transform,etc), 遮断周波数, リプル特性 (等リプル,最小二乗,バタワース,etc) などを指定すると, Filter Designer は適切なフィルタ次数を選択し, フィルタを設計, その特性を表示してくれる.
また, 設計したフィルタをより詳しく解析したい場合には, Filter Viewer と呼ばれるツールを使用することができる.
MATLAB が開発される際の基礎になった Linpack や Lapack といったソフトウェアパッケージは元々フリーなソフトウェアであるので, これらに基づいてフリーなMATLABクローンを作成しようとするプロジェクトがいくつも存在している. ここでは, その中でOctaveとScilabを紹介する (これ以外にも, Algae,Rlab,MaTX など MATLAB-like なプログラム言語は多数存在する. これらのプログラムは, Webで検索をかければ容易に見つかるはずなので, ここではこれ以上触れない). これらのクローンには MATLAB ほど多様な TOOLBOX は付属していないが, ソースコードが公開されているので MATLAB でサポートされない機種/OSの上でも利用できる可能性がある. また, MATLABが個人的に購入するにはいささか高価すぎる商品であるのに対して, これらのクローンは個人でも気軽に使用でき, 特に学生の自習用には適している.
Octave は, 基本的な言語仕様の点で MATLAB 互換となるように設計されており, Octave で作成したプログラムの多くは修正なしでそのまま MATLAB で実行可能である (MATLABの方が機能が豊富なので,逆は真でない). MATLABクローンの中でも, MATLAB互換という点では最右翼に位置するプログラムである. ただし, 入出力の GUI 化にはあまり力を入れておらず, Simulinkに相当するようなプログラムも付属していない.
本報告では, MATLAB の利用法, 信号処理との関係, MATLAB 互換プログラムなどについて紹介した. いささか散漫な内容となってしまったが, ご寛恕いただければ幸いである.
本報告を読まれて MATLAB に興味をお持ちになった方は, ぜひ実際に MATLAB を使用してみることをお勧めしたい. 始めは少々戸惑うかもしれないが, 一度 MATLAB に慣れるともう MATLAB を手放せなくなるであろうと思う.