MATLAB Lecture
				   北原鉄朗
				  2002.04.22


1. MATLAB 起動


コマンドラインから
    matlab
（統合環境が立ち上がる）

コマンドラインから
    matlab -nodesktop
とやると，kterm 内にインタプリタのみ立ち上がる（だから軽い）．

※ path を通す必要がある（当研究室の場合は /home/lab5/speech/matlab/bin/）．



2. WAV ファイルの読み込み／再生／可視化


WAV ファイルを読み込む．
>> [y, fs, bits] = wavread('test.wav');

※ matlab は，基本データ構造が行列で，関数の出力引数を複数設定できる．
上の文は，test.wav という名の WAV ファイルを読み込んで，y にサンプル列，
fs にサンプリング周波数，bits に量子化ビット数が格納される．
サンプル列は -1〜+1 の実数値に変換される．

波形をグラフ表示する．
>> plot(y);

波形を再生する．
>> sound(y, fs);

波形の一部を取り出す．
>> y1 = y(1 * fs + (0 : 4095));

※ 最初から1秒後の時点から4096点分が抽出される．

取り出した部分をグラフ表示する．
>> plot(y1);

取り出した部分を FFT する．
>> Y1 = fft(y1);

グラフ表示する．
>> plot(abs(Y1));

※ グラフから分かるように半分のところで線対称になっているため，
   4096点のうち意味があるのは2048点である．

半分だけグラフ表示する．横軸のスケールを Hz にする．
>> plot((1 : 2048) / 4096 * fs, abs(Y1(1 : 2048)));

スペクトログラムを表示する．
>> specgram(y, 4096, fs, hanning(4096), 4096 - 0.01 * fs);

※ 4096 は窓長，hanning は窓関数の種類，4096 - 0.01 * fs は
   オーバーラップさせる点の長さを示す．



3. カーブフィッティング


例題：次のデータのカーブフィッティングをしてみよう．
(x, y) = (1, 5.5), (2, 43.1), (3, 128), (4, 290.7), (5, 498.4)

まず，データを入力する．
>> x = [1  2  3  4  5]';
>> y = [5.5  43.1  128  290.7  498.4]';

データをグラフに表示してみる．
>> plot(x, y, 'o');

※ plot 関数のオプションとして 'o' を指定すると，折れ線グラフではなく
   点で表示することができる．

線形で近似する．
>> a = [ones(length(x), 1)  x] \ y;

>> a
a =
 -176.8800
  123.3400

※ これは，y 切片が -176.8800，傾きが 123.3400	であることを示す．

>> x2 = (1 : 0.1 : 5)';
>> y2 = [ones(length(x2), 1)  x2] * a;

近似直線を表示する．
>> plot(x, y, 'o', x2, y2);

次に，2次曲線で近似する．
>> b = [ones(length(x), 1)  x  x .* x] \ y;
>> x3 = (1 : 0.1 : 5)';
>> y3 = [ones(length(x3), 1)  x3  x3 .* x3] * b;

近似曲線を表示する．
>> plot(x, y, 'o', x3, y3);



4. 主成分分析


例題：国語・英語・数学・理科のテスト結果を主成分分析してみよう．

まず，データを読み込む．
>> load -ascii tensuu

※ tensuu という名のテキストファイルが読み込まれ，同じ名前の変数に格納される．

棒グラフとして表示してみる．
>> bar(tensuu);

※ 4科目とも点数の高い人，4科目とも点数の低い人，国語・英語はいいけど
数学・理科は悪い人，その逆に分かれる．すなわち，国語と英語，数学と理科
には，点数の間に相関性がある．

主成分分析の前に，平均が0，分散が1になるように正規化する．
>> N = size(tensuu, 1);
>> m = mean(tensuu);
>> s = std(tensuu);
>> tensuu2 = (tensuu - repmat(m, N, 1)) ./ repmat(s, N, 1);

主成分分析を行う．
>> [w, Z, s2, t2] = princomp(tensuu2);

第2主成分までを棒グラフで表示する．
>> bar(Z(:, 1 : 2));

各主成分の意味を考える．
>> w
w =
   -0.4873   -0.5273   -0.4990    0.4853
   -0.5105   -0.4740    0.5387   -0.4738
   -0.5083    0.4807   -0.5041   -0.5063
   -0.4935    0.5159    0.4547    0.5326

※ 第1列が第1主成分の因子負荷量を表す．この場合，各科目の点数が総合的に
低いほど，第1主成分は高い値を示す．同様に，第2主成分は，国語・英語の点数
が低いほど高い値を示し，数学・理科の値が高いほど高い値を示す．すなわち，
第1主成分が各科目を総合したテストの成績（ただし低いほどよい），第2主成分
は，文系的か理系的かを表している．

累積寄与率を見る．
>> pareto(s2 / sum(s2));

※ 棒グラフが各主成分の寄与率，折れ線グラフが累積寄与率である．第2主成分
まででほとんど100%の情報を表現できていることが分かる（情報を半分に圧縮）．



5. MATLAB らしいプログラミング


MATLAB は，同じ振る舞いをするプログラムでも，書き方によって処理速度が
大きく異なる．

例1 for 文と行列

書き方1
>> S = 0; for i = 1 : 1000, S = S + i * i; end

書き方2
>> a = 1 : 1000; S = sum(a .* a);

※ 行列を直接操作できる組み込み関数は最適化されているのに対し，
   for 文などの繰り返しは，インタプリタ言語のため遅い．

例2 メモリの確保

書き方1
>> for i = 1 : 10000, x(i) = i; end

書き方2
>> y = zeros(1, 10000); for i = 1 : 10000, y(i) = i; end

※ 書き方1は，繰り返しの度に行列のサイズが大きくなるため，
   メモリ確保をやり直しているために遅い．

所要時間を測るには，tic（計測開始）と tic（計測終了）を使う．
たとえば，
>> tic; S = 0; for i = 1 : 1000, S = S + i * i; end, toc;



6. M-ファイルプログラミング


MATLAB では，プログラムは m ファイルに書く．
ファイル名と関数名を一致させる．
入力引数はもちろん，出力引数を複数指定することができる．
関数内で，自由にローカル変数を使用することができる．

例：ファイル名「calc.m」
=================================================================
function [m, s] = calc(a, b, c)

if nargin < 3                     % 入力引数の個数が 3 未満なら
  error('引数が足りません');      % エラー
end

m = (a + b + c) / 3;
m2 = (a * a + b * b + c * c) / 3;
s = sqrt(m2 - m * m);
=================================================================



7. ヘルプ


>> help funcname
で，funcname という名の関数のヘルプが表示される．

>> helpdesk
とすると，オンラインヘルプを見ることができる．



8. 参考文献


[1] MATLAB 付属マニュアル："Getting Started with MATLAB"．
[2] 上坂吉則："MATLABプログラミング入門"，牧野書店，2000．
[3] 小林一行："MATLAB活用ブック"，秀和システム，2001．
その他，多数の本が市販されている．また，Web 上にもさまざまな情報がある．