Engee documentation
Notebook

The Fourier transform

The Fourier transform is a formula used to convert one vector of numbers into another. The first vector consists of a set of signal values taken at individual time points. The second vector consists of a discrete set of values characterizing the amplitude and frequency of harmonic signals, from which it is possible to compose the observed signal by reconstructing it in time or space.

Consider the vector , which includes measurements taken at the same time intervals. The Fourier transform of such a vector is described as follows:

Angular frequency is a complex root of единицы, – imaginary unit. For both vectors and , indexes and they take values from 0 to .

Function fft from the library FFTW performs the decomposition of a vector into a Fourier series using the fast Fourier transform (FFT) algorithm. Consider a sinusoidal signal , each taken at equal time intervals. The moments of measurement form a vector . The signal consists of two components with frequencies of 15 and 20 Hz. Let's study the signal, which was measured every 1/50 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 , the discrete components of which 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 the transformation, where the frequency of the harmonic components is derived from the abscissa, and the amplitude is derived from the ordinate.

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 fact is that the right half of the graph mirrors the left half. The discrete Fourier transform of a signal in the time domain also has its own periodicity – the first half of the vector corresponds to the region of "positive frequency values", and the second half corresponds to the region of "negative frequency values" of the components. The bursts in the region of 30 and 35 Hz actually correspond to the values -20 and -15 Hz.

In practice, it is often sufficient to use only one half of the Fourier series without performing unnecessary calculations. To do this, you can use the function rfft. But for clarity, we will output the entire Fourier series, 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, the measured signal often includes both a useful component and random noise, which makes it impossible to clearly distinguish the harmonic components. Using the Fourier transform can help cut off the random component and leave only a 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 dependence of signal power on frequency is used more often than the dependence of amplitude on frequency. The power is equal to the square of the amplitude of the component of the Fourier series divided by the number of points in the series. Let's calculate and plot the signal power spectrum, with zero frequency in the center of the axis . Despite the noise, we can clearly see both spikes 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 the entire series from vector elements Using the initial formula, you need to perform the following steps: calculations with floating-point numbers. The fast Fourier transform algorithm requires the following order of magnitude This makes this method very attractive because of the time savings, especially if the signal being processed consists of several million points. Moreover, specific implementations of the FFT algorithm can be calculated even faster in a number of specific cases. For example, in the case when is equal to 2 to some extent.

Consider an audio file bluewhale.wav – 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 vibrations 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 so that the pitch of the audio signal is accessible to our perception.

The next cell allows you to play the audio 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 );

Visualize 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 signal length to the nearest number. in the ascending direction. We will get a longer segment of the signal. And calculate the Fourier series for this new segment of the signal. If the signal length was not taken from the series (when ), then the function fft most likely, I would have supplemented the signal with zeros to the nearest length anyway. . This is done in order to significantly speed up the FFT algorithm, especially if the value of the number of points does not have 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's plot the power spectrum of this signal. According to the graph, the acoustic signal includes the main harmonic – 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]:

Phase detection of sinusoids

Using the Fourier transform, it is also possible to construct the phase spectrum of the studied signal. For the sake of example, we will create a signal in which two harmonic signals with a frequency of 15 and 20 Hz are combined. The first signal will be a cosine with a phase value , and the second signal will be a sine wave with no phase shift. The sampling frequency of the signal will be 100 Hz, the signal length is 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);

Decompose the signal into a 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 coefficients of the Fourier series, while removing the values with a 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 (amplitudes, powers, and phases) based on the decomposition of the signal into a Fourier series. The function helped us with this. fft, using the FFT algorithm, and several other functions (for example fftshiftAs a result, we were able to visually show all the properties of the input signal that we are interested in.