BLog

ImprintImpressum
PrivacyDatenschutz
DisclaimerHaftung
Downloads 

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 microhesph.

Before 1880 a lot of values are missing in the series, and missing values are denoted by -1.0 respectively. Since 1880 only 70 of in total 51726 rows contain -1.0, and the missing values can be easily interpolated from adjacent values on the days before and after. For the purpose of converting the date column to decimal years and the interpolation of missing values, I provide the C source of the command line tool sarconv with usage instructions and examples on GitHub - cyclaero/cagconv.

The Conversion

  1. Clone the cagconv project from the GitHub repository cyclaero/cagconv:
       cd ~/install
       git clone https://github.com/cyclaero/cagconv.git cagconv
       cd cagconv
  2. Compile sarconv.c:
       cc -g0 -O3 sarconv.c -Wno-parentheses -lm -o sarconv
  3. Download the daily time series of solar active regions from Solar Cycle Science:
       curl -O http://solarcyclescience.com/AR_Database/daily_area.txt
  4. Convert the YYYY MM DD date tuples to decimal years and write it out together with the daily sunspot areas to the TSV output file:
       ./sarconv daily_area.txt sar-1880-2021.tsv
       head sar-1880-2021.tsv
    # YYYY MM DD  Total  North  South
    # Time base:   1
    # Time unit:   d
    # Point count: 51728
    t/a	At/µhsp	An/µhsp	As/µhsp
    1880.0068306	107.0	107.0	0.0
    1880.0095628	173.0	110.0	63.0
    1880.0122951	387.0	279.0	108.0
    1880.0150273	337.0	274.0	63.0
    1880.0177596	478.0	412.0	66.0
  5. Open the resulting TSV file with your favorite graphing and/or data analysis application, for example with CVA:

We already can observe said 11-year solar cycle by the 13 equidistant peaks in the given time period, that is (2021.6 - 1880)/13 = 10.89 years per cycle. Anyway, let’s have a look at the SAR spectrum which may be obtained by discrete Fourier transformation of the SAR time series.

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.

Prerequisite FFTS

Clone the ffts sources from the GitHub repository anthonix/ffts:

   cd ~/install
   git clone https://github.com/anthonix/ffts.git ffts
   sed -e 's/CMAKE_COMPILER_IS_GNUCC/CMAKE_C_COMPILER_ID MATCHES "GNU|Clang"/g’ \
       -i ".orig" CMakeLists.txt

   mkdir build; cd build
   cmake -DDISABLE_DYNAMIC_CODE=on -DENABLE_SHARED=on ..
   make
   sudo make install clean

The Spectrum

Compile cyclasar.c on either of FreeBSD, Linux or macOS:

   cd ~/install/cagconv
   cc -g0 -O3 cyclasar.c -Wno-parentheses \
      -I/usr/local/include/ffts -L/usr/local/lib -lffts -lm -o cyclasar

Execute cyclasar with the spectrum method and the above output file of sarconv, i.e. sar-1880-2021.tsv, as input file:

   ./cyclasar spectrum sar-1880-2021.tsv spectral-sar-1880-2021.tsv
   cat spectral-sar-1880-2021.tsv

# 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 1/d is zero to half of the reciprocal time base. The cycle period at the frequency at 0.000019331890/d is 51728 d = 141,6 a, and the cycle period at 0.5/d is 2 d.

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:

Power spectrum of Solar Active Regions by Fourier transformation of its daily time series from 1880 to 2021

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 cyclasar, which has been introdued already above, may be used to digitally filter the daily SAR time series. Again the input file is the output file of sarconv. In addition, the frequency range which should pass the filter (here 0 to 0.001/d) and a blur value (here 10 %) shall be informed:

   ./cyclasar filter 0 0.001 10 sar-1880-2021.tsv filtered-sar-1880-2021.tsv
   cat filtered-sar-1880-2021.tsv

# 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.

Conclusion

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