Fourier analysis of the daily time series of Solar Active Regions
Sunspots were observed and recorded for centuries. They are caused by magnetic field flux and its number and extend vary by the 11-year solar cycle. In general, cycles cry for Fourier analysis. Now, why not, let’s do it.
The SAR file
A file with the daily time series of solar active regions can be downloaded from Solar Cycle Science. However, the date column of the series comes with calendar dates whereby years, months and days are separated by spaces, and converting this to decimal years, facilitates the further analysis. The original data file is http://solarcyclescience.com/AR_Database/daily_area.txt:
YYYY MM DD Total North South 1874 5 1 -1.0 -1.0 -1.0 1874 5 2 -1.0 -1.0 -1.0 1874 5 3 -1.0 -1.0 -1.0 1874 5 4 -1.0 -1.0 -1.0 1874 5 5 -1.0 -1.0 -1.0 1874 5 6 -1.0 -1.0 -1.0 1874 5 7 -1.0 -1.0 -1.0 1874 5 8 -1.0 -1.0 -1.0 1874 5 9 1298.0 446.0 852.0 1874 5 10 1150.0 343.0 807.0 1874 5 11 767.0 240.0 527.0 1874 5 12 389.0 0.0 389.0 1874 5 13 350.0 0.0 350.0 1874 5 14 312.0 0.0 312.0 1874 5 15 273.0 0.0 273.0 1874 5 16 314.0 0.0 314.0 1874 5 17 217.0 0.0 217.0 1874 5 18 119.0 0.0 119.0 ...
Besides the date column, the file comprises columns with the size of the active region in total and in addition the sizes on the north and the south hemisphere of the Sun. The area unit is millionths of a hemisphere, and I refer to this later in short as µhsp - pronounce
Before 1880 a lot of values are missing in the series, and missing values are denoted by
We already can observe said 11-year solar cycle by the 13 equidistant peaks in the given time period, that is
The Discrete Fourier transform
For performing the Discrete Fourier Transformation of a time series with equidistant samples - which btw. is 100 % the case of the SAR series, i.e. 51728 daily samples - we may choose among some open source FT libraries. Among these, the FFTW - Fastest Fourier Transform in the West seems to be the most popular one. However, I employed the FFTS - The Fastest Fourier Transform in the South, because its license is less restrictive, and FFTS may be freely used for almost any use case. And thanks to Bluestein’s Chirp-Z algorithm it is no more necessary that the number of points is a power of 2, for the Discrete Fourier transform being fast. With Chirp-Z, the DFT is fast as well for the present 51728 samples. Anyway, before proceeding further, FFTS shall be installed.
Clone the ffts sources from the GitHub repository anthonix/ffts:
# YYYY MM DD Total North South # Time base: 1 # Time unit: d # Point count: 51728 freq/1/d |At|/µhsp 0.000000000000 1555.316162109 0.000019331890 393.730560303 0.000038663780 15.943500519 0.000057995670 139.362426758 0.000077327560 129.681396484 ... ... 0.499922672440 0.446791768 0.499942004330 2.234769583 0.499961336220 0.763285697 0.499980668110 0.996469975 0.500000000000 0.095475748
The frequency range in units of
Open the resulting TSV file with your favorite graphing and/or data analysis application, for example with CVA and plot the spectrum:
Not quite impressive yet, but wait, let’s zoom in horizontally:
The peak at 0.000251315/d represents the solar cycle, since 1/0.000251315/365.245 = 10.89 years.
Now let's zoom in vertically:
Who would have expected this. There is a broad signal with satellite peaks around 0.0374265/d. That is a period of 26.7 d, and obviously that belongs to the rotation of the Sun.
But wait, there is more periodicity which can be observed here:
If this were a power spectrum of an electrochemical measurement (I am electrochemist), I would take most of these periods as real effects, perhaps except those indicated with grey numbers. The satellite peaks around 26.7 d (20.0 to 37.9 d) certainly relate to the differential rotation of the Sun. Most probably, the rotation of the Sun appears in the spectrum because of some aliasing with the lifetime of the sunspots. The spike, peak and hill from 12.9 to 5.48 d may actually be artefacts.
Digital Fourier Filter
Another useful application of the discrete Fourier transformation is the Fourier Filter. In the first step, the time series is transformed forward into the frequency domain. In the second step the complex values of the spectrum are manipulated. Usually the values for non-interesting frequencies are set to zero. In the third step the manipulated spectrum is transformed backward into the time domain.
The command line tool
# YYYY MM DD Total North South # Time base: 1 # Time unit: d # Point count: 51728 # ./cyclasar filter 0 0.001 10 sar-1880-2021.tsv filtered-sar-1880-2021.tsv t/a At/µhsp 1880.006835938 317.208587646 1880.009521484 317.810213853 1880.012329102 318.412831881 1880.015014648 319.016350177 1880.017700195 319.620326237 ... ... 2021.617797852 221.213151926 2021.620605469 221.810444637 2021.623291016 222.408759686 2021.625976562 223.007791899 2021.628784180 223.607770156
The high cut frequency of 0.001/d means a period of 1000 days, and 1000/365.245 = 2.74 years. We will see later why it is useful to cut away higher frequencies. The blur value of 10 % results in a smooth and not abrupt frequency cut, this would suppress artefacts, like minimum SAR values of the filtered curve are dragged significantly below zero.
Open the resulting TSV file with your favorite graphing and/or data analysis application, for example with CVA and plot the digitally filtered SAR time series:
All SAR cycles consist of double or even triple peaks. In most cases, a step can be observed at the lower tail of the peaks. Is this real or artefacts induced by the filter? Also the peak heights of the filtered curve (max. 3353 µhsp) are way lower than in the unfiltered one (max. 8382 µhsp).
To begin with, we are in the lucky situation of a perfect curve segment for Fourier analysis, since its start and end are almost in phase, and already this saves us from spurious phenomena caused by spectral leakage, therefore no windowing, which might have introduced other artefacts, is necessary and was not applied.
Let’s have a closer look at the peak of Solar Cycle 20:
The curve integrals are 3.60·106 µhsp·d (unfiltered - blue) and 3.61·106 µhsp·d (filtered: 0-0.001/d, 10 % kT - red). The difference is 100·0.01/3.6 = 0,28 %, and this does not leave much space for artefacts. Another interesting observation is, that the curve integral of all 13 cycles over the whole 141 years is 43.35·106 µhsp·d. 43.35·106/13 = 3,334·106 µhsp·d. Since the Sun is almost a perfect sphere, a hemi-sphere is half of this and the magnitude µ means 10-6, the total area of the Sun covered with sunspots in the course of a cycle sums up to 3.334/2·106·10-6 = 1.667 of its sphere’s area.
Anyway, it is clear now, that the peak split and the steps in the tail are not artefacts, it even seems that these are underestimated in the filtered curve for which the frequency cut was 0.001/d. This becomes quite evident, if we compare above result with a curve filtered with a cut of 0.002/d:
With that, the difference of the respective curve integrals is only 0.03 %. The new peaks let the filtered curve more closely follow the unfiltered one, and therefore these are not artefacts.
The Fourier analysis gives us an understanding about the periodicities of different effects which form a measurement curve. With the Fourier Filter it is possible to surgically tell apart the effects from different frequency ranges, and eventually accentuate the features of interest.
Copyright © Dr. Rolf Jansen - 2021-08-15 11:09:09
Discussion on Twitter: 1434159282409332738