Frequency Analysis Tools and Methods

In this notebook, you will learn tools and methods to calculate the power spectrum of a signal. A power spectrum shows the magnitude of oscillations at different frequencies in a signal.

Setup

Import Libraries

import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
from scipy.fft import fft, fftfreq

Utility Functions

Run the cell below, but feel free to ignore this code — it’s some utility functions we’ll use later for convenience.

def get_fft(data, time_vec):
    """
    Calculates the Fourier transform and frequencies of the transform for a signal
    """

    # Step 1: calculate the power spectrum by taking the fourier transform of the signal
    fft_result = fft(data)

    # Step 2: Get all frequencies for which the power spectrum was calculated in step 1
    freqs = fftfreq(len(data), d=(time_vec[1]-time_vec[0])) 

    # Step 3A: Extract positive frequencies only
    positive_freqs = freqs[freqs > 0]
    # Step 3B: Extract power spectrum for positive frequencies only
    fft_result_positive = fft_result[freqs > 0]

    # Step 4: Get real part of spectrum to get amplitude and normalize it
    psd_norm = np.abs(fft_result_positive)/np.abs(fft_result_positive).max()
    
    return positive_freqs, psd_norm

def plot_signal(sig, time_vec):
    """Makes plot of signal"""
    plt.figure()
    plt.plot(time_vec,sig)
    plt.title('Signal')
    plt.xlabel('Time')
    plt.ylabel('Amplitude')

class utils:
    get_fft = get_fft
    plot_signal = plot_signal

Section 1: Calculate the Total Frequency Content of a Signal with the Fourier Transform

The Fourier transform is a standard way of calculating the power spectrum of a signal. Here, you will apply the Fourier transform to different toy signals to get a better intuition for the relationship between peaks in the power spectrum and oscillations in your original signal.

Exercises

Example: Calculate and plot the power spectrum for the sinusoid signal oscillating at 5 Hz generated below using the Fourier transform.

# provided

# Generate sine signal
time_vec = np.linspace(0,1,50) # seconds
frequency = 5 # Hz
sine = np.sin(2*np.pi*frequency*time_vec)

utils.plot_signal(sine, time_vec)

freqs, psd = utils.get_fft(sine, time_vec)

plt.plot(freqs,psd)
plt.xlabel('Frequency (Hz)')
plt.ylabel('Power Spectral Density (normalized)')
plt.title('Power spectrum of signal')

Exercise: Calculate and plot the power spectrum for the sinusoid signal generated below oscillating at 20 Hz using the Fourier transform.

Solution

# Generate signal
time_vec = np.linspace(0,1,200) # seconds
frequency = 20 # Hz
sine = np.sin(2*np.pi*frequency*time_vec)

utils.plot_signal(sine, time_vec)

Exercise: Calculate and plot the power spectrum for the composite signal consisting of two sinusoids using the Fourier Transform. The signal is generated for you in the first cell below, you just need to run it. Note that the variable containing the composite signals is called composite_signal when you calculate the power spectrum with the utils.get_fft function.

Which frequencies do the two sinusoids that make up the composite signal oscillate at?

Hint: If the spectral peaks are all to the left in the plot, you can limit the x-axis to focus on them by adding the line plt.xlim([lower_limit, upper_limit]) (insert numbers for lower_limit and upper_limit).

# provided

# Create composite signal
sampling_frequency = 1000
time_vec = np.linspace(0,1,sampling_frequency) # seconds
frequencyA = 10 # Hz
sineA = np.sin(2*np.pi*frequencyA*time_vec)

frequencyB = 20 # Hz
sineB = np.sin(2*np.pi*frequencyB*time_vec)

composite_signal = sineA+sineB

utils.plot_signal(composite_signal, time_vec)

Solution

freqs, psd = utils.get_fft(composite_signal, time_vec)

plt.plot(freqs,psd)
plt.xlabel('Frequency (Hz)')
plt.ylabel('Power Spectral Density (normalized)')
plt.title('Power spectrum of signal')
plt.xlim([0, 30]);

Exercise: Calculate and plot the frequency spectrum for the chirp signal (chirp_signal) generated below that is gradually increasing its frequency from 1 to 15 Hz.

# provided

sampling_frequency = 1000
time_vec = np.linspace(0,1,sampling_frequency) # seconds

# generate chirp signal
chirp_signal = signal.chirp(time_vec, f0=1, f1=15, t1=1)

utils.plot_signal(chirp_signal, time_vec)

Solution

freqs, psd = utils.get_fft(chirp_signal, time_vec)

plt.plot(freqs,psd)
plt.xlabel('Frequency (Hz)')
plt.ylabel('Power Spectral Density (normalized)')
plt.title('Power spectrum of signal')
plt.xlim([0, 50]);

Exercise: Calculate and plot the frequency spectrum for the composite signal generated below consisting of a sine wave with a 30 Hz frequency and a chirp signal increasing from 1 to 15 Hz using the get_fft function.

# provided

# Create composite signal
sampling_frequency = 1000
time_vec = np.linspace(0,1,sampling_frequency) # seconds

# set frequency of sine wave and generate it
frequency = 30 # Hz
sineA = np.sin(2*np.pi*frequency*time_vec)

# generate chirp signal
chirp_signal = signal.chirp(time_vec, f0=1, f1=15, t1=1)

# create composite signal by summing sine wave and chirp signal
composite_signal = sineA+chirp_signal

utils.plot_signal(composite_signal, time_vec)

Solution

freqs, psd = utils.get_fft(composite_signal, time_vec)

plt.plot(freqs, psd)
plt.xlabel('Frequency (Hz)')
plt.ylabel('Power Spectral Density (normalized)')
plt.title('Power spectrum of signal')
plt.xlim([0, 50]);

Section 2: Welch’s method

One way to reduce the noise in your power spectrum is to use Welch’s method to calculate the spectrum. In Welch’s method, your signal is divided into segments and a power spectrum is computed for each individual segment using the Fourier transform. The final power spectrum for the whole signal period is then calculated by averaging across the power spectra for the individual segments. If you have divided your signal in to 10 segments, the final spectrum will be the average of 10 segment spectra. The length of each segment can be set with a parameter in the function from scipy called signal.welch that you will use to calculate the spectrum with Welch’s method.

Exercises

Exercise: Use Welch’s method to calculate the frequency spectrum for the signal generated below. Plot the result.

# Generate signal
sampling_frequency = 1000 # sampling frequency (Hz)
duration = 1 # seconds
time_vec = np.linspace(0,duration,sampling_frequency) # seconds
frequency = 20 #
pure_sine = np.sin(2*np.pi*frequency*time_vec)

# generate noise
noise = np.random.normal(0,1,len(time_vec))

# add noise to sine wave
sine = pure_sine + noise

utils.plot_signal(sine,time_vec)

Solution

# calculate power spectral density with Welch's method
freqs, psd = signal.welch(sine, fs = sampling_frequency)

# normalize
psd = psd/psd.max()

plt.plot(freqs, psd)
plt.xlabel('Frequency (Hz)')
plt.ylabel('Power spectral density (normalized)');

Exercise: Use Welch’s method to calculate the frequency spectrum for the composite signal below consisting of two sine waves with frequencies 10 Hz and 20 Hz and . Plot the result.

# provided

# Create composite signal
sampling_frequency = 256
time_vec = np.linspace(0,1,sampling_frequency) # seconds
frequencyA = 10 # Hz
sineA = np.sin(2*np.pi*frequencyA*time_vec)

frequencyB = 20 # Hz
sineB = np.sin(2*np.pi*frequencyB*time_vec)

composite_signal = sineA+sineB

utils.plot_signal(composite_signal, time_vec)

Solution

# calculate power spectral density with Welch's method
freqs, psd = signal.welch(composite_signal, fs = sampling_frequency)

# normalize
psd = psd/psd.max()

plt.plot(freqs, psd)
plt.xlabel('Frequency (Hz)')
plt.ylabel('Power spectral density (normalized)');

Exercise: Use Welch’s method to calculate the frequency spectrum for the composite signal generated below consisting of a sine wave with frequency 30 Hz, a chirp signal starting at 1 Hz and increasing to 15 Hz, and noise. Plot the result.

Hint: If the spectral peaks are all to the left in the plot, you can limit the x-axis to focus on them using plt.xlim([lower_limit, upper_limit]) (insert numbers for lower_limit and upper_limit).

# provided

# Create composite signal
sampling_frequency = 1000
time_vec = np.linspace(0,1,sampling_frequency) # seconds
frequency = 30 # Hz
sineA = np.sin(2*np.pi*frequency*time_vec)

chirp_signal = signal.chirp(time_vec, f0=1, f1=15, t1=1)

noise = np.random.normal(0,1,len(time_vec))

composite_signal = sineA+chirp_signal+noise

utils.plot_signal(composite_signal, time_vec)

Solution

# solution

freqs, psd = signal.welch(composite_signal, fs = sampling_frequency)

psd = psd/psd.max()

plt.plot(freqs, psd)
plt.xlabel('Frequency (Hz)')
plt.ylabel('Power spectral density (normalized)')
plt.xlim([0,100]);

Exercise: Calculate and plot the spectrum again using Welch’s method on the same signal, but this time, set the parameter nperseg to 50 (nperseg = 50) inside the welch function. Can you still tell the peaks of the two original signals that make up the composite signal apart?

Solution

freqs, psd = signal.welch(composite_signal, fs = sampling_frequency, nperseg = 50)

psd = psd/psd.max()

plt.plot(freqs, psd)
plt.xlabel('Frequency (Hz)')
plt.ylabel('Power spectral density (normalized)')
plt.xlim([0,100]);

Exercise: Change the nperseg parameter to 500 and calculate and plot the spectrum again.

Solution

freqs, psd = signal.welch(composite_signal, fs = sampling_frequency, nperseg = 500)

psd = psd/psd.max()

plt.plot(freqs, psd)
plt.xlabel('Frequency (Hz)')
plt.ylabel('Power spectral density (normalized)')
plt.xlim([0,100]);

The nperseg parameter sets the size of the segment of the signal that each frequency spectrum is calculated for in Welch’s method. If the signal contains 1000 samples, and the nperseg is 500 (and there’s no overlap between segments¹), you will have 2 segments of length 500 for which the spectrum is calculated, and the final spectrum is the average of the spectra for those 2 segments. If the signals contains 1000 samples and nperseg is 50, then you’ll have 20 segments, and the final spectrum is the average of those 20 individual spectra.

With many small segments, you get the average of many spectra, which reduces the effect of noise. However, by averaging over many spectra, you will decrease your frequency resolution (it’s harder to tell the peaks apart, as we saw in the exercise with nperseg = 50). With fewer, longer segments, the frequency resolution is higher, but the noise will be more prominent in the final spectrum. Which value to choose for nperseg will depend on the signal-to-noise ratio of your signal and how important it is to have a high frequency resolution.

¹ With Welch’s method, there is some overlap between segments. That was just not taken into account in this example to keep it simpler.

Exercise: The overlap between segments in Welch’s method can also be adjusted with the noverlap parameter. Keep nperseg = 500, but set noverlap to 20 and recalculate and plot the resulting spectrum.

Solution

freqs, psd = signal.welch(composite_signal, sampling_frequency, nperseg=500, noverlap=20)

psd = psd/psd.max()

plt.plot(freqs, psd)
plt.xlabel('Frequency (Hz)')
plt.ylabel('Power spectral density (normalized)')
plt.xlim([0,100]);

Exercise: Keep nperseg = 500, but set noverlap to 400 and recalculate and plot the resulting spectrum. What happens to the first (leftmost) peaks in the spectrum? Can you explain why?

Solution

freqs, psd = signal.welch(composite_signal, sampling_frequency, nperseg=500, noverlap=400)

psd = psd/psd.max()

plt.plot(freqs, psd)
plt.xlabel('Frequency (Hz)')
plt.ylabel('Power spectral density (normalized)')
plt.xlim([0,100]);

Section 3: Multitapering

Similar to Welch’s method, multitapering is also based on averaging across multiple spectra, though the method for calculating the individual spectra is slightly different. With multitapering, your signal is multiplied with different “tapers”, as they are called (ilustrated in panel B in the figure below). The Fourier transform is calculated for each of the differently tapered signals to get a collection of power spectra. In the end, these individual power spectra are averaged to produce the total power spectrum (panel E).

From: Analyzing Neural Time Series, Chapter 16: Multitapers, Mike Cohen.

Exercises

from nitime.algorithms import multi_taper_psd

Example: Use the multitapering method to calculate the power spectrum for the signal generated below. Plot the power spectrum against frequencies.

# provided

# Generate signal
sampling_frequency = 200 # sampling frequency (Hz)
duration = 1 # seconds
time_vec = np.linspace(0,duration,sampling_frequency) # seconds
frequency = 20 #
pure_sine = np.sin(2*np.pi*frequency*time_vec)

# generate noise
noise = np.random.normal(0,1,len(time_vec))

# add noise to sine wave
sine = pure_sine + noise

utils.plot_signal(sine, time_vec)

freqs, psd_mt, _ = multi_taper_psd(sine, Fs=sampling_frequency)

# normalize
psd_mt = psd_mt/psd_mt.max()

# Plot the result
plt.figure()
plt.plot(freqs, psd_mt)
plt.xlabel('Frequency (Hz)')
plt.ylabel('Power Spectral Density (normalized)')
plt.title('Multitaper PSD');

Exercise: Use the multitapering method to calculate the spectrum for the composite signal below consisting of two sine waves with frequencies 10 Hz and 20 Hz. Plot the result.

# provided 

# Create composite signal
sampling_frequency = 200
time_vec = np.linspace(0,1,sampling_frequency) # seconds
frequencyA = 10 # Hz
sineA = np.sin(2*np.pi*frequencyA*time_vec)

frequencyB = 20 # Hz
sineB = np.sin(2*np.pi*frequencyB*time_vec)

composite_signal = sineA+sineB

utils.plot_signal(composite_signal, time_vec)

Solution

freqs, psd_mt, _ = multi_taper_psd(composite_signal, sampling_frequency)

plt.plot(freqs, psd_mt)

plt.xlabel('Frequency (Hz)')
plt.ylabel('Power Spectral Density (normalized)')
plt.title('Multitaper PSD');

Exercise: Use the multitapering method to calculate the power spectrum for the composite signal generated below consisting of a sine wave with frequency 30 Hz, a chirp signal, and noise. Plot the power spectrum against frequencies.

Hint: If the spectral peaks are all to the left in the plot, you can limit the x-axis to focus on them using plt.xlim([lower_limit, upper_limit]) (insert numbers for lower_limit and upper_limit).

# provided

# Create composite signal
sampling_frequency = 1000
time_vec = np.linspace(0,1,sampling_frequency) # seconds
frequency = 30 # Hz
sineA = np.sin(2*np.pi*frequency*time_vec)

chirp_signal = signal.chirp(time_vec, f0=1, f1=15, t1=1)

np.random.seed(5)
noise = np.random.normal(0,1,len(time_vec))

composite_signal = sineA+chirp_signal+noise

utils.plot_signal(composite_signal, time_vec)

Solution

freqs, psd_mt, _ = multi_taper_psd(composite_signal, sampling_frequency)

plt.plot(freqs, psd_mt)

plt.xlabel('Frequency (Hz)')
plt.ylabel('Power Spectral Density (normalized)')
plt.title('Multitaper PSD')
plt.xlim([0,100]);

Section 4: Comparing frequency spectra calculated with different methods.

As you may have noticed in the previous sections, the different methods for computing power spectra tend to produce similar, but not identical results. How the power spectrum looks depends not only on the method used, but also on the particular parameters you used when calculating it. To make sure that your results aren’t just an artifact of the particular method and parameters you chose, it can be a good idea to double check with other methods.

Exercise: Use all three methods - fourier transform, the multitapering method, and Welch’s method - to calculate three different frequency spectra for the composite signal below. Make separate figures and plot and compare the results.

What do you notice? What does the differences between the plots suggest with respect to how well each method deals with noise? How good is the frequency resolution for each method?

# Create composite signal
sampling_frequency = 1000
time_vec = np.linspace(0,1,sampling_frequency) # seconds
frequency = 30 # Hz
sineA = np.sin(2*np.pi*frequency*time_vec)

chirp_signal = signal.chirp(time_vec, f0=1, f1=15, t1=1)

noise = np.random.normal(0,1,len(time_vec))

composite_signal = sineA+chirp_signal+noise

utils.plot_signal(composite_signal, time_vec)

Solution

# solution multitapering

freqs_mt, psd_mt, _ = multi_taper_psd(composite_signal, sampling_frequency)

# normalize
psd_mt = psd_mt/psd_mt.max()

# Plot the result
plt.figure()
plt.plot(freqs_mt, psd_mt)
plt.xlabel('Frequency (Hz)')
plt.ylabel('Power Spectral Density (normalized)')
plt.title('Multitaper PSD')
plt.xlim([0,100]);

Solution

# solution Welch's method

# calculate power spectral density with Welch's method
freqs_welch, psd_welch = signal.welch(composite_signal, fs = sampling_frequency)

# normalize
psd_welch = psd_welch/psd_welch.max()

plt.plot(freqs_welch, psd_welch)
plt.xlabel('Frequency (Hz)')
plt.ylabel('Power spectral density (normalized)')
plt.title('Welch\'s method')
plt.xlim([0,100]);

Solution

# solution fourier transform

freqs_ft, psd_ft = utils.get_fft(composite_signal, time_vec)

plt.plot(freqs_ft,psd_ft)
plt.xlabel('Frequency (Hz)')
plt.ylabel('Power Spectral Density (normalized)')
plt.title('Fourier Transform')
plt.xlim([0,100]);

Exercise: Make figure with subplots with the spectrum from each method in a subplot.

Solution

fig, ax = plt.subplots(ncols = 3, figsize = (15,5), sharex=True)

ax[0].plot(freqs_ft, psd_ft)
ax[0].set_title('PSD Fourier transform')
ax[0].set_xlabel('Frequency (Hz)')

ax[1].plot(freqs_welch, psd_welch)
ax[1].set_title('PSD Welch')

ax[2].plot(freqs_mt, psd_mt)
ax[2].set_title('PSD multitapering')

ax[2].set_xlim([0,100]);