メモです。余裕のあるときに後でリファインします。 前提 微小回転行列 クォータニオンによる角速度の積分 MATLABで比較プログラム 問題設定① 微小回転行列の作成 クォータニオンの計算 問題設定② 自前関数 参考文献 前提 座標系が角速度 で回転しているとする。 微小回転行列 本来三次元回転SO3は非可換であるが，微小回転を考えればこれを可換にすることができる。1 先程の角速度で微小回転行列というものが考えられており， とすれば以下のような回転行列で表せる。 これは以下に示すcross product表記を用いれば，と表すことができる。 wikipediaより引用：歪対称行列 クォータニオンによる角速度の積分 クォータニオン便利計算ノートによればクォータニオンによる積分は以下の通り。 一番下の式には1/2が入っていないと思う。見落としがなければ。 先程のcross productを用いて表すなら， $$ \dot{q}=\frac{1}{2} \begin{pmatrix} 0 & -\omega^\top\\ \omega & \omega_{\times} \end{pmatrix} q $$ と書いたほうがスッキリするだろう。 MATLABで比較プログラム P.I.corke先生のRobotics Tool boxを使用。課金できる方は自前でToolboxを買うと良い。 なお，2019年以降に東大に在籍している学生は無料で全ライセンス使えるのでぜひ活用するとよい。 petercorke.com 昔は全機能使えたと思ったのですが今入れたらいくつかの機能が削除されていたような？？？ 問題設定① 最初の角度をオール0とする。同次変換行列を作っておく。 orig_eul = [0,0,0]' orig_tr = rpy2tr(orig_eul') omegaを適当に定め，10ms積分するとする。 omega = [0.1,0.2,0.3]; dt = 0.01; 微小回転行列の作成 cross productを計算するskewという関数を用いて簡単に作成可能。 R1=skew(omega*dt)+eye(3) TR1 =r2t(R1)*orig_tr TR1 = 1.0000 -0.0030 0.0020 0 0.0030 1.000...

座標変換などで手計算をあまりしたくないのでsymsをうまく使ってこれを回避する方法をメモします。 Image Jacobianの導出 時間変化する変数を定義 時間で偏微分 特定の数式の代入による式整理 小手先の代入テク 多項式から 行列形式への変形 まとめ prettyコマンドの活用 参考文献 Image Jacobianの導出 カメラの座標系から見て位置（X,Y,Z）にある物体を考えます。焦点距離fとすると，以下の関係が成り立ちます[1]。 カメラの撮像幾何[1] カメラが6軸の速度・加速度を持つときに画像内の点がどのように動くのか，その関係性を表す行列を導出します。 何故かうさぎをアクチュエータに使うPeter先生。かわいい。 時間変化する変数を定義 時間tに関する関数として位置の各変数を定義できます。[2]の橋本先生の論文を例にとって見てみましょう。 ピンホールカメラモデルによる撮像の式（ただしf=1として省略） syms t X(t) Y(t) Z(t) image_p = [X(t); Y(t)]/Z(t) 計算結果： image_p = (X(t))/Z(t) (Y(t))/Z(t) 時間で偏微分 tで偏微分します。 商の微分 diff_p = diff(image_p,t) こうなります。 diff_p = diff(X(t), t)/Z(t) - (X(t)*diff(Z(t), t))/Z(t)^2 diff(Y(t), t)/Z(t) - (Y(t)*diff(Z(t), t))/Z(t)^2 特定の数式の代入による式整理 先程の式のXdotなどをvxに変換していくために以下の数式を代入します。 pdot=-v+w×p から計算 subs(対象式, 代入対象, 代入する式)という構文を用いて代入。 diff(X(t), t)という式を消していきます。 syms vx vy vz wx wy wz before={diff(X(t), t),diff(Y(t), t),diff(Z(t), t)}; after = { -vx-wy*Z(t)+wz*Y(t), -vy-wz*X(t)+wx*Z(t), -vz-wx*Y(t)+wy*X(t)}; diff_p_XYZ = subs(diff_p,before,after) ...

誤植はどの本にもあることです。 なお，手元にあるのは第４版です。 カルマンフィルタの基礎 作者: 足立修一,丸田一郎 出版社/メーカー: 東京電機大学出版局 発売日: 2012/10/10 メディア: 単行本（ソフトカバー） 購入: 3人 クリック: 3回 この商品を含むブログ (5件) を見る 線形カルマンフィルタの終端値とリカッチ方程式 MATLABでリカッチ方程式を解く解法 誤植箇所 p149 演習６-２の解答 線形カルマンフィルタの終端値とリカッチ方程式 線形カルマンフィルタの事前推定の共分散をPとします。 十分時間が経って定常値となる時のカルマンゲインは この時，終端値において次の式が成り立ちます。 もっと長く書くとこう。 展開すると となり，これは代数リッカチ方程式（DARE）となります。 MATLABで解けます。 jp.mathworks.com MATLABでリカッチ方程式を解く解法 Referenceによると， だそうです。 従って， [P,L,G]=dare(A',C',Q,R) と実行することで解を得られます。 誤植箇所 p149 演習６-２の解答 解答ではPの最終値が1+√３としていますが， [P,L,G]=dare(1',1',4,1) を計算すると，2(1+√2)=4.8前後が正解の値になるはずです。 なお，スカラなので手計算でリカッチ方程式を解いても同様の値になります。 ドヤ顔で書いてますが間違ってたらすみません。

結論 numpy arrayの保存 こんなデータになる MATLABでcsvを開く csvread：’，’をdelimiterにしないとおかしくなる． readtable：有能 Reference 結論 書きはじめの時は，numpyで保存したcsvがmatlabに変な感じで認識されるのを不満に思っていたのですが，途中から自分の過ちに気づきました． 要点は以下の２つです． csvread()は区切りがカンマでなくては行けない そうでない場合はreadtable()を使う はい．帰ってよし． numpy arrayの保存 自分はもっぱらsavetxtを用います． numpy.savetxt("K.csv", K, delimiter =',',fmt="%0.14f") loadする時，そのままnumpy arrayとしてloadでき，他のソフトからも開きやすいからです．（Excelはダメ） こんなデータになる 糞長くなるので人間にわかりやすいデータではないですが，一応こんな感じにsaveされます． 嘘です．delimiter =','という引数を加えないと区切りがスペースになってしまいます． 其の結果が以下の通りです． 1.554141360341085672e+09 8.172370370370369983e+02 4.427456790123457040e+02 -1.855009145168301770e-03 -4.460917310385187178e-03 1.554141360408607960e+09 8.172024691358025166e+02 4.427777777777777715e+02 -1.670469061790032239e-03 -4.653048558273605202e-03 1.554141360473935604e+09 8.172383292383292428e+02 4.427518427518427302e+02 -1.816789199264026364e-03 -4.465982169869531920e-03 1.554141360534764767e+09 8.171901234567901611e+02 4.427753086419753004e+02 -1.689939994465252445e-...

Software Installation SeDuMi SDP3 YALMIP LMIの試験プログラム 式の解説 プログラム 配置結果 やってみた所感 Software Installation ソルバはYALMIP以外にCVXgen等もありますが，どっちがいいんですかね？議論ありますが有識者の意見を聞きたいです。 Yalmipまわりのインストールは以下のスライドが比較的参考になります。2018年9月の時点ですでに情報は古いですがURL等は参考になるでしょう。 http://www.maizuru-ct.ac.jp/control/kawata/iscie/doc/install_yalmip_new.pdf SeDuMi Downloads | SeDuMi 単にダウンロード MATLAB直下に解凍。 SDP3 SDPT3 Ver3-4をダウンロード。 MATLAB直下に解凍。 その後，解凍したディレクトリに移動して以下のコマンドを実行します。 Installmex(1) YALMIP 実はYALMIP単体でもそれなりにソルバはありますが，SDP（半正定値計画問題）を解くには上記のソルバが必要です。 yalmip.github.io addpath(genpath('C:\Users\<ユーザ名>\Documents\MATLAB\SeDuMi_1_3')) addpath(genpath('C:\Users\<ユーザ名>\Documents\MATLAB\SDpt3-4.0')) addpath(genpath('C:\Users\<ユーザ名>\Documents\MATLAB\YALMIP-master')) yalmiptest 実行結果は以下のとおりです。とりあえずSDPは解けているようなのでOKです。 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ | Test| Solution| Solver message| +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ | Core functionalities| N/A| Successfully solv...

GPSの緯度経度の行列から移動距離を計算するアルゴリズムがmatlab公式に用意されてないので自分で書きました。 探せばあるのかもしれません。 今回の手法ではヒュベニの公式というのを使います。 詳細な説明等はここのサイトがよいです。 画像も借りてきました。 関数は以下の通り。整備したらGithubにでも投げます。 function D = calcGPS_dist(LON,LAT) len = length(LON); D = 0; for i = 1:len-1 D = D+hubeny([LON(i),LAT(i)],[LON(i+1),LAT(i+1)]); end end function d = hubeny(L1,L2) % input: pair of [LON,LAT] % output : distance [m] % e2 = 0.00669437999019758; a = 6378137; b = 6356752.314140; e2 = (a^2 - b^2) / a^2; mu_y = (L1(2) + L2(2))/2; W = sqrt(1-e2*(sind(mu_y))^2); M = a*(1-e2)/(W^3); N = a/W; dx = (L1(1)-L2(1))/180*pi; dy = (L1(2)-L2(2))/180*pi; d = sqrt( (dy*M)^2+(dx*N*cosd(mu_y))^2 ); end ハマリポイントは角度をradに変換するのを忘れていたところです。

matlabコードをGitで管理していますがコマンドウィンドウからもgitを叩けるのをこのサイト見て知ったのでメモします。 外部コマンドの実行の仕方 一連のgit操作 注意点・エラー 自動でGitでpushする関数 外部コマンドの実行の仕方 じつは”!”を事前につけるとコマンドプロンプトのコマンドが実行できるようです。ipconfigなども打てました。 jp.mathworks.com 一連のgit操作 一連のGitの操作に!を前置すれば済みます !git add . !git commit -m "AutoMaticUpdate" !git push origin master なんとこれをmファイルでまとめて実行させることもできます。mファイルでシェル芸ができます。 これって脆弱性なんじゃないか… 注意点・エラー matlabのユニコードの設定のせいかなにやら文字化けが発生します。これは今の所困ってないので無視。 git commitでメッセージを指定しないと通常vim等のエディタが出ますが上記の文字化けと組まれると非常にきついです。メッセージ指定しようね。 gitのversionが古いとpush時に以下のようなメッセージが出ます。（文字化け含む） fatal: AggregateException encountered. 1 縺､莉･荳翫�繧ｨ繝ｩ繝ｼ縺檎匱逕溘＠縺ｾ縺励◆縲 bash: /dev/tty: No such device or address error: failed to execute prompt script (exit code 1) fatal: could not read Username for 'https://github.com': Invalid argument ここの記事にある通り公式から最新のverをインストールすることで解決します。 自動でGitでpushする関数 で，やっぱりこういうのは関数化したいわけで， GitComiiter関数を作成しました。evalを使っている邪悪なプログラムですが引数にCommit Messageを受けて動きます。 なお，デフォルト引数としてMessageがないときは自動的に適当なコメントをつけてくれます。 function GitComiiter(messa...

序論・参考文献 よくある疑問とエラー アッカーマン法 参考資料 Outline MATLABコード MIMOの場合：EigenvalueStructure法 序論・参考文献 制御工学の授業で一度は計算させられる極配置ですが MATLABにはplaceなる関数があって，実務を始めてからは全く手計算をしなくなります。 よくある疑問とエラー 授業では忘却の彼方にありますがB行列のrankより高い極配置はできないというエラーは誰しも見たことがあるかと思います。 以下のような現象です。 A = [1 2;3 4]; B = [1; -2]; pole = [-3,-3]; place(A,B,pole) Bのランクは1ですので，このコマンドを実行すると以下のエラーを起こします。 エラー: place (line 78) "place" コマンドは、rank(B) より大きい多重度を持つ極を配置できません。 そういえば中身どうなっているんだ？Cとかで実装するとどうなるんや？という疑問がありませんか。 MATLABのは論文読んでないのでまだ実装しませんが， 本稿で紹介するアルゴリズムを用いた場合はこの場合でも，ほぼ極配置ができます。 -3.0000 + 0.0000i -3.0000 - 0.0000i ひょっとしたらちょっと誤差があるのかもしれない。 アッカーマン法 SISO系限定ですがアッカーマン法であればこのような状況でも対応できます。 なお、MATLABではackerという関数名で実装されています。（昔はReferenceがあった気がするのですが。。。？） アッカーマン（Ackermann）はこの人？某ヒロインが検索汚染してきます。 アッカーマン法の詳細な記述についてはなんとwikipediaが一番スッキリしていて見やすくなってます。（英語） 原著の論文はこれっぽいです。 J. Ackermann. Der Entwurf linearer Regelungssysteme im Zustandsraum. Regelungstechnik und Prozessdatenverarbeitung, 7:297–300, 1972. 参考資料 どこかの授業スライドですがこれが具体例込でわかりやすいです。 また，参考になりそうなスライドがありました。（要DL，MI...

MATLAB LiveScriptとは 先人の調査所感 実際に書いてみよう ファイルの作成 よく使うショートカット 図や実行結果を載せる 数式を挿入 エクスポート 正規手段で Microsoft Print to PDF 総括 MATLAB LiveScriptとは ソースコードと説明用の式や図をInteractiveに表示するMATLABの新たな表示形式です。 MATLAB公式の説明。 ライブ スクリプトとは - MATLAB & Simulink - MathWorks 日本 Pythonをやる人にはjupyterみたいなもんといえばわかるでしょうか。 MATLABをjupyterで使うということもできますがPython環境を準備したりが面倒でしたよね。 なお，公式には2016aから対応しているとのこと。 Amazonからの購入。生協にあるなら生協の方が安いかと。） MATLAB and Simulink Student Suite R2018a 出版社/メーカー: MathWorks メディア: DVD-ROM この商品を含むブログを見る MATLAB2018のStudent版はDeep Learningを含む機械学習等にも力を入れているので以前のバージョンを持っている人も試してみると良いかもしれません。 先人の調査所感 なお，興味深い調査所感が以下にあります。実行速度が半分になる可能性があるのはたしかに厄介なのと，図を含めたファイルサイズが今の所の懸念点というのは有益な情報でした。 qiita.com 実際に書いてみよう Jupyterのような使い方をわかっている前提で書きますが，要はテキスト用の箇所とコードを記述する用の箇所を切り替えているだけなのではじめての人でもすぐ馴染めると思います。 ファイルの作成 ライブスクリプトの拡張子はです。MATLAB LIVE Scriptの略なんでしょう。 GUIメニューでライブスクリプトを選択すれば良いです。 コマンドラインで実行しても作成できるそうな。 よく使うショートカット 起動してぽちぽちいじっているとコードを書くところと説明を書くところがこのように違うことがわかります。 具体的に使いそうなショートカットは Alt+Enter：コード/テキストの切り替え Ctrl+Alt+Enter : ブロッ...

流石に適当に書きすぎたか ベイズの定理などはすでにMasterしているものとします。 基礎は確率ロボティクスで勉強しよう。 確率ロボティクス (プレミアムブックス版) 作者: Sebastian Thrun,Wolfram Burgard,Dieter Fox,上田隆一 出版社/メーカー: マイナビ出版 発売日: 2016/09/21 メディア: 単行本 この商品を含むブログ (1件) を見る ちなみに英語版のドラフトは以下のURLに落ちています。 https://docs.ufpr.br/~danielsantos/ProbabilisticRobotics.pdf パーティクルフィルタとは matlab実装 状態推定 観測に依る更新 個人的な解釈 リサンプリング 尤度の設定は重み付き平均でいいのか？ コード コード置き場 実装のコツ 結果 総評・所感 パーティクルフィルタとは ベイズフィルタの亜種。状態推定に必要な確率分布を数式ではなく，多数の粒子群の集合として表現する手法です。 状態方程式とノイズの分布がわかっていればどのような関数でも同じように実装できるのが強み。 wikiを見るか以下の説明が簡潔かと。 『カルマンフィルタもパーティクルフィルタも全部ベイズフィルタで説明がつく』という説明が美しいのでやはり上記の確率ロボティクスの本を読むことをおすすめします。 パーティクルフィルタ - おべんきょうwiki - アットウィキ 実装などについては以下のサイトが参考にしやすい。 Tutorial/Practice0 – MIST Project Pythonの実装例でよければ確率ロボティクスを和訳した上田先生が公開しているモンテカルロ自己位置同定の話なんかが実装の参考になることでしょう。 github.com 後は鏡先生のスライドとかも参考になるか？ http://www.ic.is.tohoku.ac.jp/~swk/lecture/ic2012/kagami_ic20120710.pdf matlab実装 というかMATLABの実装もAtsushiさんがやってくれてます。自己位置同定に関してはこれを使えば良いでしょう。 Particle Filterを使用した自己位置推定MATLABサンプルプログラム - MyEnigma もう少し汎用的な話な...

MATLABのplot関数で作成した図の色にいつも困っているので関数作成。 テスト用データ作成 Defaultカラー抽出 どんな色がいいか？ 色変更 自動色つけ関数 関数本体 色盲用colormap テスト用データ作成 これをプロットしていきます。 t = 0:0.01:3; y = [sin(t);sin(t+0.1);sin(t+0.2);sin(t+0.3);sin(t+0.4);sin(t+0.5);sin(t+0.6);sin(t+0.7)]; Defaultカラー抽出 MATLABのデフォルトカラーはそこそこいい感じです。 RGB values for 2014b default colors - MATLAB Answers - MATLAB Central get(groot,'DefaultAxesColorOrder') 先程のデータをプロットしましょう。 0 0.4470 0.7410 % 青 0.8500 0.3250 0.0980 % 赤 0.9290 0.6940 0.1250 % 山吹色 0.4940 0.1840 0.5560 % 紫 0.4660 0.6740 0.1880 % 緑 0.3010 0.7450 0.9330 % シアン 0.6350 0.0780 0.1840 % 茶色 どんな色がいいか？ さて，そのうちどのような色を用いればいいかですが，matlabの設定はまぁ悪くないかなぁと思っています。 あと，使えそうな色はオレンジoでしょうか。 どんな色がいいかとかは論文が幾つかありますがこんなの読んでいる場合じゃねぇ。 http://www.perceptualedge.com/articles/visual_business_intelligence/rules_for_using_color.pdf 色変更 plotする時にRGBで色を指定するには以下のように書けばよいですが，正直言ってこんなの絶対書きたくない。 plot(t,y,'Color',[1., 0.5, 0.]) 一方で，Figureハンドルを入手していれば次のように代入形式で後付けで色を変更できます。 fig=plot(t,y) fig(5).Color = [0,0.7,0] fig(6).Color = [1., 0...

カルマンフィルタから信頼区間の線を引きたいというのは当然の要求だと思うのですが見渡す限り 良さげな文献が無いので自分で書きます。 書いておいて何ですが自信がないので訂正等のご意見を随時求めております。 あと，ブログの本文もMarkdownにシフトしてみました。 [tex: ]ってうつのめんどくさいな。 理論 参考文献 カルマンフィルタにおける信頼区間とは n次の状態変数における信頼区間 別解：シグマ区間からの拡張 Matlabによる実装１ カイ二乗分布の境界値 Symbolic Math Toolboxを用いた代数方程式の解法 MATLAB関数として実装 追記＆最終版：実は超簡単 疑問点 理論 参考文献 きちんとした参考文献を見つけられずに困ってます。 一番まともそうなのがStackOverflowしかないんだけど… python - Estimating confidence intervals around Kalman filter - Stack Overflow あと，確率論とかについては以下の本を。 確率ロボティクス (プレミアムブックス版) 作者: Sebastian Thrun,Wolfram Burgard,Dieter Fox,上田隆一 出版社/メーカー: マイナビ出版 発売日: 2016/09/21 メディア: 単行本 この商品を含むブログ (1件) を見る 足立先生の本も良著なのでオススメ。 カルマンフィルタの基礎 作者: 足立修一,丸田一郎 出版社/メーカー: 東京電機大学出版局 発売日: 2012/10/10 メディア: 単行本（ソフトカバー） 購入: 3人 クリック: 3回 この商品を含むブログ (4件) を見る 追記：それっぽい論文見つけました1966年の論文です。 Confidence, Prediction, and Tolerance Regions for the Multivariate Normal Distribution on JSTOR カルマンフィルタにおける信頼区間とは まずは信頼区間についてwikiで確認してください。 信頼区間 - Wikipedia 以下の図のように分布の形状と平均分散などの値がわかっていれば値がどの確率でどの範囲に収まっているかがわかります。 カルマンフィルタでは誤差...

jupyterの使い勝手がだんだんわかってきたので本格的に使うにはどうすれば良いのか模索中です。 想定した運用 ショートカットキー LatexやPDFへの変換 GUIでやる コマンドプロンプトから pandocのインストール 文献 Latexマクロをjupyterで組む。 その他のExtensionツール 想定した運用 コードの実行できる文書として使用したい。 従って，matlabコマンドと図，Texの数式の機能などが過不足なく使えるように。 Latexはマクロを使えるようにしたい。 簡単な文書として運用したい notebookで書いたものをlatexやpdfに変換すればそのまま簡易資料にできるようにしたい。 ショートカットキー マウスをなるべく使いたくないのでショートカットキーを覚えるべきですね。 Jupyterのショートカット - Qiita 実行の仕方から， 実行：Ctrl+Enter 実行&改セル：Shift+Enter マークダウンとコードの切り替え マークダウンへと移行：Esc+m コードへと移行：Esc+y そしてセルのコピペ コピー：Esc+c ペースト：Esc+p 削除：Esc+dd LatexやPDFへの変換 GUIでやる [File]から[Download as]で好きなフォーマットを選択できます。 コマンドプロンプトから jupyter nbconvert --to latex <HOGE>.ipynb と実行することで同様のことができます。 pandocのインストール そのままだとpandocがありませんという旨のエラーが出ます。 anaconda search -t conda pandoc の出力を見た感じ，conda-forgeが良さそうですので anaconda show conda-forge/pandoc とうってURLを確認します。 conda install --channel https://conda.anaconda.org/conda-forge pandoc これで，Latex形式で落とせます。 詳しい説明は別記事で。 ossyaritoori.hatenablog.com 文献 jupyter notebookをLaTeXに変換 - Qiita Jupyter nbconvert（ファイル変換）メモ - はしく...

ArduinoとMATLABを連携すると要はコンパイルなしにいろいろデバッグができて作業効率が良いということです。 Pythonでもできそうなもんだけど誰か代わりに調べてください。 この30分のビデオ見ればここにあることは全部わかります。 Using Arduino with MATLAB and Simulink - Video - MATLAB パッケージのインストール Arduinoとの接続と試運転 接続 サーボモータの制御 他のmethodの確認 連続的に制御するには？ Simukinkでもできます DCモータの制御 特定の部品を使っている場合は簡単 パッケージのインストール matlab2014以降から対応です。 Support Package InstallerからArduinoに関する追加のパッケージをインストールします。 matlabコンソールから試しに a = arduino(); とか入力してみて出たエラーからも簡単に飛ぶことが出来ます。 あとはひたすらポチポチするだけです。大体15分弱かかります。 Arduinoとの接続と試運転 Control Servo Motors - MATLAB & Simulink Example 接続 ホーム画面から「コンピュータの管理」というのを探し当ててそこからデバイスマネージャーのポートという所を探してください。 ここから得た情報を元に次のように入力してみます。 a = arduino('com4', 'Mega2560', 'Libraries', 'Servo'); サーボモータの制御 サーボモータのピンを調べます。私のはモータシールドの「OUT5」につながっているので以下のドキュメントからD5とわかります。 Arduino Motor Shield Rev3 s = servo(a, 'D5') s = Servo のプロパティ: Pins: D5 MinPulseDuration: 5.44e-04 (s) MaxPulseDuration: 2.40e-03 (s) 実際に遊んだ時のコードです。writePosition(s, angle);とcurrent_pos = readPosition(s);のようにして %% write value Base_angle = 90;...

ポ○ットモンスター 金/銀 みたいなタイトルですが、iPythonのセットアップについてです。 「なみのり」を覚えてから来てください。 前提：Python環境 Anacondaを使う場合（簡単・ストレージを割りと食う） Minicondaを使う（ストレージを沢山消費したくない人向け） 注意書き jupyter notebookをmatlabで使う 0．Python（仮想) 環境構築 １．Matlab engineのインストール エラーが出たら？ ２．必要パッケージのインストール 動作確認 jupyter notebook を GoogleDrive上で扱う。 インストール 起動 Matlabが古いけど無理やりでもこの機能を試したい！ matlab live editorという逃げの1手 さぁ、闇のゲームだ！ 第二ラウンドだ！ 前提：Python環境 私の環境ではとりあえずWindowsでやるつもりでいます。 あと、GitHubの一部コマンドが必要なので適当に調べて入れておきましょう。 私家版 Git For Windowsのインストール手順 | OPC Diary Anacondaを使う場合（簡単・ストレージを割りと食う） Anacondaは一連のPython関連の便利パッケージを落としてくれるのでPythonに今後お世話になる気のある人はこちらを使いましょう。 下の記事を参照。 ossyaritoori.hatenablog.com 画像関係に興味がある方はこちらの記事も踏んでおいたほうがいいですが、まぁ後からできるので無理にすることはないです。 ossyaritoori.hatenablog.com Minicondaを使う（ストレージを沢山消費したくない人向け） 先ほどのAnacondaの機能絞った板がこちらです。 Miniconda — Conda Pythonとcondaのパッケージ管理だけなので数百Mの容量で済むはず（試してない） どのみちやることは同じです。 注意書き jupyter in matlabは 今のところmatlab2014b以降,対応するPythonのバージョンもmatlabのバージョンに依存するようです。 jupyter notebookをmatlabで使う 以下の参考サイトに準拠してすすめます。 Matlab-based IPyth...

今回はfind関数を主に使っていきたいと思います。 find：特定の条件に一致する部分を調べる 等号条件 不等号条件 findを使わない場合の記法との比較 特定の条件に一致する部分を置換 特定の条件に一致する部分を削除 関連する文献 find：特定の条件に一致する部分を調べる 等号条件 例＞ 行列 A から値が0である箇所を探す。 find(A==0) この際の出力は１次元の配列になります。m行n列のi行j列目にある場合は，i*n+j として出力されることに注意。 以下のようなコードを書いて位置との対応を掴んでください。 mg = magic(4) find(mg==1) なお，後ろにfind(A=＝0,3)などと数値を列挙すると条件に一致する部分のうち数値の数だけ列挙してくれます。 不等号条件 例＞ 行列 A から値が10より大きい箇所を探す。 find(A>10) findを使わない場合の記法との比較 findを使わずとも A==0 という表記単体で，条件に一致するか否かの論理値を得ることができます。 A = magic(3) A = 8 1 6 3 5 7 4 9 2 >> A>5 ans = 1 0 1 0 0 1 0 1 0 findではこの論理値行列の代わりに条件を満たす行列のindexをそのまま返します。 >> find(A>5) ans = 1 6 7 8 特定の条件に一致する部分を置換 先程調べた箇所を行列に代入することで置換できます。 試しに行列の中で特定の条件を満たす部分を0にしてみましょう。 mg = magic(4) place = find(mg>10) mg(place)=0 魔法陣mgのうち，10より大きい成分の位置がplaceに入り， これは，次のように縮めて書くことが多いです。 mg = magic(4) mg(find(mg>10)) = 0 追記：先程述べた表記を活用すればfindを省略しても以下のように書くことも可能です。 mg = magic(4) mg(mg>10) = 0 特定の条件に一致する部分を削除 同様に特定の値ではなく [] を代入することで該当する箇所を削除可能です。 例＞ベクトルAのうち，0に該当する箇所を削除。 A(find(A==0)) = [...

背景 古来より我々はmatlabの作るPDFのへんてこなことに悩まされ続けてきました。 ossyaritoori.hatenablog.com ossyaritoori.hatenablog.com 他の人のツールを使おう 以下のところからツールを落としてきます。 github.com 使い方 : テンプレ このツールには主にpubfigとexpfigの2つがあります。 以下その例をば， %figure定義 hfig=figure(1); plot(～～);%なんか入れてplotしてね grid on; xlabel('frame number'); ylabel('error') legend('A','B','Location','best') title('Pyoyon') pfig = pubfig(hfig); pfig.LegendLoc = 'best'; pfig.FigDim = [15 11]; expfig('filename','-pdf'); loationをいじるところとか被ってたりするけどとりあえず何も考えなくともキレイな図が出るのが良いですね。 2018年4月追記）最後から二行目のところですが，FigDimという要素が削除されたようなので Figの大きさをいじる。 ソースコードはとても難解でしたがFigDimというクラスでFigureのサイズを司っています。 なおMATLAB公式によるとfigのうちのpositionというハンドルをいじれば変更できます。 https://jp.mathworks.com/help/matlab/ref/figure-properties.html とりあえずコンパクトなサイズとして以下を使うようにしています， pfig.FigDim = [15 11];

過去記事が地味に反応あったので調子乗って続きをば。 ossyaritoori.hatenablog.com 前回のGauss-Newton法に引き続き， Levenberg-Marquardt法についてです。 フィッティングの手法としてはGauss-Newtonより多用される王道中の王道ですね。 概要 アルゴリズムについて私的見解 知恵袋の秀逸な解答 問題設定 Gauss-Newton法との比較 補足：パラメータ空間でのアルゴリズムの様子 Levenberg-Marquardt法：コード マッピングのコード 概要 前回のGauss-Newton法を前提にしているのでそっちを先に参照してくださいな。 matlabはstudentバージョンで2015aで動作します。 アルゴリズムについて私的見解 以下の記述は私の持ってるイメージですのであしからず。 Gauss-Newton法は2次収束するNewton法において，ヘッセ行列をヤコビ行列の二乗で近似した手法ですが， ダンピングファクターλを下のように導入します。 適度にダンピングをかけることで吹っ飛ぶのを回避できるってわけ。 なお，ダンピング項のIは最も簡単には単位行列，もう少し賢く行くにはヤコビアンの二乗の対角項を使うようです。 matlabで書くならダンピング項はこんな感じか，簡単すな。 I = diag(diag(J.'*J)) 知恵袋の秀逸な解答 なお，以下の知恵袋の解答も秀逸です。 Levenberg-Marquardt法（LM法）は非線形最小二乗問題を解く有力な... - Yahoo!知恵袋 「最急降下法とGaussNewton法の組み合わせ」という言はたしかにしっくり来ます。 以下の式は数学的には超間違い!!ですが， この式は怖い人がなんか言ってきたら消します（） 問題設定 前回のと同じ，wikipediaの例題です。 ガウス・ニュートン法 - Wikipedia 関数： Gauss-Newton法との比較 λをいじった場合どのように収束速度が変わるかを示す図が以下の通り。 Legendが適当なのはごめんなさい。ダンピングと言う通り，λを大きくすると収束速度は鈍ります。 補足：パラメータ空間でのアルゴリズムの様子 変更するパラメータa,bを０～１の範囲で変えて残差の二乗和の関数の...

多数のデータ群から関数フィッティングを行う場合に非線形最小二乗法を用います。 MATLABにはOptimization Toolboxなどがあり，それを用いると簡単にフィッティングできます。 MATLABで非線形最小二乗問題を解く - Py このようなToolboxなしでも， アルゴリズムの基礎がわかっている前提で書きます。 ダンピングを入れたVerはこちら。 ossyaritoori.hatenablog.com アウトライン 問題設定 問題の定式化 ヤコビ行列の作成 補足：ヤコビ行列のサイズ 注意：subsを用いた複数の配列の代入 反復計算による求解 補足：行列のサイズ 計算結果 まとめと今後 サンプルコード 参考資料 jacobian 関数の使用 ：2017/3/7 更新 アウトライン matlabのバージョンは生協で買ったStudent version（2015a準拠)です。 偏微分」と「変数代入」を行います。 ヤコビアンを求めて，Gauss-Newton法の枠組みでこれを計算します。 こんな風にフィッティング出来ます。 なお，Gauss-Newton法は案外簡単に発散するのでその場合は以下のLevenberg-Marquardt法を試してみてください。 ossyaritoori.hatenablog.com 問題設定 例題としてwikipedeiaの問題を考えます。 ガウス・ニュートン法 - Wikipedia Fitting対象の関数は以下のようになります。 フィッティングするのはaとbの2つの変数です。 フィッティングするデータは次の様に表せます。 % data for fitting xhat = [0.038, 0.194, 0.425, 0.626, 1.253, 2.500, 3.740].'; yhat = [0.050, 0.127, 0.094, 0.2122, 0.2729, 0.2665, 0.3317].'; 通常こういう値のベクトルは縦向きですので後ろで転置してます。 今回は7×1のベクトル，すなわち7個のデータセットですね。 問題の定式化 残差と最小二乗和を以下の様に定義します。iはデータ群の各成分ってこと。 残差： 残差の二乗和： Symbolic math Toolboxを使ってこれは次の様にかけます。 ...

ダジャレみたいな軽い記事です。 NaNを0に置換 Infを0に置換 種明かし よくわからんので実験 NaNを0に置換 Aという行列の中にNanがあってそれを0にする場合，以下のように書きます。 A(isnan(A))=0 なお，matlabではfor文はべらぼうに時間がかかるので Infを0に置換 Aという行列の中にInfがあってそれを0にする場合，以下のように書きます。 A(~isfinite(A))=0 種明かし isnan(A) というコマンドは行列Aのうち，Nanがある所を1，無い所を0の「論理値」で出力します。 よくわからんので実験 以下のコードで実験してみました。 3*3の魔法陣を作ります。 Magi =[ 8 1 6; 3 5 7; 4 9 2] 代入するのは代数値でなければなりません。 Magi(diag([true true true])) 結果，ans = 8,5,2 と見事に対角の値だけにアクセスしています。 ならば，一つだけ食わせるとどうか？ Magi(true) 結果：ans = 8 そして数値を足してずらせるようです。ポインタみたいだな。 Magi(true+1) 結果：ans = 3 行 → 列の順に格納しているのがわかります。 Magi(true+3) 結果：ans = 1 ほらね。 また，こんなことも出来ます。 Magi(true,true+1) 結果：ans = 1 う～ん，true=1として処理してるだけに見えてきた。 Infの方はなんか出力が逆ですが基本やってることは同じです。 ほ～ん。

完成図 と 目的 使う関数area()について 塗りつぶし用コードまとめ 関数化（注意点あり） 使い方 困った点・改善点 2017/2/12 ：追記 補題：spline補間とPCHIP補間 完成図 と 目的 研究柄matlabで図を書くことが多いですが，稀にmatlabの背景領域に色を塗りたくなります。 要はこういうグラフが書きたいってわけ。 使う関数area()について 今回はarea()という関数を使います。 Referenceは以下のとおりです。 2 次元プロットの塗りつぶし領域 - MATLAB area - MathWorks 日本 今回は area([X1 X2],[Y1 Y2],basevalue); という形で使っていきます。 塗りつぶし用コードまとめ 他にも枠線を消したり，グリッドを前面に出したり，ゴニョゴニョした結果次のような感じになります。 gcaの'Children'の順番をいじることにより変更可能ですね） [X1 X2],[Y1 Y2],baseline); % 塗る範囲([X1 X2],[Y1 Y2])とbaselineの間 set(ha,'FaceColor',[1.0 1.0 0.9]); % 塗りつぶし色指定：今回はクリーム色 set(ha,'LineStyle','none'); % 塗りつぶし部分に枠を描かない set(gca,'layer','top'); % gridを塗りつぶしの前面に出す hAnnotation = get(ha,'Annotation'); hLegendEntry = get(hAnnotation','LegendInformation'); set(hLegendEntry,'IconDisplayStyle','off'); % legendに入れないようにする 従ってこの後に hold on plot(...) grid on grid off とやることによって完成です。 しかしこのままだとベースラインが見えてしまったり，プロット範囲を手動で変更したりが結構面倒です。 関数化（注意点あり） 同期I氏によってもっと便利に使えるよう関数化を施されました。 自動でプロット範囲を判別し，余計な所を塗らない点 Childrenのハンドラをいじることで表示の優先順を変えられるため，後付けでプロ...

突然ですがmaltabのreferenceページって結構見づらくないですか？ 特にClassとして実装されてる奴らはどういう処理をするのかかなりわかりづらいと私は思っています。 fft in MATLAB DFT用のコード（FFTも同じ） fft in MATLAB fftコマンドは非常に便利ですが， ナイキスト周波数で折り返しが生じること DFT用のコード（FFTも同じ） 以下のように関数化しました。パワースペクトル密度の表示まで意外とやることがあるのでこれで結構快適になったと思われます。 11/08 追記：直流成分の除去も中でするようにしました。単に平均を引いてやっています。 function [f power] = DFTandShow(t,x) if size(t) == 1 sampling = t; else sampling = t(2) - t(1); end % transpose so that row become longer if size(x,1) < size(x,2) x = x.'; end sf = 1/sampling % Sample frequency (Hz) L = max(size(x)) % data length(window length) % calc average average = ones(1,L)*x(:,1) / L; x = x - average; % fft N = pow2(nextpow2(L)); % Transform length N = 512; y = fft(x,N); %fft f = (0:N-1)*(sf/N); % Frequency range ( sampling fq / fft length ) power =abs(y/N) *2 ; % Power of the DFT figure; plot(f,power) title('Single-Sided Amplitude Spectrum of X(t)') xlabel('frequency (Hz)') ylabel('Power') grid on; xlim([0 f(floor(N/2))]) % cut off faster frequency than nyquist f...

随分昔に遊んだコードを掘り出して懐かしくなったのでネットに放流します。 SSD)とはまた違う画像の相互情報量なるものを計算して画像を比較する手法の紹介です。 今現在この手のコードは以下のリポジトリにて管理しています。 EvaluateImage/evaluate_image.m at master · YoshiRi/EvaluateImage · GitHub Gist版はZNCCのコードに誤りがあります。（あとで直します） 相互情報量とは 画像の相互情報量 実際の比較：地図画像と航空写真のマッチング 入力画像 結果の紹介 比較手法 実際のコード コードを組む際のコツ：軽量化 画像間の相互情報量を計算するコード 比較手法のコード プロット用コード 参考論文 相互情報量とは ある確率pで起こる事象の情報量Hは次のような式で表されます。 ある事象Aのもつ情報量をH(A)，BのをH(B)とすると，それらの間の相互情報量MI(A,B)は次の式で計算されます。 ここで，H(A,B)はAもBも起こる時の情報量です。 画像の相互情報量 では，この話を画像に持ち込むとどうなるかといいますと， 事象：画像の中である画素値が出現する確率 p_I(i)とする。すなわち，0~255の画素値のうちxという画素がどの頻度で出現するかで先程の情報量を定義します。 画像Iの持つ情報量はIの中の全ての画素値の情報量の総和 として表され，これは実際の計算では出現回数を確率として正規化したヒストグラムを用意することになります。 同様にして結合情報は， ヒストグラムを用いて同様に情報量を計算します。 注意として，例えば輝度を256階調にするとこのヒストグラムは256*256の結構大きな反復試行を要求されるので計算量が大きくなります。 最終的に画像同士の相互情報量は以下の式で計算されます。 この値が大きいほど，2つの画像の相関が高い，すなわちより似ていると評価できます。 実際の比較：地図画像と航空写真のマッチング 地図画像と航空写真をあわせるのは人間にはなんてことない作業ですが， 相互情報量はそのあたりなんとか出来るらしい[参考文献1]です。試してやんよ。 入力画像 航空写真 map_MI1 GoogleMapの地図画像 map_MI2c 結果の紹介 地図の画像を順次平行移動させな...

この記事はメモです。 https://jp.mathworks.com/matlabcentral/fileexchange/48744-lucas-kanade-tutorial-example-1/content/LucasKanadeExample1/html/LKExample1.html 関数名：OpticalFlow オプティカルフローを図としてアウトプット 参照元からの変更分 関数への変更 Conv()からfilter2()関数の使用への変更 コード %% Optical Flow % This is not my original code % https://jp.mathworks.com/matlabcentral/fileexchange/48744-lucas-kanade-tutorial-example-1/content/LucasKanadeExample1/html/LKExample1.html % reviced conv() use to filter2 function OpticalFlow(img1,img2) im1t = im2double(img1); im1 = imresize(im1t, 0.5); im2t = im2double(img2); im2 = imresize(im2t, 0.5); ww = 45; w = round(ww/2); % Lucas Kanade Here % for each point, calculate I_x, I_y, I_t fx=[-1,0,1;-1,0,1;-1,0,1]; % フィルタ行列の作成 fy = [-1,-1,-1;0,0,0;1,1,1]; Ix_m=filter2(fx,im1,'same'); % フィルタリングX方向微分 Iy_m=filter2(fy,im1,'same'); % フィルタリングY方向微分 It_m = im2 - im1; % Ix_m = conv2(im1,[-1 1; -1 1], 'valid'); % partial on x % Iy_m = conv2(im1, [-1 -1; 1 1], 'valid'); % partial on y % It_m = conv2(im1, ones(2)...

聞かれたので昼休みに作ったまでです。 MATLABでたくさんFigureを出した時に比較する手段として動画にするというのはよくあると思われます。 本プログラムのステップは 描画 getframeによる取得 Videowriterによる書き込み の３ステップになっています。 コード どうぞご自由に。コードでは2次系の粘性項を変えた場合のステップ応答を作成していますね。 %% make movie from MATLAB plot clear all; close all; Frate = 5; k = 0:0.1:3; for i = 1:size(k,2) % First, make figure sys = tf([1],[1 k(i) 1]); figure(1); step(sys); xlim([0 10]); ylim([0 2]); grid on; % Put current figure into Frame Frame(i) = getframe(1); end % % convert the Frame to movie and show % figure(2); % movie(Frame,1); % write to video v = VideoWriter('step.avi'); v.FrameRate = Frate; % Framerate open(v); writeVideo(v,Frame); close(v);

前回のコードとの違い:簡単な場合に合わせた 前回のコードよりシンプルに，Ⅰ入力Ⅰ出力系の係数のみを求めたいってなった場合，若干式が変わります。 前回の記事 ossyaritoori.hatenablog.com 目的式 なんてことはありません。1入力1出力系の傾きを求めたい時に使ってください。 もっとざっくり言うとOutputがInputに比例する場合ですね。 この式を最小化します。 計算式 ほっとんど変わらないです。 うん。 コード（MATLAB関数）～ 7/8追記 ～ とことん簡単にしました。関数を呼んでその度出力結果Ynと入力Znをわたすだけです。 MATLABにはCでいうところのstaticの代わりにpersistentという型が用意されているのでそれを用いてできるだけ これを使うのが面倒とかいうあなたはもうRLS諦めたほうがいいんじゃないかってくらいです。 コード（クラス） 入力の次元のみを使うので次元が半分になった以外なにも変わらないです。 つかいかた 同じフォルダにおいて以下のように初期化とステップごとにパラメータ更新を行えばいいです。 % 初期化 est = rls_single(1,0.98); % パラメータ更新 Estimated_val = est.estimate(Output,Input); コード（C言語） 6/8 追記 サンプルデータにあからさまにひどいのが居た場合， // Forgetting factor double lambda = 0.99; static double Pz = 1000; //入出力はまぁ適当に与えてやってください // Prepare for caluclate double Num = lambda + Zn * Pz * Zn; double Ln = Pz * Zn / Num; double En = Yn - Zn * Zest; if (Zest + Ln * En < 0 || Zest + Ln * En > 1.0 / 1000){//Updateを無視する条件の例 //No update } else{ // Update Pn Pz = 1.0 / lambda * (Pz - (Pz * Zn * Zn * Pz / Num)); // Upda...

2016/7/4 追記 背景：MATLABで使いやすい逐次最小二乗法のコードがない 組みました ざっくりした説明？ 具体的な式 その他、忘却係数について ソースコード （関数版） ソースコード （クラスファイル版） 2016/7/4 追記 1入力1出力の系での推定は次の記事プログラムを参考にしてください。 ossyaritoori.hatenablog.com 背景：MATLABで使いやすい逐次最小二乗法のコードがない 僕が無知なだけの可能性は高いですが， MATLABで逐次最小二乗法をしたいなってなった時にいい感じのコードがないんですよね。 組みました はい。今回は関数ファイルとクラスファイルの2通りで実装しました。 MATLABでオブジェクト指向の書き方をしたのは初めてで友人にいろいろ教わりながらやりました。 逐次最小二乗法は結構書き方にいろいろ流派があるように思いますが， http://yoh.ise.ibaraki.ac.jp/edu/graduateschool/Project1.pdf ある程度原理を理解すれば関数自体を組むのは簡単ですが ざっくりした説明？ ここでは次のような線形モデルを仮定します。 逐次最小二乗法の目的は，ざっくり言うと 具体的な式 推定式の導出はやや面倒で結果を示すと以下のような式で更新をしていきます。 実質，現在の出力：Yn と その他、忘却係数について 初期値は今回結構雑でもいいのですが忘却係数の扱いを少し特殊にしました。 本プログラムでは最初は忘却係数を小さめにすぐに結果をアップデートするようにし、 ふぅ。文章打つほうが大変ですね。 追記＞ ソースコード （関数版） ソースコードは今回Gistというサービスを使ってみました。 はてなブログでソースコードを貼り付ける方法を調べてみた - ウェブと食べ物と趣味のこと 毎度毎度長いのを貼ってしまうのでこれからは長い奴はこういう形式にしようかと 逐次的に保存すべき変数：共分散行列Pnとθn 最初の行を削除すれば忘却係数Rhoは消してもいいかと思います。 ソースコード （クラスファイル版） クラスで定義したソースコードがこちら。 上が定義したクラスファイルで 珍しくいろいろコメント入れましたが，やっぱりあんまり読みやすくはないですね。 コメントや修正などありましたら...

シンボリックでいろいろ作ったコードを保存するコーナーです。 ＊ コードの内容 積分器の係数を決定するところまで 特性方程式det|sI-(A-BF)|の極を自由に配置するFを自動で決定する(係数比較 Fの計算方法は僕の知る限り 中身としては solveにおいて作成するeqnsはベクトルor行列形式でもいい detの特性方程式を計算するcharpoly関数を確認した って所の2点を新たに勉強したということになります。 これがsimulinkのブロック線図。 重根極配置ではオーバーシュートが避けられず それを吟味するのは今回の趣旨じゃないのでコードがキモイのも含めてご勘弁。 %% プラント A = [0 1 0; 0 0 1; 9 7 -1] B = [0 ; 0 ; 1] C = [ 1 1 0] %% 1型サーボ系にするためのフィードバック設計 Abar = [A [0; 0; 0];-C [0]] Bbar = [B ;0] % syms f1; % syms f2; % syms f3; % syms f4; %Fbar = [f1 f2 f3 f4] Fbar = sym('f',[1 4]) abf = sym(Abar)-sym(Bbar)*Fbar % snum = simplify(eig(abf)) Seq = charpoly(abf); %% pole placement freq = 0.5; omega = 2 * pi * freq; % size(Seq) pole = [1,4*omega, 6 *omega^2,4*omega^3, omega^4]; %% solve equation eqns = pole == Seq; Ans = solve(eqns); %% outputs K1 = double(- Ans.f4) f = double([Ans.f1,Ans.f2,Ans.f3]) %% 同一次元オブザーバー Kmat = sym('K',[3; 1]) Seq2 = charpoly(A - Kmat*C); %% pole placement2 freq2 = 25; omega2 = 2 * pi * freq2; syms s; pole2 = fliplr(double(...

symbolicmath Toolboxを最近使うようになったもののいろいろとよくわからないことが多いので 環境は MATLAB2015a です。 solveで方程式を解いたのにRootOf(~,z,1)となってよくわからない symbolicで方程式を解かせるという思考停止をしていると次数の高い場合などには 【 コード入力例 】 ^4 + x^3 + a == 0, x) 【 実行結果 】 ans = RootOf(z^4 + z^3 + a, z) やる気あんのかｗ 解決策 実はsolveには変数’MaxDegree’というものが存在しそれ以上の次数の式は解こうとしません。 この場合は’MaxDegree’を4に設定することできちんと解くことができます。 【 コード入力例 】 ^4 + x^3 + a == 0, x, 'MaxDegree',4) pretty(S) %図示 どのみちかなり長い式が出ます。prettyというコマンドを打つことである程度まとまった形で見せてくれますがそれでも なっが！！！！

MATLABのシンボリック変数を用いた計算はとっても便利ですよね。 コーディングも簡単なので効率良く使えるなら積極的に用いていきたいですが， シンボリック変数の行列やベクトルへの代入 以下のように事前に用意した行列にシンボリック変数を代入したい時があります。 mat = zeros(2); syms x y mat(1,1) = x^2+y; ちなみにこのコードは以下の様なエラーを吐きます。 sym から double への変換中に、以下のエラーが発生しました: シンボリック変数は普通に定義した行列などには代入できないようです。 解決策 10分ほどハマって調べた後，行列をシンボリックに変換すれば良いことを知りました。 変換したい行列 M を sym(M) のようにしてシンボリックの行列に変換することで対応できます。 改善したコードは以下のとおりです。 mat = sym(zeros(2)); syms x y mat(1,1) = x^2+y; 参考： シンボリック変数、シンボリック式、シンボリック関数、シンボリック行列を作成 - MATLAB sym - MathWorks 日本