過去ログ一覧 |
フェムトセカンドトップへ |
http://melonsode.fem.jp/
Matlabで自己相関関数を求める方法がわかった。というか、パワースペクトルとの関係がわかった。どうも自己相関関数を一発で求める関数はないらしい。で、以前、柳川さんの言ってたとおり、フーリエ変換を使えば良かったのだ。もとの関数をまずフーリエ変換。それからそれにそれの複素共役をかけて、つまり二乗してパワースペクトルを求める、そうしてそれを逆フーリエ変換する。これで自己相関関数のできあがり。Matlabで、と書いたが実際に使っているのはoctaveというフリーソフト。これはMatlabのフリー版とも言えるソフトでかなり高機能らしい。Matlabの本を見ながらこのoctaveというソフトの使い方を勉強している。図は、試しに作ったノイズ入り正弦波（上）とその自己相関関数（下）。
以下追記(2007/9/29)：
上記のノイズ入り正弦波の自己相関関数計算のスクリプトは以下の通り。
N = 1024;
tmax = 0.5;
t = linspace(0,tmax,N);
y =sin(2*pi*100*t) + randn(1,N);
subplot(2,1,1); plot(t,y);axis([0,tmax/5,-2,2])
Y=fft(y);power=Y.*conj(Y);
auto_y=real(ifft(power)); auto_y=auto_y./max(auto_y);
subplot(2,1,2);plot(t,auto_y); axis([0,tmax/5,-2,2])
以下、フーリエ変換に関するコマンド。
フーリエ変換：Y=fft(X), Y=fft(X, n)
これはXをフーリエ変換しYに代入。X,nのnは配列数でnが２のべき乗(2^n)の場合は省略可能。
あとは同様に、
逆フーリエ変換：Y=ifft(X), Y=ifft(X, n)
二次元フーリエ変換：Y=fft2(X), Y=fft2(X, m, n)
二次元フーリエ逆変換：Y=ifft2(X), Y=ifft2(X, m, n)
ゼロ周波数位置の中心移動：Y=fftshift(X)
参考文献：「信号処理」「画像処理」のためのMATLAB入門【増補版】、高井信勝・著
octaveはなっかなか感じを掴めないでいたが、今日、すこしわかった。AとBという行列があって、普通の行列の積を計算したいときはA*Bと書けばよい。では、「A.*B」と書いたらどうなるか？そのピリオドは何を意味するか？正解は、対応する要素同士をかけ算する、ということ。これは行列のかけ算の算法とは違う。Matlabでは多次元配列を行列として扱うために、行列としての演算と、多次元配列としての演算を分ける必要がある。それで、ピリオドが必要になるのだ。また、y=sin(t)とやるだけで、各tに対するyの値を計算してくれるということに気づき、それがとても便利なことに（やっと）気付いた。もしＣ言語なりなんなりでやるんだとしたら、
for(t=0;t＜T_MAX;t++){
y[t]=sin(t)
}
みたいな風に書かなきゃならない。それが、y=sin(t)と書くだけで計算されるのはかなり簡便で、記述密度が高いと言える。Matlabではなんでも行列で表現し、一次元配列を使いたいときにも、一行Ｎ列の行列ということになる。これをちゃんと実感するまでにちょっと時間がかかった。明日からばりばりMatlab、もとへoctaveを使っていくぜい！
This web site names フェムトセカンド ／ femto second = 10-15 秒 .
Author: 野澤真一 / Nozawa Shinichi
since 2006/4/1