2007年11月21日水曜日

MATLABのFFT関数の利用方法

 高速フーリエ変換(Fast Fourier Transform, FFT)とは、離散フーリエ変換 (Discrete Fourier Transform, DFT) を計算機上で高速に計算するアルゴリズムです。
 ここで、FFTを使って、太陽黒点の周期の算出する例題を紹介します。
 太陽の黒点活動について、1700 年~1987 年までの間、1 年毎に観測したWolfer 数データが、sunspot.dat にあります。黒点の動きは、太陽黒点の活動は、約11年毎に最大になる周期をもっていることが知られていますが、フーリエ変換を用いて、黒点の動きの近似周期を求めます。
%----------------- MATLABのソース -----------------%
%--- この例題はMATLABのデモに参考したものです。---%
r = 2;
c = 2;
%過去300年にわたって太陽の黒点活動
load sunspot.dat;
year = sunspot(:,1);
wolfer = sunspot(:,2);

subplot(r,c,1);
plot(year,wolfer);
grid on;
title('太陽の黒点活動')
xlabel('年');ylabel('Wolfer数');

Y = fft(wolfer); %離散フーリエ変換
N = length(Y);
Y(1) = [];
power = abs(Y(1:N/2)).^2;
nyquist = 1/2;
freq = (1:N/2)/(N/2)*nyquist;
subplot(r,c,2);
plot(freq,power), grid on
xlabel('周期/年数')
title('ピリオドグラム')

period = 1./freq;
subplot(r,c,3);
plot(period,power), axis([0 40 0 2e7]), grid on
ylabel('パワー');
xlabel('年数/周期')

[mp,index] = max(power);
period(index)
print('-dpng','-r100','fft_sunspot.png');
%----------------- MATLABのソース -----------------%
実行の結果

0 件のコメント: