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'.
In case your have just been working on another exercise, clean up your workspace:
Add the file to your folder: 'U:\Tutorial\Matlab_scripts'
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:
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.
>> plot(n,x) % View the selected signal.
>> X = fft(x); % Perform the Fourier transform of x.
>> 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):
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:
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: