ベクトル場の可視化

最近天気予報で見る風の可視化がかっこよかったので，真似して研究で使うベクトル場の可視化をやってみたところ，細部がわかりやすくていい感じでした．

汚いですがコードはMATLABのコードは以下のような感じです．
動画のベクトル場を生成するコードはすごく長いので，単純な円を描く流れにしてあります．コードの最初の部分でお好みのベクトル場を指定してみてください．

MATLABで関数型プログラミング入門

先日まで Northeastern University の Robust Systems Lab で訪問研究員をやっていました．その時に関数型プログラミング入門のチュートリアルを行ったのですが，わりと好評だったと思うので，スライドを日本語化したものを公開します．

関数型プログラミング入門 for Matlab ユーザー from Ichiro Maruta

いわゆる手続き型の言語にも関数型言語由来の機能が追加されるようになって久しく，研究のためのコーディングでもこれらの機能は有用です．しかし，これらの機能の使い道がピンとこないまま，あるいは存在に気づかないまま活用できずにいる学生が多いように感じています．

普通であれば「関数型プログラミングについて勉強しておいて」と学生に伝えれば済むのですが，典型的な関数型プログラミングの解説は典型的な学生に関数型プログラミングの有用性を理解させられるものではなく，学生にとっての体感的な学習コスパ（=有用性/学習コスト）が著しく低いために身につかないように思われます．僕の感覚では，「手続き型言語に追加されている関数型言語由来の機能でちょっと楽をしたい」程度の低いこころざしの範囲であれば，関数型プログラミングの学習コスパは極めて高く，これは適切な教材の不在によって生じている損失のように思われました．

というわけで，このチュートリアルを行うにあたっては，学生が興味を持ちそうな卑近な（日々の宿題を解くのに役立ちそうな）例題を選び，利便性の積極的なアピールを行うアプローチを試みています．大学生という聴衆に特化しているので役に立たない場合も多々あるかと思いますが，何かの参考になりましたら幸いです．

MATLAB で figure にフォーカスを奪われない方法

MATLAB で figure(2) のように figure 関数を実行すると Figure 2 のウィンドウにフォーカスが移行し，キーボード操作は Figure 2 のウィンドウが受け付けるようになります．若干再現性に欠けるところもありますが，仕様ではそうなっているようです．
通常はこれで問題ないのですが，長時間の計算を行っていて途中経過を図で表示している場合，別の作業をしている最中に邪魔をされて不便な場合があります．
この場合

figure(2)

のかわりに

set(0,'CurrentFigure',2)

とすれば Figure 2 が描画の対象になりますがキーボードのフォーカスは移動しません． ただし， Figure 2 が存在しない場合エラーになるので，うっかりウインドウを消してしまったときのために

try
  set(0,'CurrentFigure',figidx)
catch
  figure(figidx);
end

としておくとより安全だと思います．

MATLABでQhullを使う

Qhull （http://www.qhull.org/） は計算幾何に関する有用なライブラリーで，MATLABの幾つかの関数，delaunayn や voronoin は Qhull を利用しています．

しかし，これらの関数は Qhull が持つすべての機能を実装しているわけではなく，例えばドロネー分割における法線や隣接単体の情報を出力するオプションはMATLABでの実装ではサポートされていません．

本来なら Qhull を呼び出すmexファイルを書くのが正攻法ですが，面倒なのでコンパイル済みのQhullの実行ファイルを呼び出す関数を作ってみました．

僕が欲しかったのはドロネー単体分割時の隣接単体の情報だったので，以下の関数はドロネー単体の情報に加えて単体間の隣接行列を返す仕様になっています．

function [tes,adj]=mydelaunayn(X)
%MYDELAUNAYN N-D Delaunay triangulation with adjacency matrix
%   [tes,adj]=mydelaunayn(X)
%   tes: sparse matrix which represents tessellation
%      tes(i,j)==1 if j-th simplex has vertex X(i,:)
%   adj: adjacency matrix among Delaunay simplices
%      adj(i,j)==1 if i-th and j-th simplices has common hyper plane
%
%   NOTE:
%     This routine requires excutable of Qhull.
%     http://www.qhull.org/


%% Compose input
infile=tempname;
fid=fopen(infile,'w');
d=size(X,2);
N=size(X,1);
fprintf(fid,'%d\n',d);
fprintf(fid,'%d\n',N);
fprintf(fid,strcat('%e',repmat('\t%e',1,d-1),'\n'),X');
fclose(fid);

%% Invoke Qhull
outfile=tempname;
opts={'Qt','Qbb','Qc'};
if d>3
    opts{end+1}='Qx';
end
opts=strcat(opts,{' '});
system(sprintf('qhull d s i %s Fn < %s > %s',cat(2,opts{:}),infile,outfile));

%% Read Results
fid=fopen(outfile,'r');

% compose tessellation matrix
Ntes=fscanf(fid,'%d',1);
tesmat=fscanf(fid,'%d',[d+1,Ntes])'+1;
idx=1:Ntes;
tes=sparse(reshape(tesmat',[],1),reshape(idx(ones(d+1,1),:),[],1),1,N,Ntes);

% compose adjancy matrix
Nadj=fscanf(fid,'%d',1);
adjout=fscanf(fid,'%d',[d+2,Nadj])';
assert(all(adjout(:,1)==d+1));
assert(Nadj==Ntes);
adjmat=[reshape(idx(ones(d+1,1),:),[],1),reshape(adjout(:,2:end)',[],1)];
adjmat(adjmat(:,2)<0,:)=[];
adj=sparse(adjmat(:,1),adjmat(:,2)+1,1,Ntes,Ntes);

%% Cleanup
delete(infile);
delete(outfile);
end

上の関数の実行にはパスに Qhull の実行ファイルがあることが必要で，これは http://www.qhull.org/ からダウンロードできます．他の用途への改造も比較的容易だと思うので，よければ参考にしてください．

MATLAB の無名関数で複数の文を実行する方法

通常 MATLAB では関数を定義する際に別ファイルを作成せねばならず，これは結構面倒です．

無名関数を使うと

f = @(x) x^3-x;
ezplot(f,[-2,2])

のように簡単な関数をその場で定義できますが，複数のステートメントを含む関数を書くことは簡単ではありません．
ここでは複数の文を含む無名関数を定義するために僕が使っているテクニックを紹介します．

複数の文を含む無名関数を書く場合，「任意の数の式を引数とし，最後の式の値を返す関数 progn 」を定義しておくと便利です．

progn = @(varargin) varargin{end};

f = @() progn(...
    figure(1),    ...
    ezplot(@sin), ...
    figure(2),    ...
    ezplot(@cos) );

f();

上の例では複数の図のプロットを行う関数を定義しています．

ただし，この方法では値を返さない関数や文を並べることができません．たとえば hold on などがこれに該当します．どうしてもこれらの関数を使いたい場合，邪道ですが evalc を使うことで問題を回避できます．

progn = @(varargin) varargin{end};

f = @() progn(...
    figure(1),      ...
    ezplot(@cos),   ...
    evalc('hold on'), ...
    ezplot(@sin));

f();

また，MATLABの無名関数はローカル変数を持つことができませんが，evalc と組み合わせると引数に対応する変数をローカル変数として利用することが可能です．

progn = @(varargin) varargin{end};

% 二番目に大きい固有値を返す関数を定義
eig_2nd = @(x) feval(@(x,e,idx) ...
    progn( ...
        evalc('e=eig(x);'), ...
        evalc('[~,idx]=sort(abs(e));'), ...
        e(idx(2))), ...
    x,0,0);

A=[4,0,0;
   0,1,0;
   0,0,2];
eig(A)
eig_2nd(A)

見てのとおり邪道で，書き捨てのプログラム以外で使うことは憚られますが，これを使うとかなり広範な関数をその場で定義できます．

MATLAB の無名関数内で分岐を行う方法

MATLAB には一応無名関数がありますが，関数型のプログラミング言語ほど周辺のユーティリティが充実していないので不便なことが結構あります．最近は MATLAB での無名関数の使い方もだいぶわかってきた気がするので，ちょくちょく小技を紹介していこうかと思います．

今回は僕が使っている分岐を実現する方法を紹介します．MATLABにはC言語のような３項演算子や条件分岐関数が無いために，無名関数内では手軽に分岐を行えません．

分岐を簡単に実現するひとつの方法は，あらかじめユーティリティ関数を作っておくことです．
たとえば僕は以下のような関数 switchfunc をパスに入れています．

function [ out ] = switchfunc( varargin )
%%SWITCHFUNC 条件分岐を簡潔に書くための関数
%  SWITCHFUNC( cond_1, value_1,...
%              cond_2, value_2,...
%                     :
%               true, default_value)
%  cond_1,cond_2, の中で 最初に true の値をとるものを
%  cond_k とすると，value_kを値として返します．
%  全ての条件が満たされない場合エラーを返します．

for k=1:2:nargin-1
    if varargin{k}
        out = varargin{k+1};
        return
    end
end
error('switchfunc: no match');
end

これを使うと，たとえばxの値に応じて切り替わる関数を以下のように書けます

f = @(x) ...
    switchfunc( -2 <= x && x <= -1,  x+2, ...
                -1 <  x && x <   1,    1, ...
                 1 <= x && x <=  2,    x);

ezplot(f,[-2,2]);

結果得られる関数 f  のグラフは↓のようになり

自然な表記で切り替えを持つ関数を作れていることがわかります．

しかし，実はこの switchfunc による分岐と 通常の if 文による分岐には大きな違いがあります．
それはswitchfunc を使った場合，実際の条件分岐を行う前に全てのケースについての値が計算されているという点です．
先ほどの例題のように，値の計算に副作用がない場合には速度が多少遅くなる程度の害しかありませんが，plot のような関数を条件に応じて切り替える用途にはそのままでは使えません．

通常のif 文のように条件と合致した時だけ値の評価を行いたい場合は以下のように一工夫します．

swplot = @(num) ...
            feval( ...
                switchfunc(...
                    num==1, @() ezplot(@sin),...
                    num==2, @() ezplot(@cos)));
figure(1), swplot(1)
figure(2), swplot(2)

ポイントはswitchfuncで切り替える値を無名関数にしていることで，帰ってきた結果のみを feval で実行することによって通常のif文と同様に条件に対応する値のみが評価されます． 

MATLAB の無名関数は副作用がない純粋関数しか含まないことを前提に設計されているので，無名関数のなかで完結するようプログラムを書くと自然にLISPのような関数型言語にちかい構造のプログラムになる気がします． MATLAB は元来このようなスタイルを前提とする言語ではないので表記は複雑になりがちですが，関数型言語で使われるテクニックの中にもMATLABで有用なものがあるかもしれませんね．

MATLABでサイズを指定してイメージファイルを出力する方法

MATLABの図を画像ファイルとして出力する場合，デフォルトでは印刷時と同様に適当なサイズの画像が保存されます．

しかし，プレゼン用画像の作成などでサイズを1024×768といった特別な値にしたい場合がしばしばあります．毎回やりかたを忘れて苦労するので，ここに方法をメモしておきます．


figure の PaperPositionMode を auto に設定し，print コマンドで '-r0' を付けるのがポイントのようです．ファイルを出力する際に幅が1ピクセル増える現象が発生しますが，バグなのか仕様なのか常に再現するのかもよく調べていません．バージョンによって挙動が異なる可能性もありますので注意してください．

MATLAB MEXでCUDAを使う

GPU 関連のプログラミングでは大抵、C言語＋αを使ってプログラムを書くみたいですが、データの処理や結果の可視化は普段使っている MATLAB のほうが楽です。

そこで、 CUDA な MEXファイルを作って MATLAB で使う方法を試してみました。
結構簡単ですが、ハマるところがいくつかあったので記録を残しておきます。

現在僕が使っている環境は

• Windows 7 (64bit)
• MATLAB 2010a (64bit)
• CUDA 3.1 (64bit)
• Visual Studio 2008 Professional

なので、参考にするときは注意してください。以下順にハマったところを挙げていきます。

64bitのコンパイラをインストールする

これはCUDAに限りませんが、mexファイルのコンパイルにVisual C++を使う際、64bitの環境では64bit用のコンパイラをインストールしておく必要があります。VS2008では64bit用のコンパイラはデフォルトではインストールされません。また、Express Edition には64bit用のコンパイラが含まれていないようです。

最新版のMATLAB plug-in for CUDAを入手する

実は一番ハマったのがここです（笑）

このページで右側のフォームに住所氏名などなどを記入するとダウンロードできるのが最新版です。

こっちのページでダウンロードできるのは古いバージョンです。（古いMATLABだとこっちを使うべきなのかもしれません。）最初新しいバージョンは非公開と勘違いしていました・・・

ちなみに日本語版のサイトには新しいバージョンへのリンクすら存在してないようです。

nvmex.mを修正する

MATLAB plug-in for CUDA と一緒に入手できるガイドラインにも書いてありますが、使用環境に応じてnvmex.mを修正する必要があります。ぼくの環境では

Host_Compiler_Location = '-ccbin "C:\Program Files (x86)\Microsoft Visual Studio 9.0\VC\bin"';

としました。ソースファイルで変更しろと書かれているのはここだけですが、実際にはあといくつか修正が必要で、下のほうの nvccCommandLine を作る部分で、インクルードファイルのパスを指定している部分にダブルクォーテーションが抜けてるバグがあります。Linux環境では問題が発生しなくて見逃したのでしょうか？あと、ついでにCUDAのサンプルでよく使われているcutilを使うためのヘッダファイルのパスも追加しておきます。

    ' -I"' matlabroot '\extern\include" ' ...
    ' -I"C:\ProgramData\NVIDIA Corporation\NVIDIA GPU Computing SDK\C\common\inc"'...

これに対応してリンカのオプションも mexCommandLine の作成の部分に追加しておきます。

mexCommandLine = ['mex (''' filename '.o' ''', ''' '-L' CUDA_LIB_Location ''', ''-LC:\ProgramData\NVIDIA Corporation\NVIDIA GPU Computing SDK\C\common\lib''' ', ''-lcudart'', ' '''-lcutil64'')'];

あと、cutilを使っているプログラムを実行する際には cutil64.dll をパスが通っている場所に置かなければならないので注意してください。


とりあえず、ここまでの修正でサンプルが動くところまではいけたはずですが、
何か思い出したことがあったらまた追記します。

MATLAB が出力する図を PowerPoint 向けに自動調節

MATLAB で描いた図やグラフを PowerPoint に貼り付けるときに
一番簡単な方法はFigure のメニューから [編集]-[Figureのコピー] を使って
PowerPoint 上でペーストすることだと思います．

しかし， MATLAB がデフォルト設定で作成するグラフには
以下のような問題があります．

• プレゼン用には字が小さい
• グリッドが出ない
• 無いよりある方が無難だと思うのだけど・・・
• コピー＆ペースト時にサイズが変わる
• グリッドまで変えられるので意図しない結果になることが多いですね
• 線が見にくい
• 細い＆プロジェクターで見辛い色が使われます

丁寧に設定すればこれらの問題を解決することは可能ですが，
そこそこ面倒ですし，ミーティングなどで MATLAB による実演を行う際に
これらを逐一設定するのはテンポを損ねます．

というわけで，これらの設定を自動化するプログラムを作りました．
このプログラムは MATLAB のデフォルト設定を書き換えて，
以降に描かれるグラフをプレゼン用に見やすくします．

いまのところ MATLAB を再起動する以外に元に戻す方法はありません．

ちなみに得られるグラフは以下のようになります．

使用前

使用後

コード

function [] = prettify_figure_presen()
%PRETTIFY_FIGURE_PRESEN グラフのフォントをプレゼン用に最適化
% この関数を呼び出すと
% 以降に描かれるプロットの見た目がプレゼン向けになります．
% 元に戻すにはMATLABを再起動するしかないです．
%
% ※ ボード線図は別口(prettify_bodeplot_presen参照)
%    なので注意してください．


%% フォント
set(0,'defaultAxesFontSize',20);
set(0,'defaultAxesFontName','メイリオ');
set(0,'defaultAxesFontWeight','demi'); % normal/demi/bold
set(0,'defaultTextFontSize',20);
set(0,'defaultTextFontName','メイリオ');
set(0,'defaultTextFontWeight','demi');

%% 線の幅
set(0,'defaultAxesLineWidth', 1.5); % 軸
set(0,'defaultLineLineWidth', 2); % プロット

%% プロットの色
clorder = ...
    [  0,   0, 255; % 青
       0, 128,   0; % 緑
     255,   0,   0; % 赤
     204,   8, 204; % 紫
     222, 125,   0; % 茶
       0,  51, 153; % 紺 (青と区別しづらい，注意！)
      64,  64,  64];% 濃い灰色  
set(0,'defaultAxesColorOrder',clorder/255);

%% 最初からgrid on & box on
set(0,'DefaultAxesXGrid','on');
set(0,'DefaultAxesYGrid','on');
set(0,'DefaultAxesBox','on');

%% コピペ時にサイズを変更しない
set(0,'DefaultFigurePaperPositionMode','auto');

end

2010/2/17 追記
17行目の間違いを修正しました＆reset(0)で設定を初期化できるそうです．
atsu-kan さんありがとうございます．

MATLAB で半透明の plot を作る

データが大量にあって，グラフ上で密度を表現したい場合などで
MATLAB で半透明の線を書きたいことがたまにあります．

なぜかMATLABのplotで生成される線は半透明をサポートしていないので
半透明をサポートしているpatchで線を描く関数を作ってみました．

function [ p ] = pplot(x,y,z)
%PPLOT Draw line as patch

if nargin==2
    z=0;
end

x=reshape(x,[],1);
y=reshape(y,[],1);

p=patch([x;flipud(x)],[y;flipud(y)],z*ones(2*size(x,1),1),'b');

end

使うときは以下のように使います．

x=0:1e-2:pi/2;
clf
hold on
for k=linspace(0,1,25)
    p=pplot(x,(1-k)*sin(x));
    set(p,'EdgeAlpha',0.1);
    set(p,'EdgeColor','b');
    set(p,'LineWidth',2);
    p=pplot(x,k*cos(x));
    set(p,'EdgeAlpha',0.1);
    set(p,'EdgeColor','r');
    set(p,'LineWidth',2);
end

得られる結果はこんな感じ
このように，半透明の線を描くことができ，線が集中している部分は濃い色で表示されます．
ただし，複雑な線を大量にプロットするとすごく遅いので注意してください．

論理インデックス

最近「論理インデックス(logical indexing)」というMATLABの便利機能に気づきました．
こんなふうに使います．

>> a=[1,2,3,4]
a =
    1 2 3 4
>> a(a>2)
ans =
    3 4
>>

配列のある条件を満たす部分だけ抽出する時に便利です．
通常のインデックス a(find(a>2))とかよりも高速に処理されるらしいです．

SyntaxHighlighter2でMATLABソースを表示 in Blogger

以前SyntaxHighlighterをMATLAB対応にしました．その後，ファイルを置いてあるサーバーの更新があったり，SyntaxHighlighterがバージョン2になったりで使えない状態で放置していたので，改めて作りなおしてみました．

BloggerでSyntaxHighlighter2を使う

まず最初にBloggerでSyntaxHighlighter2をつかうときの設定ですが，テンプレートに

<!--ここから -->
    <script src='http://alexgorbatchev.com/pub/sh/2.0.278/scripts/shCore.js' type='text/javascript'/>
    <script src='http://alexgorbatchev.com/pub/sh/2.0.278/scripts/shBrushPlain.js' type='text/javascript' />
    <script src='http://alexgorbatchev.com/pub/sh/2.0.278/scripts/shBrushXml.js' type='text/javascript' />
    <script src='http://alexgorbatchev.com/pub/sh/2.0.278/scripts/shBrushCss.js' type='text/javascript' />
    <script src='http://alexgorbatchev.com/pub/sh/2.0.278/scripts/shBrushJScript.js' type='text/javascript' />
    <script src='http://maru.bonyari.jp/storage/shBrushMatlab.js' type='text/javascript'/>
<!--ここまで -->
</head>
...
<!-- end outer-wrapper -->
<!--ここから -->
<script type='text/javascript'>
//<![CDATA[ 
    SyntaxHighlighter.config.bloggerMode = true;
    SyntaxHighlighter.all();
//]]>
</script>
<!--ここまで -->

のように追加するといいようです．config.bloggerModeの設定がないと改行の位置にBRタグの出来そこないが出現します．

MATLABのソースをSyntaxHighlighterで成形する

•  shBrushMatlab.js

【使い方】

上記のファイルを通常の”shBrush*.js”と同じように読み込んで，class=”brush:matlab”で使えます．上のサーバーに置いてあるのを読み込んでも構いませんが，無保証自己責任でお願いします．ちなみにSyntaxHighlighter本体も本家でホスティングされているようですね．（リンク）
上の設定例はこのホストを使う設定にしてます．

MATLABで並列計算(PSO)

素のMATLABは並列計算に弱いようで，最近のデュアルコアのPCだと，プロセッサ使用量が５０％で頭打ちします．

ただ，最近学校の MATLAB ライセンスに Parallel Computing Toolbox が加わって一部関数が置き換えられたためか，単にバージョンアップの恩恵なのかは知りませんが，一部の計算ルーチンがマルチコア用にコンパイルされているらしく，計算によっては50%を超えるプロセッサ使用率を見ることもできるようになりました．

でも依然として，大半の処理は片方のコアのみで処理されるようです．
本来は並列計算できる部分は処理系が自動的に判別して最適に処理してくれるのが理想なんですが，なかなか実現されないですね．また，最近はコア１つ当たりの性能向上は収束気味で，今後はコアの数が増える方向にPCは進化するようなので，今後MATLABでPCの進化の恩恵を享受するためには，並列計算のスキルを身につけておく必要がありそうです．

MATLAB で Parallel Computing Toolbox が使える場合，単一マシンで並列処理をするのはとても簡単で

matlabpool 2 %デュアルコアの場合
parfor idx=1:100
   % 並列化する処理
end
matlabpool close

とやれば並列に処理されます．matlabpool の後の数字は生成されるスレッド数です．

結構並列化のオーバーヘッドがあるようで，単位処理の計算量があまりに小さい場合，かえって遅くなることもあるようです．結果が芳しくないときはプロファイラでチェックしてみるべきですね．

Particle Swarm Optimization (PSO) なんかは並列計算の恩恵を最大限に受けられるアルゴリズムなので，目的関数の評価に時間がかかる問題であれば，処理速度がほぼコア数倍になります．（といってもまだデュアルコアでしか試してませんが・・・）

もうちょっと手間をかけると複数のマシンをつないでクラスタリングとかもできるようなので，人がいない夜中とかに実験室のマシンを束ねて試してみようかな〜

・・・・と思ったらクラスタリングには MATLAB Distributed Computing Server っていう別製品が必要なようですね．

また，Parallel Computing Toolbox ではスレッドの数も４つまでに制限されてるようです．

つまり現状ではクアッドコアのXeon W3570以上のCPUを使っても性能向上は期待できないということかな？

値段を考えると2.66GHzのXeon W3520か，近々更新される見込みなCore i7あたりが計算用マシンのCPUとしてはベストになるのだろうか・・・

MATLABでフィードバック系のシミュレーション

たまにうっかり忘れるのでメモ．

MATLABでプラントPと制御器Kが構成するフィードバック系のシミュレーションをするとき，

>> s=tf('s');
>> P=1/s;
>> K=1/(1+s);
>> step(P*K/(1+P*K))

とかやるのはNG．正しくは

>> step(feedback(P*K,1))

とする．一見結果は同じなのだけど

>> P*K/(1+P*K)

Transfer function:
s^2 + s
-----------------------
s^4 + 2 s^3 + 2 s^2 + s

>> feedback(P*K,1)

Transfer function:
1
-----------
s^2 + s + 1

>>

という違いがあり，高次のシステムを使うときになって初めて問題が顕在化したりしてよろしくない．
システムの割り算はMATLABでは逆システムの計算と直列接続で計算されるので，この場合P*Kと1+P*Kが共通の極を持っているにも関わらず，これらが相殺されないことによって問題が起こる．
ちなみにseriesとparallelを*と+で書くのは，少なくともSISOの場合には問題ない（中でやってる処理は同じ）．MIMOの場合は未確認．

MATLABに限らないけど，強力なツールを使うのは楽なようで難しい・・・
ソースを全部追うのは大変だし・・・

MATLABで可観測正準形/可制御正準形

制御の教科書では必ず出てくる可観測正準形/可制御正準形ですが，数値計算上は大変扱いづらく，実際に数値的にこの形式にあえて変換することなんて試験かレポート以外ではあまりないと思います．
そのあたりの事情を反映してかMATLABにはこれらの形式に変換する関数というのは存在していないようです．
ctrbf/obsvfという関数があります．

しかしまぁ研究の都合上これらの正準形を計算したい場合もあるわけで，そんな場合に使うルーチンを作りました．１入力１出力の場合にしか使えません．1入力１出力の場合だと伝達関数形式に変換した時点でこれらの形式は単純に決まるわけですが，復習をかねて教科書的な方法で計算しています．

可観測正準形

function [ csys ] = obsvcanon(sys)
%OBSVCANON Compute observability canonical form

[n,d]=tfdata(tf(sys),'v');
W=hankel(fliplr(d(1:numel(d)-1)));
G=ss(sys);
U0=obsv(G);
S=W*U0;
csys=ss(S*G.a*inv(S),S*G.b,G.c*inv(S),G.d);

end

可制御正準形

function [ csys ] = ctrbcanon(sys)
%CTRBCANON Compute controllability canonical form

[n,d]=tfdata(tf(sys),'v');
W=hankel(fliplr(d(1:numel(d)-1)));
G=ss(sys);
UC=ctrb(G);
S=UC*W;
csys=ss(inv(S)*G.a*S,inv(S)*G.b,G.c*S,G.d);

end

ダウンロード

• obsvcanon.m
• ctrbcanon.m

MATLABのグラフに目盛りを追加

MATLABでグラフを作るとき，自動でつく目盛り以外に特定の値の目盛りがほしい場合があります．真値のところに目盛りがほしいとか収束先のところに目盛りがほしいとかそんなときです．

じつはこれをプログラムからやるのは結構面倒だったりします．
なので専用の関数を作ってみました．

こんなグラフに

>> addxtick(2/pi)

とすると

こんな風に目盛りがつきます．
ラベルもいじりたいときは

>> addxtick(pi/2,'pi/2')

と入れれば

となります．便利ですよ．
短いプログラムなのでコピーしてaddxtick.mという名前で保存して使ってください．

function [] = addxtick( tk, tklbl )
%ADDXTICK --- Add new tick to x-axis
%   Ex.    addxtick(3.14)
%   Ex2.  addxtick(3.14, 'Pi')

if nargin==1
    set(gca,'XTick',unique(sort([get(gca,'XTick'), tk ] )));
else
    temp = cellstr(get(gca,'XTickLabel'));
    temp = {temp{:},tklbl};
    [temp2,idx] = sort([get(gca,'XTick'),tk]);
    [temp2,idx2]=unique(temp2);
    idx=idx(idx2);
    set(gca,'XTick',temp2);
    set(gca,'XTickLabel',temp(idx));   
end

end

あと，ほとんど同じですがaddytick.mも置いておきます．

function [] = addytick( tk, tklbl )
%ADDYTICK --- Add new tick to y-axis
%   Ex.    addytick(3.14)
%   Ex2.  addytick(3.14, 'Pi')

if nargin==1
    set(gca,'YTick',unique(sort([get(gca,'YTick'), tk ] )));
else
    temp = cellstr(get(gca,'YTickLabel'));
    temp = {temp{:},tklbl};
    [temp2,idx] = sort([get(gca,'YTick'),tk]);
    [temp2,idx2]=unique(temp2);
    idx=idx(idx2);
    set(gca,'YTick',temp2);
    set(gca,'YTickLabel',temp(idx));   
end

end

xPC Target で CNT24-4D(PCI)H を使う

大学のMATLABのライセンスが統合された関係で，xPC Targetが使えるようになりました．研究室は古いPCがあまりがちなので，これを実験機に回せば安上がりにラピッドプロトタイピングな環境を構築できそうです．

しかし，xPC Targetに対応してるIOボードというのは意外と少なくて

http://www.mathworks.com/products/xpctarget/supportedio.html

に列挙されてるやつだけです．しかも互換品で型番が違うだけみたいなボードもPCIボードの VendorID / DeviceID で完全にはじいてしまうようです．研究室にあるContecのボードのいくつかは上のリストに挙げられてるんですが，エンコーダ用のカウンタは，すでに廃品種の CNT24-4D(PCI) しか対応して無くて研究室にある現行の CNT24-4D(PCI)H はDeviceIDが違うので認識されません．

幸いなことに xPC Target の IOボードのドライバのコードは書き換えることが出来て，mexなコンパイラの設定をすれば（とりあえず Visual C++ 2008 Express editionが使えた・・）コンパイルして使うことが出来ます．

とりあえず，CNT24-4D(PCI)Hに関してはCNT24-4D(PCI)用のソースの簡単な書き換えだけで対応できることを確認しました．二慣性共振系の実験装置に早速つないで使ってみましたがなかなかいい感じです．実験は本業じゃないし，簡単に実験できるのはいいことですね〜

実は研究室には Interface社のボードの方が多いので，ドライバを書いて使えるようにしたいところです．DAとADくらいならドライバを書くのもそんなに難しくなさそうですしね．というかInterfaceさんが出してくれたりしないかなぁ・・

MATLAB対応のSyntaxHighlighter

！！　ここの内容はもう古いです．新しいバージョンを作りました ！！

MATLABのソースファイルに対応したBloggerで使える整形ツールを探してたんですが，無いようなのでSyntaxHighlighterを改造してみました．

やっぱり人の書いたプログラムもたまには見てみるものですね．勉強になります．
ちなみに追加したコードはこんな感じ．

dp.sh.Brushes.Matlab = function()
{
var keywords = 'break case catch classdef continue else elseif ' +
   'end for function global if methods otherwise ' +
   'persistent properties return switch try while';
this.regexList = [
{ regex: new RegExp('%.*', 'gm'),    css: 'comment' },   // one line comments
{ regex: dp.sh.RegexLib.DoubleQuotedString,     css: 'string' },   // strings
{ regex: dp.sh.RegexLib.SingleQuotedString,     css: 'string' },   // strings
{ regex: new RegExp(this.GetKeywords(keywords), 'gm'),  css: 'keyword' }   // c# keyword
];

this.CssClass = 'dp-c';
this.Style = '.dp-c .vars { color: #d00; }';
}

dp.sh.Brushes.Matlab.prototype = new dp.sh.Highlighter();
dp.sh.Brushes.Matlab.Aliases = ['matlab', 'MATLAB'];

Firefoxでは折り返しがうまく行きませんね．そのうち改善されるかな？
PSOのプログラムとかをこれを使って密かに追加したので良かったら見てくださいませ．

MATLAB用PSOクラス

(English version is HERE.)

MATLABの2008aからオブジェクトの新しい書き方が導入されました．

あんまりオブジェクト指向なプログラミングの経験はないのですが，新しい書き方はだいぶ簡単で使いやすいと思います．

というわけで，とりあえず PSO のルーチンを2008a仕様に書き換えてみました．

pso.m -- PSOクラス
psotest.m --- pso.mの使用例，ベンチマーク問題を解きます．psotest2.m --- pso.mの使用例，簡単な制約付き最適化問題を解きます．

使い方はpsotest2.mを見ればだいたいわかるはず‥

アルゴリズムに関する詳細はIFAC2008の発表予稿，もしくはAutomaticaを見られる環境であれば

• Ichiro Maruta, Tae-Hyoung Kim, Toshiharu Sugie,
  Fixed-structure H∞ controller synthesis:
  A meta-heuristic approach using simple constrained particle swarm optimization,
  Automatica, In Press, Corrected Proof,
  Available online 1 January 2009

を参照してくださいませ．

MATLABで無名関数

ずっと先輩のコードとかを見て何となく
MATLABのプログラムを書いていたんだけども、
改めてチュートリアルとか読んでみると知らなかった機能が結構ある。

中でも無名関数があったのを知らなかったのは大分損してたかも。
ちなみに使い方はこんな感じ。

>> f=@(x,y) x^2+y^2

f =

  @(x,y)x^2+y^2

>> f(1,3)

ans =

  10

>>