User Tools

Site Tools


analysis:nsb2015:week6

Fourier series, transforms, power spectra

Goals:

  • Refresh your memory of basic trigonometry and become comfortable plotting periodic functions in MATLAB
  • Understand the central idea of the Fourier transform and implement your own decomposition and reconstruction code
  • Appreciate the fundamental issue of spectral leakage, and how analysis parameters and signal properties affect it
  • Gain a basic familiarity with different methods for spectral estimation, and use some nonparametric estimators to obtain power spectral densities (PSDs) for some interesting neural signals

Resources:

Introductory remarks

From the very first electrical recordings of the human brain (Berger, 1929) it has been clear that oscillations are a salient feature of brain activity. The functional relevance, underlying mechanisms and clinical applications of oscillatory activity continue to be active areas of research; all of these endeavours rely critically on accurate characterizations and quantifications of the data. A fundamental tool in the analysis of oscillations is Fourier analysis, which is introduced in this module.

Generating and plotting basic sinusoids

To study oscillations we require periodic functions of time, such as sin(t), which repeat its values in regular intervals.

Recall that to plot this sort of function in matlab, we first define the time axis (commonly a variable with the name tvec), pass this to the function we wish to plot as an argument, and plot the result:

%% plot a simple sinusoid
Fs = 100; % in samples per second (Hz)
t0 = 0; t1 = 1; % start and end times
tvec = t0:1./Fs:t1; % construct time axis
 
f = 2; % frequency of sine to plot
y = sin(2*pi*f*tvec); % note sin() expects arguments in radians, not degrees (see sind())
 
stem(tvec,y);

Apart from its frequency (2Hz in the above example), a sinusoid is further characterized by its phase and amplitude, which are readily manipulated as follows:

phi = pi/2;
 
figure;
y = sin(2*pi*f*tvec + phi); % a phase shift
stem(tvec,y);
hold on;
plot(tvec,cos(2*pi*f*tvec),'LineWidth',2) % notice, cosine is simply phase shifted sine
legend('Sin (Phase Shifted)', 'Cos');
hold off;
 
a = 2;
 
figure;
y = a.*sin(2*pi*f*tvec + phi); % amplitude change
stem(tvec,y);

As an aside, sinusoids are the basis for wireless transmission, such as that for FM radio, where a message is encoded in the frequency modulation of a “carrier” signal:

%% frequency modulation (FM) signal
f2 = 10;
m = 2;
 
subplot(311)
s1 = sin(2*pi*f*tvec);
plot(tvec,s1); title('message');
 
subplot(312);
s2 = sin(2*pi*f2*tvec);
plot(tvec,s2); title('carrier');
 
subplot(313);
s3 = sin(2*pi*f2*tvec + m.*sin(2*pi*f*tvec - pi/2));
plot(tvec,s3); title('FM signal');

Notice how the sinusoid in the bottom plot smoothly modulates its frequency. Many oscillations in the brain, such as hippocampal theta and ventral striatal gamma, show such modulations as well.

Sums of sinusoids and harmonic series

What happens if we sum some sinusoids together? A compact representation for this is a harmonic series:

$$S(t)=\sum_{n=1}^{N}{a_{n} cos(2 \pi n f_{0} t+\phi_{n})}$$

Note that each sinusoid of this series has a frequency that is an integer multiple $n$ of some base frequency $f_0$. This is a compact representation because we only specify the amplitude $a_n$ and phase $\phi_n$ for each of the series, with the frequencies fixed.

Figure 4.7 in the Weeks chapter is a sum of four sinusoids, with a base frequency of 2:

%%Harmonic series example
mag = [0.1 0 1.3 0.5]; % magnitudes for each term
pha = [-pi/6 0 pi 2*pi/3]; % phases for each term
f = 2; % base frequency
 
signal_out = zeros(size(tvec));
for ii = 1:numel(mag) % note, the book chapter uses i, not best practice!
 
    this_signal = mag(ii)*cos(2*pi*f*ii*tvec + pha(ii));
    plot(tvec,this_signal,'r:'); hold on;
    signal_out = signal_out + this_signal; % build the sum
 
end
figure;
plot(tvec,signal_out,'LineWidth',2);

☛ Why is it not a good idea to use the variable name i? (Hint: for the same reason you should not define a variable named pi…)

It looks like we can create some interesting signals by summing sinusoids.

In fact, the central insight underlying Fourier analysis is that we can use a sum (series) of sinusoids to approximate any signal to arbitrary precision. An important corollary to this is that any signal can be decomposed into a series of sinusoids. Let's demonstrate this.

Decomposing and reconstructing a signal

We will use the MATLAB function fft() (for Fast Fourier Transform; a technical point discussed in the Weeks chapter is that this in fact computes something called the Discrete Fourier Transform or DFT) to obtain the amplitudes and phases ($a_n$ and $\phi_n$ in the equation above) of a randomly generated signal, and then piece the signal back together by plugging them in to a harmonic series:

x = round(rand(1,8)*10); % generate a length 8 vector of integers between 0 and 10
xlen = length(x);
 
% get magnitudes and phases of Fourier series
X = fft(x);
Xmag = abs(X); % magnitudes, a_n
Xphase = angle(X); % phases, phi_n
 
n = 0:xlen-1;
t = 0:0.05:xlen-1; % a finer timescale to show the smooth signal later
 
for iH = xlen-1:-1:0 % reconstruct each harmonic
    s(iH+1,:) = Xmag(iH+1)*cos(2*pi*n*iH/xlen + Xphase(iH+1))/xlen;
    sm(iH+1,:) = Xmag(iH+1)*cos(2*pi*t*iH/xlen + Xphase(iH+1))/xlen;
    % detail: xlen appears here because the fundamental frequency used by fft() depends on this
end
 
ssum = sum(s);
smsum = sum(sm);
 
figure;
plot(n, x, 'go', t, smsum, 'b', n, ssum, 'r*');
legend({'original','sum - all','sum - points only'});

This gives something like:

Notice that the reconstruction is essentially perfect. We retrieved the original 8 points by plugging in the coefficients returned by fft() into a series of sinusoids.

☛ What is the purpose of running the for loop backwards in the above code? (Hint: recall that MATLAB slows down substantially if it needs to allocate memory for expanding variables.)

The set of magnitudes and phases that describe the harmonic series that reconstruct the signal are known as the magnitude spectrum and phase spectrum respectively. The magnitude spectrum can be interpreted as signal power at different frequencies.

Interpreting the output of MATLAB's fft() function

You might have noticed the magnitude and phase spectra above have a peculiar form. Let's explore this a little further using a signal with known frequency and phase content:

Fs = 20; % in samples per second (Hz)
t0 = 0; t1 = 1; % start and end times
tvec = t0:1/Fs:t1-(1/Fs); % construct time axis; generate exactly 20 samples
 
f = 2; % signal frequency
y = sin(2*pi*f*tvec); % construct signal, a 2Hz sine wave sampled at 20Hz for 1s
 
yfft = fft(y,length(y));
yfft_mag = abs(yfft); yfft_ph = angle(yfft);
stem(yfft_mag)

The result:

Some issues are apparent:

  • fft() did not return a frequency axis, so the output is not straightforward to interpret.
  • Whatever the coefficients returned, there is not a single peak as we would expect from the pure input signal: in fact there are two peaks.

To understand this, recall that the largest frequency that can be detected in a signal sampled at Fs is Fs/2, the Nyquist frequency. Thus, we would expect the frequency axis from fft() to go from 0 to Fs/2 at most. In addition, it turns out the Fourier transform is defined for real as well as for complex signals, and it returns spectra for both these components. Since we are only interested in real signals, we get a symmetric spectrum back (note that if we did not use abs() on the output of fft() we would get some imaginary components as well).

To construct a frequency axis that takes both these ideas into account, we can do:

Npoints = length(y);
F = [-Npoints/2:Npoints/2-1]./Npoints; % construct frequency axis
 
yfft_mag = fftshift(yfft_mag); % align output, see note below
stem(F,yfft_mag);
 
xlabel('Frequency (Fs^{-1})');

fftshift() cuts the second (complex) half of the spectrum and pastes it backwards at the beginning, so that our frequency axis is now correct; it is in units of 1 / Fs so 0.1 corresponds to the 2Hz we put in.

For analysis purposes we don't care about the complex part of the spectrum, and this is usually ignored to yield what is referred to as the “single-sided spectrum”. Because we would never actually use fft() directly to estimate spectral content (see the section on spectral estimation below) we won't do more work to get rid of it now.

Notice also the superscript in the xlabel: MATLAB can interpret basic LaTeX math symbols for figure text.

A final point about the output of fft() is that it contains the magnitude spectrum at specific frequencies. For instance, we have the 0.1*Fs point, but not 0.125*Fs! This will become important later.

Zero-padding the FFT

The above example was constructed nicely so that our signal contained exactly two full periods. This will not be true for real world signals. What happens if we don't have this perfect input?

☛ Change the tvec variable above to contain one more sample, like so:

tvec = t0:1/Fs:t1;

This means that our signal is now no longer an integer number of periods. The resulting two-sided spectrum is:

Notice that:

  1. The peaks now appear at a frequency not exactly equal to the true frequency
  2. Other non-zero components have appeared

To explain (1), recall that fft() evaluates the magnitudes of specific frequencies only. Inspection of the spectrum above indicates that we don't have a frequency bin that is exactly 0.1 (the true frequency of our signal). Let's fix that using the second, optional, argument of fft().

tvec = t0:1/Fs:t1;
nPoints = [length(tvec) 64 256 1024];
 
for iP = 1:length(nPoints) % repeat fft with different numbers of points
 
    nP = nPoints(iP);
    subplot(2,2,iP);
 
    y = sin(2*pi*f*tvec);
    yfft = fft(y,nP);
    yfft_mag = abs(yfft); yfft_ph = angle(yfft);
 
    F = [-nP/2:nP/2-1]./nP;
    yfft_mag = fftshift(yfft_mag);
    plot(F,yfft_mag,'kx',F,yfft_mag,'k');
 
    title(sprintf('%d point FFT',nP));
    xlabel('Frequency (Fs^{-1})');
 
end

This gives:

As we increase the number of points to evaluate the FFT over, we get increased frequency resolution, and indeed the peaks of the spectrum converge to 0.1 as expected. Under the hood, this is in fact accomplished by padding the input signal with zeros before the DFT is computed; doing this does not change the spectral content (the zeros are not periodic), but allows a much longer harmonic series, with a smaller fundamental frequency, to be evaluated.

A typical value to use for the number of points to evaluate the FFT is the next power of 2 (after however many samples your signal contains). This is because part of what makes the FFT fast is that it can easily divide the signal in half. But non-power of 2 values also work.

☛ What happens if you try to evaluate the FFT using a number of points smaller than that in your signal?

As we increase the number of points to evaluate the FFT, we can obtain coefficients for frequencies of arbitrary precision. For this reason, the power spectrum is often referred to as power spectral density or PSD.

FIXME Demonstrate that you get the same effect by zero-padding (rather than making the signal longer)

Spectral leakage

It is clear from the above figure that even if we increase the frequency resolution of our FFT with zero-padding, we still have an imperfect estimate of the true frequency content of our signal. In particular, the estimate around the true frequency has a nonzero width – this part of the magnitude spectrum is referred to as the main lobe – and we have nonzero spectral content for other frequencies as well (the side lobes).

To explore the source of this spectral leakage, let's make our signal longer while evaluating the FFT over the same number of points:

tvec = t0:1/Fs:t1-(1/Fs);
nRepeats = [1 2 4 8];
 
nP =  1024;
 
for iP = 1:length(nRepeats)
 
    subplot(2,2,iP);
 
    y = sin(2*pi*f*tvec);
    y = repmat(y,[1 nRepeats(iP)]); % repeat the signal a number of times
 
    yfft = fft(y,nP);
    yfft_mag = abs(yfft); yfft_ph = angle(yfft);
 
    F = [-nP/2:nP/2-1]./nP;
    yfft_mag = fftshift(yfft_mag);
    plot(F,yfft_mag,'kx',F,yfft_mag,'k');
 
    title(sprintf('%d repeats',nRepeats(iP)));
    xlabel('Frequency (Fs^{-1})');
 
end

The result:

Notice that the magnitude spectrum converges to the true frequency as we increase the length of the signal. However, for any signal of finite length, there will always be some spectral leakage. This occurs because we are effectively cutting off the signal abruptly; other than changing t