1. 程式人生 > >通過Matlab 使用 FFT 分析週期性資料

通過Matlab 使用 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')

單獨的傅立葉係數難以解釋。計算係數更有意義的方法是計算其平方幅值,即計算冪。由於一半的係數在幅值中是重複的,因此您只需要對一半的係數計算冪。以頻率函式的形式繪製功率頻譜圖,以每年的週期數為測量單位。

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

太陽黑子活動發生的最大頻率低於每年一次。為了檢視更易解釋的週期活動,以周期函式形式繪製冪圖,以每週期的年數為測量單位。該繪圖揭示了太陽黑子活動約每 11 年出現一次高峰。