在数据中查找峰值

在数据中查找峰值

使用 findpeaks 求出一组数据中局部最大值的值和位置。

% 文件 spots_num.mat 包含从 1749 年到 2012 年每年观测到的太阳黑子的平均数量。
% 这些数据可以从美国国家航空航天局获得。
% 
% 求出最大值及其出现的年份。将它们与数据一起绘制出来。
load('spots_num.mat')

[pks,locs] = findpeaks(avSpots);
% disp([pks,locs]);

figure
plot(year,avSpots,year(locs),pks,'or')
xlabel('Year')
ylabel('Number')
axis tight

  一些峰值彼此非常接近。有些峰值不会周期性重复出现。每 50 年大约有五个这样的峰值。

为了更好地估计周期持续时间，请再次使用 findpeaks，但这次将峰间间隔限制为至少六年。 计算最大值之间的间隔均值。

% 一些峰值彼此非常接近。有些峰值不会周期性重复出现。每 50 年大约有五个这样的峰值。
% 
% 为了更好地估计周期持续时间，请再次使用 findpeaks，但这次将峰间间隔限制为至少六年。
% 计算最大值之间的间隔均值。
[pks,locs] = findpeaks(avSpots,'MinPeakDistance',6);

figure
plot(year,avSpots,year(locs),pks,'or')
xlabel('Year')
ylabel('Number')
title('Sunspots')
axis tight

legend('Data','peaks','Location','NorthWest')

cycles = diff(locs);
meanCycle = mean(cycles);
disp(meanCycle); % meanCycle = 10.8696

 众所周知，太阳活动周期大约为 11 年。我们使用傅里叶变换来检验这一周期。

% 众所周知，太阳活动周期大约为 11 年。我们使用傅里叶变换来检验这一周期。
% 减去信号的均值以重点关注其波动。
% 采样率以年为单位。使用从零到奈奎斯特频率的频率。
Fs = 1;
Nf = 512;

df = Fs/Nf;
f = 0:df:Fs/2-df;
trSpots = fftshift(fft(avSpots-mean(avSpots),Nf));% fftshift：将零频分量移到频谱中心

dBspots = 20*log10(abs(trSpots(Nf/2+1:Nf)));

yaxis = [20 85];
figure
plot(f,dBspots,1./[meanCycle meanCycle],yaxis)
xlabel('Frequency (year^{-1})')
ylabel('| FFT | (dB)')
axis([0 1/2 yaxis])
text(1/meanCycle + .02,25,['<== 1/' num2str(meanCycle)])

% 傅里叶变换确实在预期频率处出现峰值，证实了 11 年的估计值是准确的。
% 您还可以通过定位傅里叶变换的最高峰值来求得周期。
[pk,MaxFreq] = findpeaks(dBspots,'NPeaks',1,'SortStr','descend');

Period = 1/f(MaxFreq) % Period = 10.8936

hold on
plot(f(MaxFreq),pk,'or')
hold off
legend('Fourier transform','1/meanCycle','1/Period')

% 这两个估计值非常吻合。

 这两个估计值非常吻合。