Exploring cyclic data using FFT (FFT)

In this demonstration, we analyse data containing a chronicle of sunspot number observations on the surface of the Sun and identify cyclic components using the Fast Fourier Transform algorithm.

Introduction

The Fourier transform allows us to analyse the pattern of change in data, and of course the earliest application was to analyse physical and natural processes.

The data we will look at is a compilation of nearly 300 years of observational history. Table sunspot summarises measurements made by astronomers who have monitored the number and size of sunspots on the surface of the Sun.

Examining the data

Let's load the necessary libraries and then display a graph showing the number of sunspots recorded on the surface of the Sun between 1700 and 2000.

Pkg.add(["CSV"])

using CSV, DataFrames
gr()

sunspot = CSV.read("sunspot.csv", DataFrame)
year = sunspot[:,1];
relNums = sunspot[:,2];
plot( year, relNums, xlabel="Год наблюдения",
    ylabel="Количество пятен", title="Данные о пятнах на солнце", leg=false )

Let's zoom in on the graph and try to study more closely the cyclical nature of sunspot occurrences. To do this, we will plot the graph for the first 50 years of observations.

plot( year[1:50], relNums[1:50], c=1, marker=(:circle,3),
    leg=false, xlabel="Год надлюдения", ylabel="Количество пятен" )

As we can see, the number of sunspots varies significantly, but we can assume that there is a periodic component in the data, because the peaks on the graph are located at intervals of about 10-15 years.

Fourier transform

The Fourier transform is often a basic tool in signal processing. It allows to identify harmonic components in data (oscillatory processes). Using the function fft from the pre-installed library FFTW we convert sunspot data from the temporal domain to the frequency domain.

From the resulting vector of values we need to remove the first element, which stores the sum of the data. Graph the rest of the resulting vector. It contains the real (Re) and imaginary (Im) parts of the complex coefficients of the Fourier transform.

using FFTW

y = fft( relNums )
y[1] = NaN
scatter( y, marker=(:o, 4, :white), markerstrokecolor=:red,
    xlabel="Re(y)", ylabel="Im(y)", title="Коэффициенты Фурье", leg=false )

Fourier coefficients are difficult to interpret directly. A more informative indicator of the contribution of a particular coefficient to the process is the square of their amplitude, which can be interpreted as the power of the signal. Since in the vector of power values of each component each value occurs twice (for positive and negative frequency scales), we need to calculate the power of any of the two halves of the coefficient vector. Let us plot the power spectrum of components of different frequencies, showing the number of observed cycles in 1 year.

N = length(y);

power = abs.( y[1:Int32(floor(N/2))] ).^2; # power of first half of transform data
maxfreq = 1/2;                      # maximum frequency
freq = (1:N/2)/(N/2)*maxfreq;       # equally spaced frequency grid

plot( freq, power, xlabel="Частота циклов (количество в году)", ylabel="Мощность", leg=false )

Recall that the peaks in the spot count plot are observed much less frequently than once per year. To make the graph easier to read, let us plot it as a function of the period of the harmonic component (cycle length in years).

The graph shows that the most noticeable oscillatory process, expressed in the change in the number of sunspots, is repeated approximately once every 11 years.

period = 1 ./ freq;
plot( period, power, xlimits=(0,50), xlabel="Длина цикла (годы)", ylabel="Мощность", leg=false )

This is especially noticeable after constructing an approximate graph, on which we were able to see the most powerful components of the signal closer.

Conclusion

We have studied the most frequent and simple scenario of applying the Fast Fourier Transform (FFT) through the fft function of the FFTW library to search for harmonic components of a signal, that is, to find repetitive, cyclic processes in a signal and determine their power in the overall signal.