Signal basics, fft, and aliasing in Matlab

This lab explores basic aspects of sin/cos waves and plotting in Matlab. It then explores the Fourier Transform (FFT) a bit. The FFT is probably the most important transformation in signal processing.
More detailed explanations can be found on Michael X. Cohen's excellent Youtube channel.
Here is a screen shot that shows the relationship between time domain signals (top row) and frequency domain signals (bottom row) for some simple sine waves:
Cohen_sin_waves.png

Signal basics

First, let's create a sine wave. The formula for this is:
sin wave
Fs = 1000; % Sampling frequency, e.g., Hz (samples/sec)
T = 1/Fs; % Sampling period (time between samples, e.g., sec)
L = 1000; % Length of signal (in samples)
t = (0:L-1)*T; % Time vector
 
% Sine wave parameters
theta = 0; % Phase, in radians
F = 10; % Frequency of sin wave (cycles/sec)
 
y = sin(2 * pi * F * t + theta);
 
create_figure('sin wave');
plot(t, y)
xlabel('Time (sec)')
Try changing F to 20. What would you expect to happen? What happens?
Theta (θ) is often used to reference an angle, as it is here. θ in the equation above is the offset, or phase shift. The other part of the equation, sin wave, provides the angle based on frequency (F) and time (t).
Try changing theta, θ. What happens with positive θ? Negative θ?
θ is specified in radians, with radians equal to 360 degrees (a full circle), or πequal to 180 degrees, a semicircle.
Try replacing sin with cos above, to generate a cosine wave. What is the value of cos(θ), where theta is the angle, at time 0 (0 radians?)
Try typing sin(0) and cos(0).
What theta would shift the sin wave by 90 degrees, so it looks like the cosine wave?

The FFT

The function fft( ) estimates the coefficients of the Fourier transform, transforming a time-domain signal (i.e., an observed fMRI time series) into a series of sine waves with different amplitudes and phases.
Fourier coefficients are complex numbers, with values in the real and imaginary planes. If you are intrigued, check out this video for a great walk-through. The magnitude of this signal is what we're mainly interested in plotting here, and it relates to the power at a particular frequency. That's the length of the coefficient in complex space, as shown in the image below, and is extracted using abs(yf), where yf are the estimated Fourier coefficients. The phase is the angle of the vector as it moves through the complex plane.
Cohen_complex_numbers.png
We can think of sine waves as projections onto one axis of a vector traveling around a circle, like the hand of a clock. If we're looking at the position of the clock hand only on one axis, it will make a sine wave-type pattern. In the case of the Fourier transform, the two dimensions are the real and imaginary number lines. We can think of Fourier coefficients describing a 3-d sine wave that spirals through time (t) with frequency f. The higher the frequency, the tighter the spiral. The larger the magnitude, the larger the spiral's radius is (its path is farther from the origin).
Cohen_complex_sine_waves.png
Image credits: Michael X. Cohen
create_figure('sin wave', 1, 2);
plot(t, y)
xlabel('Time (sec)')
 
subplot(1, 2, 2)
 
yf = fft(y); % Fourier transform
yf = abs(yf/L); % When X is complex, ABS(X) is the complex modulus (magnitude) of the elements of X
 
f = Fs*(0:(L-1))/L; % define frequencies for each element of yf
 
plot(f, yf)
 
nyquist = Fs ./ 2; % Any signal faster than Nyquist freq will be aliased
 
h = plot_vertical_line(nyquist);
set(h, 'LineStyle', ':')

Aliasing

When a signal is sampled too slowly, it produces a phenomenon known as aliasing, in which an underlying higher-frequency signal appears to be a lower-frequency one. This occurs because the sample points capture values of the true high-frequency signal at systematic points in the phase of the true signal, as illustrated in the figure below.
Let's say we have a true signal that consists of a sine wave at F = 10 Hz, originally sampled at Fs = 1,000 Hz. This is shown in blue. What happens if we sample the signal at Fs = 20 Hz? We observe a sample every 1/2 cycle of the original wave, resulting in a constant value for every sample! This is shown by the black line in the figure. The 10 Hz signal has been aliased to a constant value.
Let's say we sample the signal at 15 Hz, shown in the center panel in the figure. Now, we have still not captured the frequency of the actual signal -- the observed samples appear to be moving more slowly in time. They are still aliased to a lower temporal frequency.
Finally, let's sample at 40 Hz. Now we capture the underlying waveform and its frequency very accurately, without aliasing.
What is happening is that the frequency of the signal is being reflected around a particuar value -- 1/2 the sampling frequency. This is called the Nyquist frequency, or "folding frequency". If we sample at 15 Hz, for example, the Nyquist frequency is 1/15 = 7.5 Hz. Any signal faster than that, e.g., our 10 Hz signal, will be reflected back to a lower frequency. As the 10 Hz signal is 2.5 Hz above the Nyquist, the aliased signal will have a frequency of 2.5 Hz below the Nyquist, or 5 Hz.
Fs = 1000; % Sampling frequency, e.g., Hz (samples/sec)
T = 1/Fs; % Sampling period (time between samples, e.g., sec)
L = Fs * 5; % Length of signal (in samples) for 5 sec signal
t = (0:L-1)*T; % Time vector
 
y = sin(2 * pi * F * t + theta);
 
figure;
 
subplot(3, 1, 1);
plot(t, y); hold on
 
Fs = 20; % New sampling frequency
downsample_by = round(1000/Fs);
 
plot(t(1:downsample_by:end), y(1:downsample_by:end), 'k.-', 'LineWidth', 2);
Nyquist = Fs/2;
title(sprintf('Sampling rate = %3.1f, Nyquist limit = %3.2f', Fs, Nyquist))
 
subplot(3, 1, 2);
plot(t, y); hold on
 
Fs = 12; % New sampling frequency
downsample_by = round(1000/Fs);
 
plot(t(1:downsample_by:end), y(1:downsample_by:end), 'k.-', 'LineWidth', 2);
Nyquist = Fs/2;
title(sprintf('Sampling rate = %3.1f, Nyquist limit = %3.2f', Fs, Nyquist))
 
subplot(3, 1, 3);
plot(t, y); hold on
 
Fs = 40; % New sampling frequency
downsample_by = round(1000/Fs);
 
plot(t(1:downsample_by:end), y(1:downsample_by:end), 'k.-', 'LineWidth', 2);
Nyquist = Fs/2;
title(sprintf('Sampling rate = %3.1f, Nyquist limit = %3.2f', Fs, Nyquist))
 
This example shows aliasing in time, but aliasing can occur in space as well. For example, when a moving wheel on a car is sampled at a lower frame rate in a movie, the spokes of the wheel sometimes appear to be moving backwards. Another example is the Moire effect. When gratings of parallel lines and spaces are offset from one another, a lower-frequency pattern of light and dark bars appears due to spatial aliasing.