Calculating the period of the sunspot cycle

Two weeks ago, I used the sunspot number data provided by the Solar Physics Group at NASA’s Marshall Space Flight Center to demonstrate positioning plots in window. This week, I’d like to show how to calculate the period of the sunspot cycle.

If you haven’t already done so, download the sunspot numbers file and place it in your IDL path. Read it with the astrolib READCOL procedure:

file = file_which('spot_num.txt', /include)
readcol, file, year, month, sunspots

Next, transform the sunspot series to the frequency domain and compute magnitude and power spectra:

mspec = abs(fft(sunspots))
pspec = mspec^2

(Aside: FFT: it’s all you need.)

I’d like to display the power spectrum as a function of frequency. This requires a few statements to set up a frequency vector based on the time data from the sunspot numbers file:

sampling_interval = 1/12.0 ; years, from file
freq_nyquist = 0.5/sampling_interval ; years^{-1}
n_sunspots = n_elements(sunspots)
freq = findgen(n_sunspots/2)/(n_sunspots/2-1)*freq_nyquist

The sampling interval is one month. This allows me to determine the Nyquist frequency (the maximum resolvable frequency), which along with the number of sunspots in the file, gives me the information to define a frequency vector. Note that I’m defining only positive frequencies up to the Nyquist frequency; this is because I’m going to display a one-sided power spectrum. Fold the negative Fourier modes into the spectrum, omitting the zero (DC) mode:

freq_nodc = freq[1:n_sunspots/2-1]
pspec_nodc = 2*pspec[1:n_sunspots/2-1]

and display the result on logarithmic axes:

p = plot(freq_nodc, pspec_nodc, /xlog, /ylog, $
   xtickunits='exponent', $   ; IDL 8.2.1
   xtitle='Frequency ($yr^{-1}$)', $
   ytitle='Spectral Density', $
   title='Power Spectrum of Sunspot Numbers, 1749-2010')

Note that I can set the style of tick units on the axes: ‘exponent’ and ‘scientific’ are new in IDL 8.2.1.

The power spectrum has a peak near 0.1 yr^{-1}. Locate the peak frequency with MAX and mark it on the plot with POLYLINE:

pspec_nodc_peak = max(pspec_nodc, i_peak)
freq_nodc_peak = freq_nodc[i_peak]
!null = polyline([1.0,1.0]*freq_nodc_peak, p.yrange, color='orange', /data)

The inverse of the peak frequency gives the period of the sunspot cycle: approximately 11 years. Mark this on the plot with TEXT:

sunspot_peak_period = 1.0/freq_nodc_peak
str = '$T_{peak}$ = ' + strtrim(sunspot_peak_period,2) + ' years'
!null = text(1e-1, 1e-4, str, /data)

Here’s my result:

Power spectrum of sunspot numbers time series. Data courtesy NASA.

Finally, just because it’s interesting, test Parseval’s theorem (energy is conserved in the time- and frequency-domain representations of the signal) in several forms:

IDL> print, total(pspec)*n_sunspots, total(sunspots^2)
1.46437e+007 1.46437e+007
IDL> print, mspec[0], mean(sunspots)
52.0133      52.0134
IDL> print, total(pspec[1:*]), stddev(sunspots)^2*(n_sunspots-1)/n_sunspots ; convert from sample variance
1965.67      1965.66

Conservation of greatness (of FFT)!


About Mark

I solve scientific programming and visualization problems with IDL.
This entry was posted in data access, data analysis, visualization and tagged , , , , , . Bookmark the permalink.

5 Responses to Calculating the period of the sunspot cycle

  1. Pingback: Calculating the period of the sunspot cycle | The IDL Data Point | Solar Flare 2012

  2. Reblogged this on asubsetofdaves and commented:
    Nice to see solar data being used to illustrate a fundamental tool of the astronomer: the FFT :o)

  3. Paul Hiemstra says:

    Thanks for this excellent post. I think, however, that there is a bug in the code, shouldn’t:

    freq = findgen(n_sunspots/2)/(n_sunspots/2-1)*freq_nyquist


    freq = findgen(n_sunspots/2-1)/(n_sunspots/2-1)*freq_nyquist
    ^^^ ; extra -1

  4. Simon says:

    Perhaps you ought to use a periodogram, not an fft. Try Lomb’s.

  5. mohammed says:

    please possible i know any things about subcycle for daily sunspots data >>thanks

Leave a Reply

Fill in your details below or click an icon to log in: Logo

You are commenting using your account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s