Matlab中的FFT函数示例

这个例子展示了如何使用FFT函数来分析太阳黑子活动在过去300年中的变化。太阳黑子活动是周期性的，大约每11年达到一个最大值。让我们确认一下。这是一个叫做苏黎世黑子相对数的数量图，它测量了黑子的数量和大小。天文学家已经把这个数字列了将近300年了。

工具/原料

电脑
Matlab软件

方式/步骤

1
在号令行窗口，输入如下号令：
load sunspot.dat
year = sunspot(:,1);
relNums = sunspot(:,2);
plot(year,relNums)
title('Sunspot Data')
如图1所示。
2
按“Enter键”，获得Figure1。
如图2所示。
3
以下是前50年的现状。
在号令行窗口，输入如下号令：
plot(year(1:50),relNums(1:50),'b.-');
如图3所示。
4
旌旗灯号处置的根基东西是快速傅立叶变换（FFT）。要获取太阳黑子数据的FFT，请键入以下内容。
Y的第一部门Y（1）只是数据的和，可以删除。
在号令行窗口，输入如下号令：
Y = fft(relNums);
Y(1) = [];
复平面上Fourier系数（由Y给出）的分布图很标致，但很难诠释。我们需要一种更有效的方式来查抄Y中的数据。
在号令行窗口，输入如下号令：
plot(Y,'ro')
title('Fourier Coefficients in the Complex Plane');
xlabel('Real Axis');
ylabel('Imaginary Axis');
如图4所示。
5
Y的复震级平方称为功率，功率与频率的关系图称为“周期图”。
在号令行窗口，输入如下号令：
n = length(Y);
power = abs(Y(1:floor(n/2))).^2;
nyquist = 1/2;
freq = (1:n/2)/(n/2)*nyquist;
plot(freq,power)
xlabel('cycles/year')
title('Periodogram')
如图5所示。
6
周期/年的比例有些不便利。
我们可以用年/周期来作图，估量一个周期的长度。
在号令行窗口，输入如下号令：
plot(freq(1:40),power(1:40))
xlabel('cycles/year')
如图6所示。
7
为了便利起见，我们绘制了功率与周期的关系图（此中period=1./freq）。正如预期的那样，有一个很是显著的周期，其长度约为11年。
在号令行窗口，输入如下号令：
period = 1./freq;
plot(period,power);
axis([0 40 0 2e+7]);
ylabel('Power');
xlabel('Period (Years/Cycle)');
如图7所示。
8
最后，我们可以经由过程选择最强的频率来更切确地确定周期长度。红点定位这一点。
在号令行窗口，输入如下号令：
hold on;
index = find(power == max(power));
mainPeriodStr = num2str(period(index));
plot(period(index),power(index),'r.', 'MarkerSize',25);
text(period(index)+2,power(index),['Period = ',mainPeriodStr]);
hold off;
如图8所示。