Engee documentation
Notebook

Fourier transform

The Fourier Transform is a formula that is used to translate one vector of numbers into another. The first vector consists of a set of signal values taken at discrete moments in time. The second vector consists of a discrete set of values characterising the amplitude and frequency of harmonic signals, from which the observed signal can be made up by reconstructing it in time or space.

Consider a vector $x$, including $n$ measurements taken at equal time intervals. The Fourier transform of such a vector is described as follows:

$$y_{k+1} = \sum_{j=0}^{n-1}\omega^{jk}x_{j+1}$$

The angular frequency $\omega = e^{-2\pi i/n}$ is complex root of one, $i$ is imaginary unit. For both vectors $x$ and $y$, the indices $j$ and $k$ take values from 0 to $n-1$.

The function fft from the library FFTW decomposes a vector into a Fourier series using the fast Fourier transform (FFT) algorithm. Consider a sinusoidal signal $x$, each taken at equal time intervals. The moments when the measurements are taken form a vector $t$. The signal consists of two components with frequencies 15 and 20 Hz. Let's study the signal, the measurement of which was made every 1/50th of a second for 10 seconds.

In [ ]: 
Pkg.add(["WAV"])
In [ ]: 
Ts = 1/50;
t = 0:Ts:10-Ts;                     
x = sin.(2pi * 15t) .+ sin.( 2pi * 20t );
plot( t, x, xlabel="Время (с)", ylabel="Амплитуда", leg=false )
Out[0]:

Decompose the signal into a Fourier series and create a vector $f$, whose discrete components form a sample of the signal in the frequency domain.

In [ ]: 
using FFTW

y = fft(x);
fs = 1/Ts;
f = (0:length(y)-1)*fs/length(y);

Let's plot the result of transformation, where the abscissa shows the frequency of harmonic components, and the ordinate shows the amplitude.

In [ ]: 
plot( f, abs.(y), title="Амплитудный спектр", xlabel="Частота (Hz)", ylabel="Амплитуда" )
Out[0]:

The graph shows four bursts, although the original signal had two components with frequencies of 15 and 20 Hz. The point is that the right half of the graph mirrors the left half. The discrete Fourier transform of the signal in the time domain also has its periodicity - the first half of the vector corresponds to the area of "positive frequency values", and the second half - to the area of "negative frequency values" of the components. The bursts in the region of 30 and 35 Hz actually correspond to the values of -20 and -15 Hz.

In practice, it is often sufficient to use only one half of the Fourier series without performing unnecessary calculations. For this purpose you can use the function rfft. But for clarity, let's derive the whole Fourier series by shifting the right half to the region of negative frequencies using the function fftshift.

In [ ]: 
n = length(x);
fshift = (-n/2:n/2-1)*(fs/n);
yshift = fftshift(y);
plot( fshift, abs.(yshift), title="Амплитудный спектр", xlabel="Частота (Uw)", ylabel="Амплитуда" )
Out[0]:

Noisy signal

In a real-world scenario, a measured signal often includes both a useful component and random noise, which makes it impossible to unambiguously distinguish between the harmonic components. Using the Fourier transform can help cut off the random component and leave only the systematic, oscillatory process. Consider a new measurement vector xnoise, the result of adding the signal x and random Gaussian noise.

In [ ]: 
using Random
Random.seed!(0)
xnoise = x .+ 2.5 .* randn( size(t) );

In practice, the frequency dependence of signal power is used more often than the frequency dependence of amplitude. The power is equal to the square of the amplitude of the Fourier series component divided by the number of points in the series. Let's calculate and plot the power spectrum of the signal, with zero frequency at the centre of the axis $x$. Despite the presence of noise, we can clearly see both bursts on the signal power graph.

In [ ]: 
ynoise = fft( xnoise );
ynoiseshift = fftshift( ynoise );    
power = abs.( ynoiseshift ).^2 ./ n; 
plot( fshift, power, title="Спектр мощности", xlabel="Частота (Гц)", ylabel="Мощность сигнала" )
Out[0]:

Conversion rate

To calculate the Fourier transform of an entire series of $n$ elements of a vector $y$ using the original formula would require about $n^2$ floating point calculations. The fast Fourier transform algorithm requires about $n log n$ operations, which makes this method very attractive because of the time savings, especially if the signal to be processed consists of several million points. Moreover, specific implementations of the FFT algorithm can compute even faster in a number of specific cases. For example, in the case where $n$ is 2 to some degree.

Consider the audio file bluewhale.wav - a recording of acoustic communications of blue whales recorded on an underwater microphone.

In [ ]: 
using WAV
x, fs = wavread( "bluewhale.wav" );

Blue whales communicate in infrasound, the frequency band of acoustic oscillations almost does not overlap with the frequencies to which human hearing is sensitive. Therefore, the time vector in the presented data is compressed by a factor of 10 to make the height of the sound signal available to our perception.

The next cell allows you to play back the sound signal read from the file.

In [ ]: 
using Base64

buf = IOBuffer();
wavwrite(x, buf; Fs=fs);
data = base64encode(unsafe_string(pointer(buf.data), buf.size));
markup = """<audio controls="controls" {autoplay}>
              <source src="data:audio/wav;base64,$data" type="audio/wav" />
              Your browser does not support the audio element.
              </audio>"""
display( "text/html", markup );

Let's visualise this signal:

In [ ]: 
whaleMoan = x[24500:31000];
t = 10 * (0:1/fs:(length(whaleMoan))/fs);
plot( t, whaleMoan, xlabel="Время (с)", ylabel="Амплитуда", label=false, xlimits=(0,t[end]) )
Out[0]:

Now let's increase the length of the signal to the nearest number $2^N$ upwards. We will get a longer segment of the signal. And calculate the Fourier series for this new signal segment. If the signal length had not been taken from the series $2^N$ (at $N = 1,2\dots$), the function fft, would probably still augment the signal with zeros to the nearest length $2^N$. This is done to speed up the FFT algorithm considerably, especially if the value of the number of points has no small divisors.

In [ ]: 
m = length( whaleMoan ); 
n = nextpow(2, m);

whaleMoan_zeros = zeros( n )
whaleMoan_zeros[1:length(whaleMoan)] = whaleMoan

y = fft( whaleMoan_zeros );

Let us plot the power spectrum of this signal. According to the graph, the acoustic signal includes the main harmonic - the sound at a frequency of about 17 Hz, and several other noticeable harmonic components.

In [ ]: 
f = (0:n-1) .* (fs/n)/10; # вектор значений частоты
power = abs.(y).^2/n;     # спектр мощности
plot( f[1:Int(floor(n/2))],power[1:Int(floor(n/2))],
    xlabel = "Частота (Гц)", ylabel = "Мощность", leg=false )
Out[0]:

Determining the phase of sinusoids

The Fourier transform can also be used to construct the phase spectrum of the signal under study. For the sake of the example, we will create a signal in which two harmonic signals with frequencies of 15 and 20 Hz are summed. The first signal will be a cosine with a phase value of $-\pi/4$, and the second signal will be a sinusoid with no phase offset. The sampling frequency of the signal will be 100 Hz and the signal length will be 1 s.

In [ ]: 
fs = 100;
t = 0:1/fs:1-1/fs;
x = @. cos(2*pi*15*t - pi/4) - sin(2*pi*40*t);

Let's decompose the signal in Fourier series and plot the amplitude as a function of frequency.

In [ ]: 
y = fft(x);
z = fftshift(y);

ly = length(y);
f = (-ly/2:ly/2-1)/ly*fs;

plot( f, abs.(z), st=:stem, marker=(5, :white, :o),
    markerstrokecolor=1,
    xlabel="Частота (Гц)", ylabel="|y|", grid=true )
Out[0]:

Now let's look at the phase part of the Fourier series coefficients, removing the values with low amplitude. The graph shows the phase spectrum of the signal as a function of frequency.

In [ ]: 
tol = 1e-6;
z[abs.(z) .< tol] .= 0;

theta = angle.(z);

plot( f, theta/pi, st=:stem, marker=(5, :white, :o),
    markerstrokecolor=1,
    xlabel="Частота (Гц)", ylabel = "Фаза / π", grid=true )
Out[0]:

Conclusion

We have constructed several different spectra (amplitude, power and phase) based on the Fourier series expansion of the signal. The function fft, which uses the FFT algorithm, and several other functions (e.g. fftshift) helped us with this. As a result, we were able to clearly show all the properties of the input signal that we are interested in.