www.gusucode.com > signal 案例源码程序 matlab代码 > signal/FindLocalMaximaInDataExample.m
%% Find Peaks in Data % Use |findpeaks| to find values and locations of local maxima in a set of % data. % % The file |spots_num.mat| contains the average number of sunspots observed % every year from 1749 to 2012. The data are available from NASA. % % Find the maxima and their years of occurrence. Plot them along with the % data. % Copyright 2015 The MathWorks, Inc. load(fullfile(matlabroot,'examples','signal','spots_num.mat')) [pks,locs] = findpeaks(avSpots); plot(year,avSpots,year(locs),pks,'or') xlabel('Year') ylabel('Number') axis tight %% % Some peaks are very close to each other. The ones that are not recur at % regular intervals. There are roughly five such peaks per 50-year period. % % To make a better estimate of the cycle duration, use |findpeaks| again, % but this time restrict the peak-to-peak separation to at least six years. % Compute the mean interval between maxima. [pks,locs] = findpeaks(avSpots,'MinPeakDistance',6); 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) %% % It is well known that solar activity cycles roughly every 11 years. Check % by using the Fourier transform. Remove the mean of the signal to % concentrate on its fluctuations. Recall that the sample rate is measured % in years. Use frequencies up to the Nyquist frequency. Fs = 1; Nf = 512; df = Fs/Nf; f = 0:df:Fs/2-df; trSpots = fftshift(fft(avSpots-mean(avSpots),Nf)); dBspots = 20*log10(abs(trSpots(Nf/2+1:Nf))); yaxis = [20 85]; 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)]) %% % The Fourier transform indeed peaks at the expected frequency, confirming % the 11-year conjecture. You also can find the period by locating the % highest peak of the Fourier transform. [pk,MaxFreq] = findpeaks(dBspots,'NPeaks',1,'SortStr','descend'); Period = 1/f(MaxFreq) hold on plot(f(MaxFreq),pk,'or') hold off legend('Fourier transform','1/meanCycle','1/Period') %% % The two estimates coincide quite well.