Matlab Tips

Links:

Matlabでアニメーションgifを作成

Matlabでは直接アニメーションgifは作れないようなので,一旦tif形式で1枚1枚のファイルに書き出して, image magic付属の(gradsにもついている)convertでファイルの結合・gifへの変換・アニメーション化,を行う. サンプルプログラムを示そう.

tsz=100
dim=ones(20,10);
xs=[1:1:20]'*ones(1,10);
ys=ones(20,1)*[1:10];
k=2*pi/10; 
l=2*pi/5; 
omega=2*pi/10; 

for t=1:tsz
  dim=sin(k * xs + l * ys - omega * t);
  tifffile=sprintf('d:\\tmp\\%02d.tif',t)
  contour(dim',[-0.8 -0.6 -0.4 -0.2 0.2 0.4 0.6 0.8])
  print(tifffile,'-dtiff','-r75')
  hold off;
end
cd d:\tmp
pwd
pause
!convert -loop 0 -delay 20 *.tif anim.gif % loopを指定しないとpower pointでanimationしない.

3D graphで鉛直上向きに減少する座標を使う方法

海洋の深度や大気の気圧は,鉛直上向きに減少ししかもマイナス記号をつけないのが一般的だ.これをMatlabの3D Graphで実現するには,
1) 妥協してマイナス記号をつける
2) 図の下向きが,物理空間の鉛直上向きで一旦図を描き,Figureをクリックしてz軸のプロパティの「降順」をクリックする.
の2通りがある.(1)の方が安直なので自分で見るには適しているが,最終的な図を作るには(2)がよい.

マニュアルに載っていないlineの属性

以下を実行して見て取れるように,マニュアルには載っていないけれど(と思う),線の太さなどもplotコマンドで指定できる. Matlabの線はデフォルトでは細すぎて,よく見えないことが多いので,適宜太くするとよいだろう.

plot(t,sin(2*t),'-mo',...
                'LineWidth',2,...
                'MarkerEdgeColor','k',...
                'MarkerFaceColor',[.49 1 .63],...
                'MarkerSize',12)

テンソルっぽい計算

Matlabでは行列の計算は定義されているが,テンソルになるとダメである.例えば,3次元配列(2×2×2)と1次元配列の掛け算はやってくれない(エラーになる).もっと簡単な,

x=ones(2,2,2); y=ones(2,1);
z=x(1,1,:)*y;
つまり(1×1×2)と(2×1)をやらせようとしても
??? エラー: ==> *
行列入力はサポートしていません.
となってしまう.このような場合は,reshapeを使って,行列の演算にしてやることで計算ができる.つまり,
z=reshape(x(1,1,:),1,2)*y;
とすれば良い.

%cと%sはsscanfでは異なり、sprintfでは同じ意味

sscanfで空白を含む文字を読み出すと、%cであれば空白が保存されるが、%sでは除去される。 通常空白を保存したいだろうから、%c(%sではなく)を使う方が良い。書き出しはどっちでも良い。次の例で動作の違いが分かるだろう。

cstr='a b c'; 
sscanf(cstr,'%s')
sscanf(cstr,'%c')
sprintf('%s',cstr)
sprintf('%c',cstr)

ans =
abc

ans =
a b c

ans =
a b c

ans =
a b c

NaNからユーザー指定の欠損値への変更

Matlabの欠損値 NaN を含むデータを、他のソフトウェアで利用するには(例えばGradsでの描画)、NaN を適当なユーザーが指定した定数に変換してやらなくてはならないことがしばしばある。この変換は、次のようにしてmatlab の isnan関数を使うことで可能である。

x=[1 NaN 3];
x(isnan(x))=0

x =
     1     0     3

ちなみにmatlabではNaNとの==,~=などの評価は機能しない。たとえば、

x=[1 NaN];
x==NaN
x~=NaN

ans =
     0     0
ans =
     1     1

のように、NaN との==による比較は常に偽、~=での比較は常に真が返される。したがって、 NaN以外の値であれば、x(x==3)=0 などとして変換ができるが、NaNについてはそうはいかずisnanを使わざるを得ないことを記憶しておくべきである。なお、無限大になる可能性のある演算を含んでいる場合は、NaNとともにinfを除外するために、isfinite関数を使うこともできる。ただし、無限大になる可能性のある演算は本来行うべきではないので、isnanをもちいる方がゼロ割の間違いを発見しやすい。

複素行列の転置

単純な転置'では、複素数は転置に加えて共役が取られる。転置だけしたい場合は、.'を用いること。例えば

x=[1+i 2+2*i; 3+3*i 4+4*i];
x'
x.'

を実行すると

ans =
   1.0000 - 1.0000i   3.0000 - 3.0000i
   2.0000 - 2.0000i   4.0000 - 4.0000i
ans =
   1.0000 + 1.0000i   3.0000 + 3.0000i
   2.0000 + 2.0000i   4.0000 + 4.0000i
となる。

列ベクトルに簡単に変換する方法

xが列ベクトル、行ベクトルのどちらでも、以下の操作で必ず列ベクトルになる。

x = x(:);               

ただし、(xsz,yxz)の行列も、(xsz*ysz,1)の列ベクトルに変換されるので、行列への 利用は慎重に。

file名の流用

[path,name,ext] = fileparts(outctlname)
outctl.grdname=sprintf('^%s.grd',name)

関数で欠損値が定義されていない場合になんとかする方法

関数の引渡し変数にundefを入れておいて,それが定義されていないことを次のようにしてチェックし,適宜な値を入れる.

if exist('undef')~=1; undef=-1.e+10; end

簡単なエラーメッセージの出力方法

重要なのはメッセージを出してそこで処理を止めること.

たとえば、
if (strcmpi(inctl.tstp(2:3),'MO') & strcmpi(inctl.tstp(2:3),'YR')
    '*testfun* unit of tstp is not MO nor YR',var,pause
end

書式

引数が整数の場合
d引数を10進値として出力する
u引数を符号無し10進値として出力する
o引数を符号無し8進値として出力する
x引数を符号無し16進値として出力する
引数が浮動小数点の場合
f引数をddd.ddd(dは10進数の1桁)の形式で出力する
e引数をd.ddde+ddの書式で出力する

配列の複製

配列の複製を作成するには、
A(1:nsiz,1)*ones(1,msiz) 
あるいは
ones(nsiz,1)*A(1,1:msiz)
というようように、onesを用いて行列の掛け算をする方法と、
A(1:nsiz,ones(1,msiz))
あるいは
A(ones(nsiz,1),1:msiz)
のように行列の中にonesを書く方法とがある。一般に後者の 法が、速度が速いようである。ちなみに
A(ones(n1siz,n2siz),1:msiz)      %n1siz,n2sizを別に指定しても
A(ones(n1siz*n2siz,1),1:msiz)    %n1siz*n2sizのみが意味あるなら、別々に指定するべきではない
の二つで得られる結果は変わらない。最初の例は、プログラムの可読性を損なうので、行うべきではない。

関数の操作方向

多くの関数は第一要素について演算を行う。

sum,prod,mean,std,median

第一要素について演算を行うことが、

size(sum(randn(3,2)))
size(prod(randn(3,2)))
size(mean(randn(3,2)))
size(median(randn(3,2)))

の答えがすべて、[1 2] となることで確認できる。

sort

第一要素について演算を行うことが、
x=[3 7 5; 0 4 2]; sort(x)

とすると答えが、

     0     4     2
     3     7     5

となることで確認できる。

FFTの操作方向

第一要素について演算を行う。 fft のhelp には「FFTは、行列に対しては列単位で操作を行います。」となっている。 これは、dim(xsz,ysz)という配列があれば、xに関して(yを固定して)演算を行うということである。 これはtaper関数(hamming等)が[nsiz,1]の大きさの結果を返すことと整合的である。

% 列について変数
x=randn(1,8);           
x(2,:)=x;               
f=fft(x);      
max(abs(f(1,:)-f(2,:)))     %ゼロならf(1,:)とf(2,:)は同一

% 行について変数
x=randn(8,1);
x(:,2)=x;
f=fft(x);
max(abs(f(:,1)-f(:,2)))      %ゼロならf(:,1)とf(:,2)は同一

出力引数が無い場合に、絵を出す関数を作る

関数中で

if (nargout==0)
  plot(x1,y1)
end

のようにすると、出力引数が無い場合に画面に図を出す関数ができる。

大計Matlab toolbox

/opt/matlab/toolbox/stats

岡田 直資 作,アニメーションtiff生成スクリプト (4/13受け)

岡田です。

まず、Matlab で 3D の図が書かれている tiff ファイル群を作るスクリプト
(2つ)をお送りします。

・図を描くのに必要なデータ
    データを収めていたマシンが修理から帰ってきたばかりなので、データを
    引き出せません。今日の夜にはデータを出せると思います。
・tiff ファイルから GIF Animation を作るプログラム。

は必要でしょうか。お知らせください。

----------------------------------------------------------------------
スクリプト 1

%cd d:\kanamoto\matlab\bd200
clear
for I=121:174
  J=0;
  II=28
  JJ=35
  isz=20; i1=1+2*(II-1); i2=i1+isz-1; 
  jsz=20; j1=1+2*(JJ-1); j2=j1+jsz-1; 
  ksz=20; k1=1 ; k2=k1+ksz-1; 
  imsh=2;jmsh=2;kmsh=2;
  cmax= 2.5E-04;
  casename='bd20b'  
%  filename=sprintf('%s%03d','ut.',I);
  filename=sprintf('%s%s%03d',casename,'ut.',I);
  fid=fopen(filename,'r','ieee-be');rdsize=[128*128*20];
  ur=fread(fid,rdsize,'float32'); ur=reshape(ur,128,128,20); 
  fclose(fid);
%  filename=sprintf('%s%03d','vt.',I);
  filename=sprintf('%s%s%03d',casename,'vt.',I);
  fid=fopen(filename,'r','ieee-be');
  vr=fread(fid,rdsize,'float32'); vr=reshape(vr,128,128,20); 
  fclose(fid);
%  filename=sprintf('%s%03d','wt.',I);
  filename=sprintf('%s%s%03d',casename,'wt.',I);
  fid=fopen(filename,'r','ieee-be');
  wr=fread(fid,rdsize,'float32');
  wr=reshape(wr,128,128,20); 
  fclose(fid);
%  filename=sprintf('%s%03d','o.',I)
  filename=sprintf('%s%s%03d',casename,'o.',I)
  fid=fopen(filename,'r','ieee-be');
  omegar=fread(fid,rdsize,'float32'); omegar=reshape(omegar,128,128,20); 
  fclose(fid);
  omegzr=zeros(size(ur));
  omegzr=omegaz(ur,vr);
  for K=k1:k2
    u(:,:,K)=ur(i1:i2,j1:j2,K)'; 
    v(:,:,K)=vr(i1:i2,j1:j2,K)'; 
    w(:,:,K)=wr(i1:i2,j1:j2,K)'; 
    omega(:,:,K)=omegar(i1:i2,j1:j2,K)'; 
    omegz(:,:,K)=omegzr(i1:i2,j1:j2,K)'; 
  end
  clear *r
  omga=zeros(size(omega)); omgb=zeros(size(omega));
  omga(omegz>0)=omega(omegz>0);   omga(omegz<0)=-1; 
  omgb(omegz<0)=omega(omegz<0);   omgb(omegz>0)=-1; 
  cmin= cmax;
  imsh1=i1*0.05;imsh=imsh*0.05;imsh2=i2*0.05;
  jmsh1=j1*0.05;jmsh=jmsh*0.05;jmsh2=j2*0.05;
  kmsh1=k1*0.05;kmsh=kmsh*0.05;kmsh2=k2*0.05;
  x=((i1:i2)-0.5)*0.05;y=((j1:j2)-0.5)*0.05;z=((k1:k2)-0.5)*0.05;
  clf
  p = patch(isosurface(x, y, z, omgb, cmin));
  set(p, 'FaceColor', [0,0,1], 'EdgeColor', 'none');
  isonormals(x,y,z,omgb, p)
  p = patch(isosurface(x, y, z, omga, cmax));
  set(p, 'FaceColor', [1,0,0], 'EdgeColor', 'none');
  isonormals(x,y,z,omga, p)
  daspect([1 1 1]);
  [cx cy cz] = meshgrid(imsh1:imsh:imsh2,jmsh1:jmsh:jmsh2,kmsh1:kmsh:kmsh2);
  h2=coneplot(x,y,z,u,v,w,cx,cy,cz,2);
  set(h2, 'FaceColor', 'green', 'EdgeColor', 'none');
%  xlim([min(x) max(x)]); ylim([min(y) max(y)]); zlim([min(z) max(z)]); 
  xlim([(i1-1.00000001)*0.05 i2*0.05]); ylim([(j1-1.00000001)*0.05 j2*0.05]); zlim([0.0 1.0]); 
  grid on;
  for ID=1:1
%    dtheta=15*(ID-1)
%    dphi=15*sin((ID-1)*pi/12)
%%    dtheta=37.5
%%    dphi=60
    view(3)
    camva(11);
    lighting gouraud;
    lit=[1/ID^2,1/ID^2,1/ID^2];
    light('Position',campos,'Color',lit);
%    camorbit(dtheta,dphi);
%    stringx=sprintf('%s%d%s%s%+10.1E%s','\bf ',I,' hour  ','渦度=',cmax,'s^-^1');
    stringx=sprintf('%s%d%s','\bf ',I,' hours');
    title(stringx,'Fontsize',12)
    xlabel(' x (km)'); ylabel(' y (km)'); zlabel(' z (km)'); 
%    axis tight;
    box on;
%    camva(11);
    J=J+1;
    tifffile=sprintf('%s%03d%03d%s','bd20butca',I,J,'.tif')
    print(tifffile,'-dtiff','-r75')
  end
%  cla
%  close all
%  clear p
%  clear o*
%  clear
end

----------------------------------------------------------------------
スクリプト 2

以下のスクリプトは omegaz.m という名前にしておく。

function omegz=omegaz(u,v);
  omegz=zeros(size(u));
  omegz(2:end-1,2:end-1,:)=v(3:end,2:end-1,:)-v(1:end-2,2:end-1,:)-(u(2:end-1,3:end,:)-u(2:end-1,1:end-2,:)); 
  omegz( 1 ,2:end-1,:)=v(2,2:end-1,:)-v( end ,2:end-1,:)-u( 1 ,3:end,:)+u( 1 ,1:end-2,:);
  omegz(end,2:end-1,:)=v(1,2:end-1,:)-v(end-1,2:end-1,:)-u(end,3:end,:)+u(end,1:end-2,:);
  omegz(2:end-1, 1 ,:)=v(3:end, 1 ,:)-v(1:end-2, 1 ,:)-u(2:end-1,2,:)+u(2:end-1, end ,:);
  omegz(2:end-1,end,:)=v(3:end,end,:)-v(1:end-2,end,:)-u(2:end-1,1,:)+u(2:end-1,end-1,:);
  omegz=omegz/100.0; 

二つのパネルでcontourf levelを同じにする方法

-----------Minobe の質問 6 Mar 1998 --------------- Hi all: I have a problem to produce two panels drawn by contourf routine in one figure using the same color convention. In order to use the same color convention I specified color levels in each panels. The color for the minimum become common for two plots, but colors for the maximum is different. This may cause a wrong visual impression. The problem can be reproduced by the following commands. I use Matlab 5.1 for Window 95.
peak1=peaks;
peak2=peak1./3.;
subplot(2,1,1);
cint=1;
lvmn=floor(min(min(peak1))/cint)*cint ;
lvmx=ceil(max(max(peak1))/cint)*cint;
clevs=lvmn:cint:lvmx;
contourf(peak1,clevs);
colorbar;
shading flat;
subplot(2,1,2);
contourf(peak2,clevs);
colorbar;
shading flat;

Thank you very much in advance. Shoshiro Minobe minobe@neptune.sci.hokudai.ac.jp Graduate School of Science, Hokkaido University, Japan ---------------------回答 6 Mar 1998--------------------------------- Hello, I haven't found a perfect solution to your problem, but here is a possible workaround:

peak1=peaks;
peak2=peak1./3.;

% make a copy of the last column
lastrow=peak2(:,end);

%append four copies of the last column to the data

peak2(:,end+1:end+4)=lastrow(:,ones(1,4));
 peak2(end,end)=max(max(peak1));

subplot(2,1,1);
cint=1;
lvmn=floor(min(min(peak1))/cint)*cint ;
lvmx=ceil(max(max(peak1))/cint)*cint;
clevs=lvmn:cint:lvmx;
contourf(peak1,clevs);
colorbar;
shading flat;
subplot(2,1,2);
[cout hand]=contourf(peak2,clevs);

% get the XLIM property of the subplot
xlims=get(gca,'xlim');

%hide the last four columns by changing the xlim property
set(gca,'xlim',[xlims(1) xlims(2)-4]);

colorbar;
shading flat;

  If you have any questions, please send email
  (日本語OK) to me 
  or for general MATLAB support to techmat@cybernet.co.jp

Regards,

Rob Henson
Cybernet Systems
MATLAB Distributor in JAPAN

rob@cybernet.co.jp

研究用ソフトウェア に戻る