[研究報告の目次へ戻る]

MATLABと信号処理

山梨大学工学部電気電子システム工学科
大木 真
ohki@es.yamanashi.ac.jp
http://www.sp.es.yamanashi.ac.jp/‾ohki/

  1. まえがき
  2. MATLAB の紹介
  3. MATLAB と信号処理
  4. MATLAB類似のフリーソフトウェア
  5. まとめ
  6. 参考文献

1. まえがき

本報告では, 総合的な科学技術計算環境の1つであるMATLABについて解説する.

MATLAB は, 行列計算を中心とした数値計算, 数式処理, ビジュアリゼーション, GUI構築などの機能を総合した科学技術計算環境であり, またプログラミング言語でもある. この種のプログラムとしてはMathematicaが有名であるが, Mathematicaが数学や物理学の研究者に特に人気が高いのに対し, MATLABは工学研究者, 特にシステム制御や信号処理の研究者に利用者が多い. システム制御や信号処理の分野では, シミュレーションレベルまでの研究ならFORTRANやCを使用せずに MATLABでプログラミングを行うことが常識化しつつあり, 教育に関しても座学にMATLABを組み合わせて 直観的に分かりやすい講義を行うことが主流になりつつある.

MATLAB の一般的な利用法や, システム制御における MATLAB の重要性については既に多くの報告が存在するので, 以下では信号処理におけるMATLABの利用について解説して行きたい. なお, 山梨大学総合情報処理センターに MATLAB が導入されたのは 前回の機種更新時 (1998年) であるが, センター研究報告にMATLABに関する報告が掲載されるのは 今回が初めてのようなので, 以下ではまずMATLABの一般的な紹介から始める.


2. MATLAB の紹介

2.1 概要

MATLAB は米 MathWorks 社の製品であり, 日本ではサイバーネット社が販売している. MATLABは, LinpackやLapackといったFortran言語の数値計算サブルーチンパッケージを 対話的に使用することから出発したプログラムである. 現在では, 数式処理やビジュアリゼーション, 対話型(GUI)プログラム作成を含む総合的な科学技術計算環境となっている.

MATLAB の最大の特徴は, 行列計算をほぼ数式そのままの形でプログラム化できることであり, これによってアルゴリズムのコーディングの能率を非常に高くすることが可能である.

また, MATLAB は基本的にはインタープリターなので, 新たに作成したプログラムをそのままコマンドとして システムに付け加えて行くことが可能である. この機能を利用して多数のアプリケーションパッケージ (MATLAB では TOOLBOX と呼ぶ) が開発されており, 特に制御理論関係の TOOLBOX は充実している.

2.2 インタープリターの使用法

MATLAB を起動すると, コンソールから直接命令を入力するインタープリターとして使用することができる. まず, このインタープリターの使い方について紹介する.


図-1 MATLAB 起動画面

2.2.1 ベクトルと行列の入力

MATLAB では, 手書きでベクトルや行列を表記する場合と同じように, 列をスペース (空白) で分離し行を改行で分離するという書き方をすることができる. また, 列をコンマ (,) で行をセミコロン (;) で分離するという書き方も可能である.
  >> 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.2.2 基本演算

行列とスカラー, 行列とベクトル, および行列と行列の加減算や乗算は, 普通の数のように +,-,* で計算できる.
  >> 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

2.2.3 行列操作関数

行列式, 階数 (ランク), トレース, 固有値分解, 特異値分解などの基本的な行列操作はそれぞれ関数が用意されている.
  >> 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 の関数は複数の返り値を返すことができ, 返り値はスカラーだけでなくベクトルや行列であっても良い.

2.2.4 数値計算関数

MATLAB では, sin,cos などのような一般的な数値計算関数も引数としてベクトルや行列をとることができる. その場合, 自動的に返り値も適切なサイズの行列となる.
  >> 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 のようなループ制御命令をほとんど使用することなくプログラムを書くことができる.

2.2.5 入出力関数

MATLAB では, コマンドラインから単に変数名を入力すると, その内容が画面に表示される.
  >> 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言語と同じである.

2.2.6 制御構造

MATLAB では, 他のプログラミング言語と同様に様々な制御構造を利用することができる.
  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
という記法によって等差数列を要素とするベクトルを簡単に入力することが可能である.

2.2.7 ビジュアリゼーション

MATLAB には処理結果をグラフ化する機能があり, 作図されたグラフィックを取り扱うGUIも用意されている. 例えば,
  >> x = randn(1,100);
  >> plot(A)
と入力すると下の図のように, xのグラフが描かれる.


図-2 2次元グラフ

plot の代わりに mesh を使うと, 3次元グラフが描かれる.

  >> A=cos(2*pi*[0:0.05:1])' * cos(2*pi*[9:0.05:1]);
  >> mesh(A)


図-3 3次元グラフ

2.3 プログラム言語としての使用法

MATLAB は, カレントディレクトリおよび予め指定されたディレクトリに .m あるいは.mexという拡張子のついたファイルが存在すると, それらをMATLABのプログラムファイルとみなす. プログラムファイルには, 以下の3種類のものがある.
  1. script M-file
  2. function M-file
  3. MEX-file

2.3.1 script M-file

script M-file は, MATLAB のコマンドラインから入力する命令と同じものを書き込んだテキストファイルである. 例えば, カレントディレクトリにtest1.mというファイルを作成して, その中に以下のような内容を書き込んだとする.
  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];
のように, 各コマンドの終りにセミコロン ; を付ける. セミコロンを付けると実行結果を表示しなくなる.

2.3.2 function M-file

function M-file は, 他の言語でサブルーチンや関数と呼ばれる機能を実現するためのものである. 例を取って, functionM-fileの文法を説明する. 以下は, 離散確率分布を表すベクトルpを与えて, 対応するエントロピーを返す関数の例である.
  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
となる.

2.3.3 MEX-file

MATLAB には, function M-file を C 言語に変換し, これをコンパイルして MATLAB から直接呼出し可能な実行オブジェクトに変換する機能がある. これらの実行オブジェクトのことをMEX-fileと呼ぶ. プログラムの内容にもよるが, MEX-fileは元のfunctionM-fileよりも高速になるのが普通である.

前節のエントロピー関数を例にとって, MEX-file の作成法を示す.

  >> mcc -x entropy
  >> ls
  ans =
  entropy.c
  entropy.h
  entropy.m
  entropy.mexsg
  entropy_mex.c
mcc が,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 の方が実行される.

2.4 TOOLBOX

MATLAB には TOOLBOX と呼ばれる多様なアプリケーションが用意されている. TOOLBOXと呼ばれるものの実態は, functionM-fileとMEX-fileによって構成された関数群である. 山梨大学総合情報処理センターの教育用計算機には, 以下の TOOLBOX がインストールされている (2002年3月現在. 2002年4月に総合情報処理センターの機種が更新されるので, 以降は若干の変更が有り得る). TOOLBOX は次々と新しいものが開発されており, 特にシステム制御関係の TOOLBOX は最新の理論を逸早く取り入れることで定評がある.

また, TOOLBOX ではないが MATLAB と協調して動作するプログラムとして以下のものがある.

Simulink については, 3節で述べる.

3. MATLABと信号処理

3.1 積和演算

信号処理においては, 積和演算 (連続時間信号では積分であるが, 計算機上で扱う場合には連続時間信号であっても離散化して扱わざるをえないので, 以下では離散の場合のみについて述べる) がきわめて重要な役割を果たす. たとえば, 信号 x(k) のフーリエ変換を求めようとすれば,
           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 のプログラムになることが分かる.

3.2 フィルタ設計

目的とする仕様を与えて必要とするフィルタを設計することは, 信号処理システムの設計において中心的な作業の1つである. しかし, 実はフィルタ設計を詳しく述べた参考書はきわめて少ないのが現状である. 一般的な信号処理の教科書には, フィルタ設計に関することはごく入門的な内容しか書かれていないのが普通であるし, 原論文を読んでフィルタ設計法を理解することは フィルタ設計を専門としない人間の手に余る.

そこで, フィルタ設計に関する CAD ツールが重要となる. MATLAB の Signal Processing TOOLBOX には, インタラクティブにフィルタを設計することができる Filter Desinger が付属している. フィルタの種類 (lowpass,highpass,bandpass,Hilbert transform,etc), 遮断周波数, リプル特性 (等リプル,最小二乗,バタワース,etc) などを指定すると, Filter Designer は適切なフィルタ次数を選択し, フィルタを設計, その特性を表示してくれる.


図-4 Filter Designer

また, 設計したフィルタをより詳しく解析したい場合には, Filter Viewer と呼ばれるツールを使用することができる.


図-5 Filter Viewer

3.3 信号とシステムの解析

信号のスペクトル解析や考えている信号処理システムの特性解析を行い, その結果をグラフで示すことは, 実際の信号処理において中心的な作業の1つである. これらに関しては, 一般的な信号処理の教科書にもかなり詳しく述べられている. しかし, 実際にその処理をプログラムしようとすると結構面倒なことが多い. MATLAB の Signal Processing TOOLBOX には, これらの処理をインタラクティブに処理することができる SPTool が付属している. 前節でのべた Filter Desinger も, この SPTool から呼び出すことができる. SPTool には, 波形表示 (Signal Browser) やスペクトル解析 (Spectrum Viewer) などが含まれている.


図-6 SPTool


図-7 Signal Browser


図-8 Spectrum Viewer

3.4 Simulink

Simulink はブロックダイアグラムシミュレータである. 信号処理・システム制御に関する様々な機能ブロック (積分器,遅延器,加算器,乗算器,非線形関数,信号発生器,etc) をビジュアルに組み合わせてブロックダイアグラムを作成し, そのブロックダイアグラムのシミュレーションを行うことができる. また, MATLAB側で定義したシステムオブジェクト(状態空間モデルなど) を Simulink のブロックダイアグラムの中に取り込むことができる.


図-9 Simulink


4. MATLAB類似のフリーソフトウェア

MATLAB が開発される際の基礎になった Linpack や Lapack といったソフトウェアパッケージは元々フリーなソフトウェアであるので, これらに基づいてフリーなMATLABクローンを作成しようとするプロジェクトがいくつも存在している. ここでは, その中でOctaveとScilabを紹介する (これ以外にも, Algae,Rlab,MaTX など MATLAB-like なプログラム言語は多数存在する. これらのプログラムは, Webで検索をかければ容易に見つかるはずなので, ここではこれ以上触れない). これらのクローンには MATLAB ほど多様な TOOLBOX は付属していないが, ソースコードが公開されているので MATLAB でサポートされない機種/OSの上でも利用できる可能性がある. また, MATLABが個人的に購入するにはいささか高価すぎる商品であるのに対して, これらのクローンは個人でも気軽に使用でき, 特に学生の自習用には適している.

4.1 Octave

Octaveは, Wisconsin 大学化学工学科の先生たちによって講義用の補助プログラムとして作られた. 現在では, Free Software Foundation の GNU プロジェクトの一部として GPL (GNU General Public License) にしたがって配布されている.

Octave は, 基本的な言語仕様の点で MATLAB 互換となるように設計されており, Octave で作成したプログラムの多くは修正なしでそのまま MATLAB で実行可能である (MATLABの方が機能が豊富なので,逆は真でない). MATLABクローンの中でも, MATLAB互換という点では最右翼に位置するプログラムである. ただし, 入出力の GUI 化にはあまり力を入れておらず, Simulinkに相当するようなプログラムも付属していない.


図-10 Octave 起動画面

4.2 SciLab

Scilabは, フランスの研究機関 INRIA で開発されている. Scilab の文法は MATLAB と良く似てはいるが, 細部には違いがある. また, scriptfileやfunctionfileの作成の仕方もMATLABとは若干異なっている. 一方, ビジュアリゼーションやシミュレーション環境には力を入れており, 独自の GUI や Simulink 類似のブロックダイアグラムシミュレータ(Scicos)を提供している.


図-11 独自の GUI を持つ Scilab 起動画面


図-12 ブロックダイアグラムシミュレータScicos


5. まとめ

本報告では, MATLAB の利用法, 信号処理との関係, MATLAB 互換プログラムなどについて紹介した. いささか散漫な内容となってしまったが, ご寛恕いただければ幸いである.

本報告を読まれて MATLAB に興味をお持ちになった方は, ぜひ実際に MATLAB を使用してみることをお勧めしたい. 始めは少々戸惑うかもしれないが, 一度 MATLAB に慣れるともう MATLAB を手放せなくなるであろうと思う.


6. 参考文献

6.1 一般的な参考文献

  1. G. J. Borse: MATLAB 数値解析,オーム社.
  2. D. M. Etter,D. C. Kuncicky: MATLAB ビギナーズガイド,山海堂.
  3. 芦野 隆一,レミ・ヴァイアンク−ル:はやわかりMATLAB,共立出版.
  4. 芦野 隆一・他: MATLAB による微分方程式とラプラス変換,共立出版.
  5. 大石 進一: MATLAB による数値計算,培風館.
  6. 高谷 邦夫: Matlab の総合応用,森北出版.
  7. 上坂 吉則: MATLAB プログラミング入門,牧野書店.
  8. 小国 力: MATLAB と利用の実際,サイエンス社.
  9. 小国 力: MATLAB グラフィクス集,朝倉書店.
  10. 小国 力: MATLAB による線形計算ソフトウェア,丸善.
  11. 小国 力: MATLAB によるカオス/フラクタル,朝倉書店.
  12. 小林 一行: MATLAB ハンドブック,秀和システム.
  13. 小林 一行: MATLAB 活用ブック,秀和システム.

6.2 システム制御関連の参考文献

  1. 足立 修一: MATLABによる制御のためのシステム同定,東京電機大学出版局.
  2. 足立 修一: MATLABによる制御工学,東京電機大学出版局.
  3. 川田 昌克,西岡 勝博: MATLAB/Simulinkによるわかりやすい制御工学,森北出版.
  4. 野波 健蔵・編著:MATLABによる制御理論の基礎,東京電機大学出版局.
  5. 野波 健蔵・編著:MATLABによる制御系設計,東京電機大学出版局.
  6. 尾形 克彦: MATLABによる制御工学問題の解法,新技術開発センター.

6.3 信号処理関連の参考文献

  1. J. H. McClellan, R. W. Schafer and M. A. Yoder: DSP First, Prentice-Hall (マクレラン,シェーファー,ヨーダー (荒實・訳): MATLABによるDSP入門,ピアソン・エデュケーション).
  2. 青木 由直: オペレ−タ法ディジタル信号処理   C言語とMATLABによる処理の実際,コロナ社.
  3. 足立 修一: MATLABによるディジタル信号とシステム,東京電機大学出版局.
  4. 金城 繁徳,尾知 博: 例題で学ぶディジタル信号処理,コロナ社.
  5. 高井 信勝: 「信号処理」「画像処理」のためのMATLAB入門,工学社.
  6. 真田 幸俊: MATLAB/Simulink による CDMA,東京電機大学出版局.
  7. 樋口 龍雄,川又 政征: MATLAB対応 ディジタル信号処理,昭晃堂.

6.4 Octave/Scilab 関連の参考文献

  1. 大石 進一: Linux数値計算ツール,コロナ社.
  2. 浪花 智英: Octave/Matlabで見るシステム制御,科学技術出版.

6.5 ソフトウェア開発元

  1. MathWorks: http://www.mathworks.com/
  2. サイバーネットシステム: http://www.cybernet.co.jp/
  3. Octave開発元(ウィスコンシン大学): http://www.che.wisc.edu/octave/
  4. Scilab開発元(INRIA): http://www-rocq.inria.fr/scilab/


 [ 研究報告の目次へ戻る ]