Fourier transformation is used to decompose time series signals into frequency components each having an amplitude and phase. Using the inverse Fourier transformation the time series signal can be reconstructed from its frequency-domain representation. Fourier transformation is one of the most important concepts in digital signal processing and is not only used for estimating the spectral distribution of a signal in the frequency domain (the power spectrum). Fourier transformation is also the foundation of coherence analysis and certain types of so-called surrogate signals. Finally, the Fourier transformation is implemented in many DSP (Digital Signal Processing) routines because any mathematical operation in the time domain has an equivalent operation in the frequency domain that is often computationally faster. Thus, Fourier transformation is occasionally implemented solely to speed up algorithms. Using the inverse Fourier transformation, the time-domain signal is reconstructed from its frequency domain representation.
In this tutorial we will study the formula for computing the discrete Fourier transform (DFT) and numerically study the DFT on a short signal (only a few samples) in order to keep track on the indices in the FT formula (which most people consider complicated and abstract when working with long signals). One “problem” with large time series is that you get as many frequency components in the frequency domain as you have data samples in the time domain (!).
Once you have understood this simple exercise, you can easily invent a slightly more complicated example (longer signal and multiple frequency components) or intuitively accept that the idea also applies to longer signals.
As in the other tutorials, we ask you to type the commands instead of copying directly from this page, because part of learning Matlab and the content of the formulas introduced here is to get “hands-on experience”…!
Preparation for the tutorial
The commands in this box ensure that the following exercise has all the necessary scripts and does not interfere with variables or graphics from previous exercises. If you come directly from another exercise it should be clear which commands you can skip.
Start MATLAB and find your current working directory using 'pwd'.
In case your have just been working on another exercise, clean up your workspace and remove the figures:
clear all % Deletes all variables in your workspace.
close all % Closes all figure windows.
Go to a directory where you can find you the files you need to download and save, e.g., 'U' using 'cd'.
Make a directory “Tutorial” using 'mkdir'.
Go to this directory using
and make a directory 'Matlab_scripts'.
To run the rest of this exercise, you have to download a MATLAB function (right-click, and select “Save link as”)
For generating a sine wave of a certain frequency, length, sampling frequency and phase, you may use 'SignalGenerator_sin'. Generate one and a half cycle of a 0.5 Hz oscillation:
clear all % Clean your workspace just in case…
close all % Close any figure windows.
frq = 0.5; % The frequency.
len = 3; % The length of the wave in seconds.
Fs = 1000; % Sampling frequency of the wave.
phs = 90/360; % the phase of the wave.
[wave] = signalgenerator_sin(frq, len, Fs, phs); % Generate a sine wave.
'wave' is a sine wave sampled at high-frequency, which you may consider an approximation of “the original continuous process”, and is used here for visualization only.
Recall that the sampling theorem allows us to sample a signal at just two times per cycle of the highest frequency component in the signal. Thus, we sample the 'wave' at times 0, 1, 2, and 3 s giving a digitized 'signal' of N = 4 samples:
n = [0 1 2 3]; % Sampling times (in seconds) of the 0.5 Hz signal.
x = [1 -1 1 -1]; % Values of 'wave' at times 'n'.
N = 4; % Number of samples in 'x'.
figure(1) % Activate figure window "1" with the sine wave.
hold on % Keep the sine wave in subsequent plotting to figure 1.
plot(n, x, 'r') % Plot the digitized 'x' to see that it tracks the 'wave'.
plot(n, x, 'rs') % The red squares mark the digitized samples.
axis([-0.5 3.5 -1.5 1.5]) % Define axes to better visualize signals.
We now apply the Discrete Fourier Transform (DFT) to the signal in order to estimate the magnitude and phase of the different frequency components. The DFT formula is:
X : the frequency domain representation of signal time-series signal 'x'. The formula yields one complex number X[k] for every k. k : the k'th frequency component; k = 0,1,…N-1. N : the total number of samples in signal 'x'. x : the time series signal (the data); n = 0,1,…N-1. n : the n'th sample (in the time domain). j : the imaginary unit.
is the k'th basis function.
To understand this equation better - you need to know about complex numbers, and that a signal can be described by a complex number, z.
z consists of : A the amplitude
Note that in the summation over n = 0, 1, … N-1, the value of the basis function () is computed (“sampled”) at the same times 'n' as your recorded signal x[n] was sampled.
Thus, the DFT formula basically states that the k’th frequency component is the sum of the element-by-element products of 'x' and '', which is the so-called inner product of the two vectors and , i.e.,:
Because is a unit-length vector, the inner product of x and may be seen (somewhat abstractly…) as a “projection of the data x onto the complex basis function“.
In order to directly apply the DFT formula to our 'x', we therefore compute the basis vectors ('') corresponding to the value of the basis functions at times 'n' and our value of 'N' (use the “arrow-up key” to speed up typing):
We distinguish between “basis functions” and “basis vectors” such that the functions refer to the mathematical formula whereas the vectors refer to the output of the functions when having entered a given set of input variables ('n', 'k', and 'N'). The elements of the basis vectors are complex and for any element (time point or sample):
1) 'real(e_k)' (the real part of e_k) corresponds to the cosine wave.
2) 'imag(e_k)' (the imaginary part of ek) corresponds to the sine wave.
3) 'angle(e_k)' corresponds to the phase of e_k.
4) 'abs(e_k)' being the modulus (or norm or magnitude).
The basis vectors should be thought of as the “basic wave shapes” that we project our measured signal onto and, thus, the time series signal is reconstructed by superposition of the basic waves with appropriate weighting factors (amplitudes).
Try to visualize the basis vectors or “basic wave shapes” [Use “arrow-up key” to retrieve previous MATLAB commands and speed up plotting]:
plot(n, real(e_0)) % The real part of e_0.
plot(n, imag(e_0),'r') % The imaginary part of e_0.
In case you prefer a different scaling you can run this little loop (once you get used to MATLAB, this is a good example of a visualization action that would take a lot more time in other software or if you used a graphical user interface):
for i=1:4; subplot(4,1,i); grid; axis([-0.5 3.5 -1.5 1.5]); end
The blue lines are the real part of the waves; the red lines are the imaginary part.
Note that the lowest frequency (k = 0) cosine wave, makes zero complete cycles over the N samples, i.e., it is a DC signal.
The Fourier transform ('X') of our sampled signal 'x' at frequency components k = 0, 1, 2, 3, is now obtained as the sum of the element-by-element products of 'x' and the sampled basis functions 'e_k', which you just computed above:
'X' contains four complex numbers—one for each of the four frequency components (k = 0, 1, 2, 3). The complex numbers have an amplitude (absolute value or “norm”) and a phase relative to the first sample in the time series. The amplitudes of the four frequency components (k = 0, 1, 2, 3) are:
Questions: Do you see the logic of which basis vector has non-zero amplitude?
Plot a small amplitude spectrum:
k = [0 1 2 3];
The power spectrum is more commonly used in neurophysiology than the amplitude spectrum and is simply “the spectrum of squared absolute values of the Fourier transform vectors” you calculated above (although in practice some additional normalization is applied, which we shall not be concerned with at this stage):
Question: In comparing the amplitude and the power spectrum, do you see advantages and disadvantages of using either of the two?
FT was developed more than two hundred years ago, but is slow in its original form as implemented above. In the 1960s, the Fast Fourier Transform (FFT) was developed, which speeds up computations by a factor of 100–1000 times. All Fourier transformations in MATLAB are based on FFT, we shall not cover the mathematical tricks to make Fourier Transform a Fast Fourier Transform, but just use it to cross-check that the build-in MATLAB function 'fft' produces the results we obtained above:
X_fft = fft(x)
The formula for transforming a signal in the frequency-domain back into the time domain is similar to the DFT formula, except that now we calculate the samples in the time domain one-by-one as a summation over all the frequency components 'k':
Let's reconstruct the original signal from the frequency-domain representation:
x_reconstruct = ifft(X_fft)
or on the basis of the “manually computed 'X':
x_reconstruct = ifft(X)
As you will notice, the original signal is recovered! This was a very simple and short signal.