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.
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:
Signal basics
First, let's create a sine wave. The formula for this is:
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, 
, provides the angle based on frequency (F) and time (t).
, 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.
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.
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).
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.
Spatial aliasing produces ghosting or 'wrap-around' artifacts in MRI.
Explore: Effects of sampling rate on aliasing
Now let's look at the effect of down-sampling the signal to a new, lower-frequency sampling rate Fs.
You can change the slider values and it will re-generate plots with the new sampling rate.
Fs_high = 1000; % Sampling frequency, e.g., Hz (samples/sec)
T = 1/Fs_high; % Sampling period (time between samples, e.g., sec)
L = Fs_high * 5; % Length of signal (in samples) for 5 sec signal
t = (0:L-1)*T; % Time vector
y = sin(2 * pi * F * t + theta);
Fs = 20; % New sampling frequency
Fs = 12
Fs = 12
plot_signals_with_downsampling(y, t, Fs, F)
The top panel shows the original signal at 1000 Hz (black) and signal sampled at the new, downsampled rate (pink). If the lines are on top of one another, there's no aliasing due to the reduced sampling rate. But if the pink line appears to be moving at a slower frequency, it is aliased.
Explore: Effects of phase on aliasing
Now let's look at the effect of changing the phase (theta).
You can change the slider values and it will re-generate plots with the new phase.
Fs_high = 1000; % Sampling frequency, e.g., Hz (samples/sec)
T = 1/Fs_high; % Sampling period (time between samples, e.g., sec)
% length of signal in sec
num_sec = 1
num_sec = 1
% theta in radians
theta = 1.7
theta = 1.7000
% New sampling frequency
Fs = 12
Fs = 12
L = Fs_high * num_sec; % Length of signal (in samples)
t = (0:L-1)*T; % Time vector
y = sin(2 * pi * F * t + theta);
plot_signals_with_downsampling(y, t, Fs, F)
What do you notice about the aliasing effects as you change the phase (theta)? Why?
Which of the values that you can adjust in the sliders impact aliasing? Describe the impact of each.
What happens as you truncate the timeseries to make it short, as opposed to very long? Why does this happen?
More Exercises
The Nyquist frequency is 1/2 the sampling rate. Any signal with frequencies higher than that (faster oscillation) will be aliased back into the low frequency range.
In these examples, the sampling frequency is very high (e.g., 1000 Hz), and the underlying frequency of the sin wave is very low. So there is no problem with aliasing. Now let's look at a case where we're sampling a faster-moving underlying frequency with a lower sampling rate.
Fs = 50; % Rate we're sampling at: Sampling frequency, e.g., Hz (samples/sec)
T = 1/Fs; % Sampling period (time between samples, e.g., sec)
L = 150; % Length of signal (in samples)
t = (0:L-1)*T; % Time vector
% Sine wave parameters
theta = 0; % In radians
F = 20; % Frequency of underlying signal we're detecting (cycles/sec, e.g., Hz)
y = sin(2 * pi * F * t + theta);
create_figure('sin wave', 1, 2)
subplot(1, 2, 1)
plot(t, y)
xlabel('Time (sec)')
title('Time domain')
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
yf = 2 * yf(1:L/2+1); % Take only the half of frequencies below the nyquist
f = Fs*(0:(L/2))/L; % define frequencies for each element of yf
plot(f, yf)
xlabel('Frequency (Hz)')
title('Frequency domain')
nyquist = Fs ./ 2; % Any signal faster than Nyquist freq will be aliased
h = plot_vertical_line(nyquist);
set(h, 'LineStyle', ':')
What is the frequency of the underlying signal?
What is the sampling rate for this signal?
What is the Nyquist frequency?
Would we expect aliasing?
How can I tell from the FFT plot that the signal is/is not aliased?
What signal frequencies (Fs) would you expect to produce aliasing based on the Nyquist frequency? 10, 15, 20, 25, 30, 40...?
At what sampling frequency would the ability to detect a 30 Hz signal break down and start to produced aliased lower-frequency signals? Modify the code to try it and test your answer.
Subfunctions
These functions simplify the code above, making some specialized plots, etc. You can run them by running code blocks in the Live Script editor.
function plot_signals_with_downsampling(y, t, Fs, F)
downsample_by = round(1000/Fs);
figure;
plot(t, y, 'k-'); hold on
set(gcf, 'Position', [0 644 1358 231])
plot(t(1:downsample_by:end), y(1:downsample_by:end), '.-', 'Color', [1 .3 .6], 'LineWidth', 2);
Nyquist = Fs/2;
title(sprintf('Sampling rate = %3.1f, Nyquist limit = %3.2f', Fs, Nyquist))
figure; hold on;
[myfft, freq, line_han] = fft_calc(y, 1/1000);
set(line_han, 'Color', 'k', 'LineWidth', 2, 'LineStyle', '-');
[myfft, freq, line_han] = fft_calc(y(1:downsample_by:end), 1/Fs);
set(line_han, 'Color', [1 .3 .6], 'LineWidth', 2, 'LineStyle', '-');
set(gca, 'XLim', [0 (max(Fs, 1.2 * F))]);
end