MATLABマニュアル（超整理版） by 小笠原 2016/3/28
目 次
１．はじめに
２．基本
３．グラフィクス具体例
3.1.マッピング
3.2.時系列
3.3.図のレイアウト・デコレーション
3.4.色
3.5.視点を変える・光をあてる
3.6.一部をプロットしたくないとき（欠測値）
4.1.時系列データ解析（例：地震波形）
4.2.インバージョン（例：最小自乗法）
5.0.超簡単アスキーファイル読み込み法
5.1.一般的な読み込み方
5.2.バイナリファイルの読み込み方
5.3.キーボードやマウスからの入力
９．無くて七癖
MATLABを使うと、エクセルよりもかなり速く仕事ができる。
例：MATLABのt=0:0.01:1;plot(t,sin(2*pi*t));をエクセルでやると。。。
MATLABを使うと、エクセルよりもかなり高級なことができる。
例：MATLABのt=0:0.01:1;y=sin(2*pi*t)；fft(y);plot(abs(y))
MATLABを使うと、一般の言語よりもかなり楽に科学計算プログラムが組める場合がある。
例：2X2の２つの行列a,bのかけ算をMATLABで書くとa*b。言語で組むと。。。
困った時は、 ヘルプデスク （かゆいところに手が届かないときもあるが）
初めての人で、余裕がある人は、helpdeskのGetting Startedを通読することをお薦めする。
便利な違い
配列の型もサイズも宣言不要
プログラムを一行ずつ打ち込んでその場で実行させられる
一般言語だと一々細々 ＜---＞ MATLABは「これやって」と書く程度
複数行を書いてファイル名.mで保存→ファイル名＋改行→実行
グラフの読みとりや、計算パラメータのＧＵＩ調整も楽。
注意すべき違い
基本的に、変数は行列扱い、演算も行列演算。
* と .* の違い。/ と ./ あるいは ＾ と .^ との違い。
; の意味
a(1,:)の意味
t=0:0.01:1とすると
その他の演算子に関する詳細は、下記を参照のこと
file:////Ch1tp033_mino/研究室共通/MATLAB/jhelp/techdoc/refbycat/group2.html
３．グラフィクス具体例
MATLABのヘルプ → 例題とデモ → ビジュアライゼーション。
3.1.マッピング関係
例１：等間隔全格子点上に値→等高線 contour(X,Y,Z)
：等間隔全格子点上に値→３次元立体表示 mesh(X,Y,Z) または surface(X,Y,Z)
：等間隔全格子点上に値→ベクトルを矢印表示 quiver(X,Y,U,V) または quiver3(X,Y,Z,U,V,W)
demo_contour
MATLABのヘルプ → 例題とデモ
→ ビジュアライゼーション → ３次元プロット。
MATLABのヘルプ → 例題とデモ
→ 言語/グラフィクス → ３次元サーフェスプロット
例２：非格子上のデータから格子上のメッシュデータを求める
linspace：全観測点をカバーする領域の格子座標を求める
meshgrid：観測地点をカバーするメッシュを作る
griddata：観測値からメッシュ上の値を求める
x,y：観測点の座標、z：観測値 → X,Y：等間隔格子の座標、Ｚ：格子上の推定値
xlin = linspace(min(x),max(x),33);
ylin = linspace(min(y),max(y),33);
[X,Y] = meshgrid(
xlin,ylin);
Z = griddata(x,y,z,X,Y,'cubic');
mesh(X,Y,Z)
hold on
plot3(x,y,z,'.','MarkerSize',15) %一様でないデータ
例３：３角形要素分割一発表示
tridemo
他の例の詳細は、
file:////Ch1tp033_mino/研究室共通/MATLAB/jhelp/techdoc/umg/umg.html
file:////Ch1tp033_mino/研究室共通/MATLAB/jhelp/techdoc/basics/getting6.html#1030
を参照のこと
例１：１成分の波形表示（linear vslinear)
plot(x,y)
例２：１成分の波形表示（linear vslog)
semilogy(x,y)
例３：１成分の波形表示（log vslog)
loglog(x,y)
例４：複数成分の波形同時表示例
subplot(2,2,1);plot(x1,y1)
subplot(2,2,2);plot(x2,y2)
subplot(2,2,3);plot(x3,y3)
subplot(2,2,4);plot(x4,y4)
file:////Ch1tp033_mino/研究室共通/MATLAB/jhelp/techdoc/plotedit/plotedit.html
または、
file:////Ch1tp033_mino/研究室共通/MATLAB/jhelp/techdoc/umg/umg.html
の「グラフのラベリング」を参照のこと。3.3.1. すべての文字サイズを14ポイントに変えたい時
set(0,'defaultuicontrolfontsize',14);
set(0,'defaultaxesfontsize',14);
set(0,'defaulttextfontsize',14);
3.3.2.横軸をdate、ラベリングを年月日にする場合の例
% create data to label date on Jan.1 and Jun.1for everyyear
yr=1988:1990;
t = sort([datenum(yr,1,1),datenum(yr,6,1)])
plot(t,t) % draw a straight line, y=t, as an example
datetick('x')
|
RGB 値 |
省略名 |
完全名 |
|
[1 1 0] |
ｙ |
yellow |
|
[1 0 1] |
m |
magenta |
|
[0 1 1] |
c |
cyan |
|
[1 0 0] |
ｒ |
red |
|
[0 1 0] |
g |
green |
|
[0 0 1] |
b |
blue |
|
[1 1 1] |
w |
white |
|
[0 0 0] |
k |
black |
|
[1 0.5 0] |
|
橙 |
|
[0.875 0.625 0.25] |
|
金 |
camlight[right left]; camlight(az,el)
ligntingflat: 表面が平坦面のモザイク状； lighting phong:表面がなめらか
view(az,el)
for j = 1:n
plot_command
M(j) = getframe
end
movie(M)
a(20:30)=NaN 配列Ａの20番目から30番目までが、プロットされない。
例1：フィルター係数を求めてフィルターをかける
x2=filter(b,a,x1） b,aはfilterの係数
フィルター係数の計算例
４点の移動平均
a=1; b=[1/4 1/41/4 1/4]
５次バタワースＬＰＦ(fc10Hz、1000Hzサンプル）
[b,a]=butter(5,10/500);
５次バタワースＨＰＦ(fc10Hz、1000Hzサンプル）
[b,a]=butter(5,10/500,'high');
５次バタワースＢＰＦ(fc10Hz、1000Hzサンプル）
[b,a]=butter(5,[10/500 100/500]);
５次バタワースＢＳＦ(fc10Hz、1000Hzサンプル）
[b,a]=butter(5,[10/500 100/500],'stop');
例２：積分をする
X = 0:pi/100:pi;
Y = sin(x);
Z = trapz(X,Y) 答 1.9998
または、
q = quad('sin',0,pi) 答 2.0000
例３：コンボリューション・デコンボリューション
x=conv(h,x1)
[x ｒ］=deconv(h,x2) rは剰余
例４：fftをかける
X = fft(x); plot(abs(X)); plot(angle(X))
X = fft(ｘ,n) ： n点の FFT。 Xの長さが n より小さい場合は、後ろに 0 を加えて計算する。
Xの長さが n より長い場合は、後ろがうち切られる。
線型最小二乗法
観測方程式: ｛ｙ｝ = [ X ］{p} + {error} のときの
norm（y-X*p）最小の解{p}を求める方法
p = X \ y （注意： ｙ \ Xではない！）
residual = ｙ - X*p
５．ファイルの読込
例１
：各行に数値が１個のアスキーファイル
load hts.txt % hts.txtは、１行１個の値が複数行あるファイル
plot(hts) % 波形が表示される。
例２
：各行に数値が複数個のアスキーファイル
load niz.txt % niz.txtは、各行に時刻、
x,y,z成分の値が複数行あるファイル
subplot(3,1,1);plot(hts(:,1),hts(:,2)) % x成分が時間に対して表示される。
subplot(3,1,2);plot(hts(:,1),hts(:,3)) % y成分が時間に対して表示される。
subplot(3,1,3);plot(hts(:,1),hts(:,4)) % x成分が時間に対して表示される。
例３：各行に数値と文字が混在するアスキーファイル（下記例）
Type112.34 45 Yes
Type2 23.54 60 No
[names,types,x,y,answer]= textread('mydata.dat','%s%s %f%d %s');
他のvariationは、helptextreadを参照。
textscanもお試しあれ。
例４：コンマ区切りのテキストファイル a=dlmread('file.name',',','A5..D15'); Excelで読んだらセルA5〜D15に表示されるデータだけを読む。
5.1.1. ファイルのオープン・クローズ
[fname,pname] = uigetfile('*.txt')
fid=fopen([pname fname],'r');
注記： fidは、例えば、Fortranでのwrite(1,*)の1に対応する変数。
でオープンし、読み書きが終わったら
status=fclose(fid);
注記：クローズに失敗したらstatusは負の値になる。
以下のプログラムでは、ファイルの開閉の部分は省略。
5.1.2. ファイルの第１行を読み飛ばし、第２行から最後まで読むとき
b=fgetl(fid)
a=fscanf(fid,'%g %g %g %g',[4,inf]);
5.2.1. 通常
[fname,pname] =uigetfile('*.txt')
fid=fopen([pname fname],'r');
a=fread(fid,2,'float')% 最初の浮動小数点（4byte）を２つ読んでaに入れる。
b=fread(fid,'uint32') % 以下、最後まで、符号無し整数として読み、bに入れる
% その他、int32（4byte符号付き整数）、char（1byte文字）等あり。
5.2.2. Windows上のMATLABで、Unixのバイナリファイルを読むとき
fid=fopen([pname fname],'r','ieee-be');% 違うのは'ieee-be'だけ。
i =input(' 1 or 0 ?');
[x,y] = ginput: マウスをクリックした場所の座標
k = menu('Choose a color','Red','Green','Blue')
6．ＧＵＩを使ってグラフの値を読む・計算パラメータを調節する
6.1.初心者版：スペクトルのフラットレベルとコーナー周波数を読み、画面に出力する
clear
clf
echo off
fc=10;
fo=1;
n=1:256;
f=fo.*n;
amp=(1+(f./fc).^4).^(-0.5);
ans='n';
loglog(f,amp)
set(gca,'YLimMode','manual')
set(gca,'YLim',[min(amp)/10,max(amp)*10])
set(gca,'XLimMode','manual')
set(gca,'XLim',[min(f)/10,max(f)*10])
axis tight
disp('Click the corner to determine flat level and corner frequency')
disp('Hit any key to finish it')
while ans=='n'
w=waitforbuttonpress;
if w==0
cp=get(gca,'CurrentPoint');
x=cp(1,1);y=cp(1,2);
loglog(f,amp,[1,x],[y,y],[x,256],[y,max(y.*(256./x)^(-2),0.001)])
else
ans='y'
end
end
x
y
6.2.プロ（本学物理ＯＢ熊沢氏）が作ったプログラム
t=0:0.1:10;
y=sin(t);
h=plot(t,y);
getp(h);
サインカーブとカーソル（縦線）が表示される。
カーソルをドラッグすると移動ができ、しかも値が表示される
ただし、下記の内容をgetp.mというファイル名で保存しておくこと。
function getp(action)
if nargin==0;
error('Too few input arguments');
end
if(ishandle(action))
handle.h=action;
%findobj('type','line');
x=get(handle.h,'xdata');
y=get(handle.h,'ydata');
ax=axis;
xi=(ax(2)-ax(1))/2;
[temp,i]=min((x-xi).^2);
xi=x(i);
yi=y(i);
handle.h2 = line('xdata',[xi xi],'ydata',[ax(3) ax(4)],'ButtonDownFcn','getp(''down'')','erasemode','xor');
handle.h3 = line('xdata',xi,'ydata',yi,'Marker','o','erasemode','xor');
handle.h4 = uicontrol('position',[10 40 60 20],'Style','text','string',...
['x=' sprintf('%0.3g',xi)],'HorizontalAlignment','left');
handle.h5 = uicontrol('position',[10 30 60 20],'Style','text','string',...
['y=' sprintf('%0.3g',yi)],'HorizontalAlignment','left');
set(gcf,'userdata',handle);
elseif(ischar(action))
switch action
case 'down'
set(gcf,'WindowButtonMotionFcn','getp(''move'')');
set(gcf,'WindowButtonUpFcn','getp(''up'')');
case 'move'
handle= get(gcf,'userdata');
x=get(handle.h,'xdata');
y=get(handle.h,'ydata');
pt=get(gca,'CurrentPoint');
[temp,i]=min((x-pt(1)).^2);
xi=x(i);
yi=y(i);
set(handle.h2,'xdata',[xi xi]);
set(handle.h3,'xdata',xi,'ydata',yi);
set(handle.h4,'string',['x=' sprintf('%0.3g',xi)]);
set(handle.h5,'string',['y=' sprintf('%0.3g',yi)]);
drawnow
case 'up'
set(gcf,'WindowButtonMotionFcn','');
set(gcf,'WindowButtonUpFcn','');
end
end
channel = ddeinit('excel','forecast.xls'); % Excelとの通信開始
mymtx = ddereq(channel,'r1c1:r10c10'); % Excelの指定範囲から、
MATLABの配列mymtxに読み込む
rc = ddeterm(channel);% Excelとの通信終了；正常終了ならばrc=0
逆に、Excelにデータを送る例は、下記の通り：
rc = ddepoke(channel,'r1c1:r5c5', eye(5));
a=min(x(:)) 配列ｘが２次元でも３次元でも最小値スカラーを求めてくれる。
min(x)は不可。２次元のときはベクトル、３次元配列のときは行列が出力されてしまう。
10. Memory管理
問題：２ＧＢにメモリを増やしたのにOut of Memoryとなる。
原因：パソコンは、ばかでかい行列をメモリー上の連続領域に確保する。
空き容量のtotalが十分であっても、小さい空き領域が分散している場合は不可。
また、計算中に、一度確保した行列領域を広げねばならなくなった場合、そこに
別の使用中メモリがあった場合も不可。
大地主であっても、小面積の土地が分散している場合は、自分の土地に
巨大な建物を建てることができないのと同じ。
行列保存のために使える最大領域を知る方法：
\\Ch1tp033_mino\研究室共通\MATLAB\メモリ問題\dumpmemmex.dll
上記のファイルをC:\MATLAB6p1\workにコピーしてdumpmemmexを実行させる。
ごちゃごちゃと沢山表示された最後の行に、最大連続領域の容量が表示される。
メモリの虫食いを軽減するための対策１：
使わない他のプログラムは停止。
MATLABもquitして再起動させる。
MATLABをNo Javaモードで開始。
方法：ショートカットを右クリック
→プロパティ→リンク先＝「... matlab.exe -nojvm」
その後ショートカットから開始
メモリの虫食いを軽減するための対策２：
コマンドsaveで変数を一時退避。
コマンドclearの実行。たった一つのスカラー変数が、障害になっているかも。
コマンドpackの実行。メモリの虫食いを軽減してくれる。
コマンドhuge_matrix=zeros(big,big)で巨大行列領域をメモリ上に確保
コマンドloadで変数を空き領域に再配置。
上記を適宜利用すべし。