The Neurophysiological Biomarker Toolbox (NBT)

# Fourier transform and reconstruction of real signals

Approximate time needed for this tutorial: 20 minutes

The commands on this page ensure that the following exercise has all the necessary scripts and data, and that the variables or graphics from previous exercises do not interfere with the present tutorial. 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'.

>> pwd

In case your have just been working on another exercise, clean up your workspace:

>> clear all
>> close all

Go to the directory 'U' using 'cd'.

>> cd U:\

Make a directory “Tutorial” using 'mkdir'.

>> mkdir('Tutorial')

Go to this directory using

>> cd('Tutorial')

and make a directory 'Matlab_scripts'.

>> mkdir('Matlab_scripts')

and place it in the directory 'U:\Tutorial'.

You also need an extra function Viewer_1ch.m : viewer_1ch.m

Then use the 'addpath’ command to inform Matlab where to find the scripts:

>> addpath U:\Tutorial\Matlab_scripts

Go to the directory with the data and load the data using the 'load' command.

>> load MEG_data

The term “workspace” refers to all the variables currently available for viewing, manipulating, saving or clearing. To view the content of your workspace use the command 'whos', which should return the following:

>> whos
Name               	Size	                    Bytes  Class

LSM		1x360000              2880000  double array
LSM_noise	1x360000              2880000  double array
RSM		1x360000              2880000  double array
RSM_noise	1x360000              2880000  double array
Trig		1x360000              2880000  double array

Grand total is 1800000 elements using 14400000 bytes

The variables 'RSM' and 'LSM' (“right and left sensorimotor cortex”) are vectors containing the 360.000 samples obtained with an MEG sensor over right and left sensorimotor cortex, respectively, during a 20 min recoding at a sampling frequency of 300 Hz. The median nerve was stimulated every 3 s with a 0.2-ms current pulse, which caused the thumb to twitch and is known to evoke pronounced neuronal processing in somatosensory cortex. 'RSM_noise' and 'LSM_noise' are the same as 'RSM' and 'LSM' except that the recordings were performed without a subject in the MEG device to estimate the noise level in the laboratory and instruments. The variable 'Trig' is a trigger channel containing a square pulse at times of stimulation, which we do not need for the present tutorial.

It is STRONGLY recommended to work in units of seconds in variables even though the actual computations are performed on samples. Working with correct units allows your intuition to detect mistakes in processing data. Therefore, we define the sampling frequency 'Fs' (a widely used abbreviation for “sampling frequency”) and use units of seconds when defining time variables:

>> Time_window = 3;			% Size of viewing window in seconds.
>> Fs = 300;				% The sampling frequency.
>> Viewer_1ch(RSM, Fs, Time_window);	% Brings up a figure window with a 3-second segment of the signal RSM.

Pressing an arbitrary key will advance the viewed signal to the next segment of length 'Time_window'.
IMPORTANT: Press 'Ctrl-c' to terminate the viewer and return to the command line. Moreover, you may want to stretch the figure window with the mouse to use the full width of the monitor or maximize the window to see more details.

Let's select a 1-second signal that looks “interesting” in terms of analyzing its frequency content and reconstructing it from its frequency-domain representation.

>> x = (RSM(32.6*Fs:33.6*Fs-1));	% Define signal of interest.
>> n = 0:1/Fs:1-1/Fs;			% Define a time vector.
>> figure(1)
>> plot(n,x)				% View the selected signal.
>> X = fft(x);				% Perform the Fourier transform of x.
>> figure(2)
>> plot(abs(X(1:end/2)).^2)		% Plot the power spectrum.
>> N = length(x);			% Define 'N' to be the number of samples.

Now we can identify the dominant frequency components from the power spectrum in Figure 2 and manually try to reconstruct the signal in Figure 1. Use the “zoom button” in the figure window to see the frequency components with peaks. A quick left-mouse-button double-click returns the original plot.

Clearly, component k=12 is important. It corresponds to the strong oscillatory component that you can visually identify in the signal 'x'. Let's calculate the k=12 basis vector (i.e., the value of the basis function at the times that 'x' is sampled):

>> k=12; e_12 = exp(j*2*pi*n*Fs*k/N);
>> figure
>> plot(real(e_12))

This is a sine wave, but to see its contribution when reconstructing 'x' from 'X' we need to weigh (multiply) the basis vector with the 12th component of 'X', cf. the inverse DFT formula:

We do this and inspect the wave:

>> figure
>> plot(1/N*real(X(12)*e_12))

Note that the phase has changed because of the non-zero phase of Fourier component 'X(12)'.

The exercise is now to reconstruct as much of the time-series signal 'x' using as few frequency components 'k' as possible. We can do this by defining a 'size(X)' vector of zeros and adding more and more elements from 'X'. Plotting the real part of the inverse Fourier transform of 'X_select' will reveal how well a few selected frequency components suffice in reconstructing the signal:

>> X_select = zeros(size(X));		% Initialize vector of zeros.
>> X_select([12]) = X([12]);		% Add Fourier component 12 to 'X_select'.
>> figure
>> plot(real(ifft(X_select)),'r')

Note, the signal is the same as the “manually” computed signal using the basis vector above. Now add the next component, e.g., k = 1:

>> X_select([1 12]) = X([1 12]);	% Add Fourier components 1 and 12 to 'X_select'.
>> figure
>> plot(real(ifft(X_select)),'r')

Question 1: What happened to the reconstructed signal?

Question 2: Can you understand why?

Now compete with your neighbor: Who can get the best match between the reconstructed signal and the original signal using a maximum of 8 components?

Question 3:
Add a title (use 'help title' to see how) with your name and the numbers of the k components you chose. Save the figure in jpeg format and e-mail it to: klaus.linkenkaer<>cncr.vu.nl

If there is time left:
Question 4:
Can you think of an objective measure to quantify how well your reconstructed signal resembles the original signal?