Main Content

FFT を使用した周期データの解析

ライブ スクリプトを開く

フーリエ変換を使用して、一定期間にわたる自然現象など、データの変動を解析できます。

およそ 300 年の間、天文学者はチューリッヒ相対黒点数を使用して黒点の数と大きさを表にしてきました。1700 ～ 2000 年前後におけるチューリッヒ数をプロットします。

load sunspot.dat
year = sunspot(:,1);
relNums = sunspot(:,2);
plot(year,relNums)
xlabel('Year')
ylabel('Zurich Number')
title('Sunspot Data')

黒点活動の周期性をさらに詳しく調べるために、最初の 50 年間のデータをプロットします。

plot(year(1:50),relNums(1:50),'b.-');
xlabel('Year')
ylabel('Zurich Number')
title('Sunspot Data')

フーリエ変換は、データの周波数成分を特定する信号処理の基本ツールです。関数 fft を使用して、チューリッヒ データのフーリエ変換を行います。データの和が保存されている、出力の最初の要素を削除します。出力の他の要素をプロットします。これには、実軸に関するフーリエ複素係数の鏡像が含まれます。 

y = fft(relNums);
y(1) = [];
plot(y,'ro')
xlabel('real(y)')
ylabel('imag(y)')
title('Fourier Coefficients')

フーリエ係数をそのままで解釈するのは困難です。より意味のある係数の尺度は振幅の二乗です。これはべき乗の尺度です。係数の半分は振幅内で繰り返されるため、係数の半分のべき乗を計算するだけで済みます。パワー スペクトルを、年あたりのサイクル数で測定した周波数の関数としてプロットします。 

n = length(y);
power = abs(y(1:floor(n/2))).^2; % power of first half of transform data
maxfreq = 1/2;                   % maximum frequency
freq = (1:n/2)/(n/2)*maxfreq;    % equally spaced frequency grid
plot(freq,power)
xlabel('Cycles/Year')
ylabel('Power')

最大の黒点活動は、年に 1 回より少ない頻度で行われています。周期活動をよりわかりやすく表示するには、パワーをサイクルあたりの年数で測定した期間の関数としてプロットします。このプロットによって、黒点活動が 11 年に一度ピークを迎えることがわかります。 

period = 1./freq;
plot(period,power);
xlim([0 50]); %zoom in on max power
xlabel('Years/Cycle')
ylabel('Power')

参考

| |

関連するトピック