Fast Fourier Transform and Power Spectral Density

Fast Fourier Transform (FFT)

The most fundamental procedure to analyse a signal in the frequency domain is the Fourier transform, which can be computed using the algorithm Fast Fourier Transform (FFT).

Given x, a discrete data set (vector) with length N, its forward Discrete Fourier Transform (DFT) is another vector X, also with length N with elements:

$$ X[k] = \sum_{n=0}^{N-1} x[n] \mathrm{e}^{-j2\pi kn/N} \;\;\;\;\;,\;\;\;\;\; 0 \leq k \leq N-1 $$

The Inverse Discrete Fourier Transform (IDFT) inverts this operation and gives back the original vector x:

$$ x[n] = \frac{1}{N} \sum_{k=0}^{N-1} X[k] \mathrm{e}^{j2\pi kn/N} \;\;\;\;\;,\;\;\;\;\; 0 \leq n \leq N-1 $$

The relationship between the DFT and the Fourier coefficients a and b in

$$ x[n] = a_0 + \sum_{k=0}^{N/2-1} a[k]cos\left(\frac{2\pi kt[n]}{Ndt}\right)+b[k]sin\left(\frac{2\pi kt[n]}{Ndt}\right) \;\;\;\;\;,\;\;\;\;\; 0 \leq n \leq N-1 $$

is:

$$ \begin{array}{l} a_0 = X[0]/N \\ \\ a[k] = 2 \cdot real(X[k+1])/N \\ \\ b[k] = -imag(X[k+1])/N \end{array} $$

where x is a length N discrete signal sampled at times t with spacing dt.

In [1]:
from __future__ import division, print_function  # for version compatibility
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
import sys
sys.path.insert(1, r'./../functions')
In [2]:
A = 2 # amplitude
freq = 2 # Hz
phase = np.pi/4 # radians (45 degrees)
fs = 100
dt = 1 / fs
time = np.arange(0, 400) / fs
x = A * np.sin(2 * np.pi * freq * time + phase) + 1
x = np.asarray(x)
N = x.shape[0]
t = np.arange(0, N) / fs

plt.figure(figsize=(8, 3))
plt.plot(time, x, linewidth=2)
plt.xlabel('Time [s]')
plt.ylabel('Amplitude')
plt.tight_layout()
In [3]:
import scipy.fftpack

X = scipy.fftpack.fft(x, N)

# Fourier coefficients (only the positive-frequency terms)
a0 = np.real(X[0]) / N
a = +np.real(X) / N
b = -np.imag(X) / N

# Fourier synthesis based on the Fourier coefficients
y = np.zeros((N, N))
for k in np.arange(0, N):
    w = 2 * np.pi * k / (N * dt)
    y[:, k] = a[k] * np.cos(w * t) + b[k] * np.sin(w * t)
In [4]:
plt.figure(figsize=(8, 3))
plt.plot(time, x, linewidth=2)
plt.plot(t, np.sum(y, axis=1), 'r--', linewidth=2)
plt.xlabel('Time [s]')
plt.ylabel('Amplitude')
plt.tight_layout()
In [5]:
freq = 100.0
t = np.arange(0, 5, .01)
y = 2*np.sin(5*np.pi*2*t) + np.sin(2*np.pi*20*t) + np.random.randn(t.size) 

fig, ax = plt.subplots(1, 1, squeeze=True, figsize=(8, 3))

ax.set_title('Temporal domain', fontsize=18)
ax.plot(t, y, 'b', linewidth=2)
ax.set_xlabel('Time [s]')
ax.set_ylabel('y')
ax.locator_params(axis = 'both', nbins = 5)
plt.tight_layout()
In [6]:
# frequency content
N = y.size
yfft  = scipy.fftpack.fft(y, N)
# let's take only the positive frequencies and normalize the amplitude
yfft  = np.abs(yfft)/N
freqs = scipy.fftpack.fftfreq(N, 1./freq)
freqs = freqs[:np.floor(N/2)]
yfft  = yfft[:np.floor(N/2)]
In [7]:
fig, ax = plt.subplots(1, 1, squeeze=True, figsize=(8, 3))
ax.set_title('Frequency domain', fontsize=18)
ax.plot(freqs, yfft, 'r',  linewidth=2)
ax.set_xlabel('Frequency [Hz]')
ax.set_ylabel('FFT(y)')
ax.locator_params(axis = 'both', nbins = 5)
plt.tight_layout()

Power Spectral Density (PSD)

The function psd.py (code at the end of this text) estimates power spectral density characteristcs using Welch's method. This function is just a wrap of the scipy.signal.welch function with estimation of some frequency characteristcs and a plot. The psd.py returns power spectral density data, frequency percentiles of the power spectral density (for example, Fpcntile[50] gives the median power frequency in Hz); mean power frequency; maximum power frequency; total power, and plots power spectral density data.

Let's exemplify the use of psd.py.

In [8]:
from psd import psd
help(psd)
Help on function psd in module psd:

psd(x, fs=1.0, window='hanning', nperseg=None, noverlap=None, nfft=None, detrend='constant', show=True, ax=None, scales='linear', xlim=None, units='V')
    Estimate power spectral density characteristcs using Welch's method.
    
    This function is just a wrap of the scipy.signal.welch function with
    estimation of some frequency characteristcs and a plot. For completeness,
    most of the help from scipy.signal.welch function is pasted here.
    
    Welch's method [1]_ computes an estimate of the power spectral density
    by dividing the data into overlapping segments, computing a modified
    periodogram for each segment and averaging the periodograms.
    
    Parameters
    ----------
    x : array_like
        Time series of measurement values
    fs : float, optional
        Sampling frequency of the `x` time series in units of Hz. Defaults
        to 1.0.
    window : str or tuple or array_like, optional
        Desired window to use. See `get_window` for a list of windows and
        required parameters. If `window` is array_like it will be used
        directly as the window and its length will be used for nperseg.
        Defaults to 'hanning'.
    nperseg : int, optional
        Length of each segment.  Defaults to half of `x` length.
    noverlap: int, optional
        Number of points to overlap between segments. If None,
        ``noverlap = nperseg / 2``.  Defaults to None.
    nfft : int, optional
        Length of the FFT used, if a zero padded FFT is desired.  If None,
        the FFT length is `nperseg`. Defaults to None.
    detrend : str or function, optional
        Specifies how to detrend each segment. If `detrend` is a string,
        it is passed as the ``type`` argument to `detrend`. If it is a
        function, it takes a segment and returns a detrended segment.
        Defaults to 'constant'.
    show : bool, optional (default = False)
        True (1) plots data in a matplotlib figure.
        False (0) to not plot.
    ax : a matplotlib.axes.Axes instance (default = None)
    scales : str, optional
        Specifies the type of scale for the plot; default is 'linear' which
        makes a plot with linear scaling on both the x and y axis.
        Use 'semilogy' to plot with log scaling only on the y axis, 'semilogx'
        to plot with log scaling only on the x axis, and 'loglog' to plot with
        log scaling on both the x and y axis.
    xlim : float, optional
        Specifies the limit for the `x` axis; use as [xmin, xmax].
        The defaukt is `None` which sets xlim to [0, Fniquist].
    units : str, optional
        Specifies the units of `x`; default is 'V'.
    
    Returns
    -------
    Fpcntile : 1D array
        frequency percentiles of the power spectral density
        For example, Fpcntile[50] gives the median power frequency in Hz.
    mpf : float
        Mean power frequency in Hz.
    fmax : float
        Maximum power frequency in Hz.
    Ptotal : float
        Total power in `units` squared.
    f : 1D array
        Array of sample frequencies in Hz.
    P : 1D array
        Power spectral density or power spectrum of x.
    
    See Also
    --------
    scipy.signal.welch
    
    Notes
    -----
    An appropriate amount of overlap will depend on the choice of window
    and on your requirements.  For the default 'hanning' window an
    overlap of 50% is a reasonable trade off between accurately estimating
    the signal power, while not over counting any of the data.  Narrower
    windows may require a larger overlap.
    If `noverlap` is 0, this method is equivalent to Bartlett's method [2]_.
    
    References
    ----------
    .. [1] P. Welch, "The use of the fast Fourier transform for the
           estimation of power spectra: A method based on time averaging
           over short, modified periodograms", IEEE Trans. Audio
           Electroacoust. vol. 15, pp. 70-73, 1967.
    .. [2] M.S. Bartlett, "Periodogram Analysis and Continuous Spectra",
           Biometrika, vol. 37, pp. 1-16, 1950.
    
    Examples (also from scipy.signal.welch)
    --------
    >>> import numpy as np
    >>> from psd import psd
    #Generate a test signal, a 2 Vrms sine wave at 1234 Hz, corrupted by
    # 0.001 V**2/Hz of white noise sampled at 10 kHz and calculate the PSD:
    >>> fs = 10e3
    >>> N = 1e5
    >>> amp = 2*np.sqrt(2)
    >>> freq = 1234.0
    >>> noise_power = 0.001 * fs / 2
    >>> time = np.arange(N) / fs
    >>> x = amp*np.sin(2*np.pi*freq*time)
    >>> x += np.random.normal(scale=np.sqrt(noise_power), size=time.shape)
    >>> psd(x, fs=freq);

In [9]:
#Generate a test signal, a 2 Vrms sine wave at 1234 Hz, corrupted by
# 0.001 V**2/Hz of white noise sampled at 10 kHz and calculate the PSD:
fs = 10e3
N = 1e5
amp = 2*np.sqrt(2)
freq = 1234.0
noise_power = 0.001 * fs / 2
time = np.arange(N) / fs
x = amp*np.sin(2*np.pi*freq*time)
x += np.random.normal(scale=np.sqrt(noise_power), size=time.shape)

plt.figure(figsize=(8, 3))
plt.plot(time, x, linewidth=2)
plt.xlabel('Time [s]')
plt.ylabel('Amplitude')
plt.tight_layout()
In [10]:
fpcntile, mpf, fmax, Ptotal, f, P = psd(x, fs=freq)

Periodogram

In [11]:
freq = 100.0
t = np.arange(0, 5, .01)
y = 2*np.sin(5*np.pi*2*t) + np.sin(2*np.pi*20*t) + np.random.randn(t.size) 
N = y.shape[0]

from scipy import signal, integrate
fp, Pp = signal.periodogram(y, freq, window='boxcar', nfft=N)
fw, Pw = signal.welch(y, freq, window='hanning', nperseg=N, noverlap=0, nfft=N)
# quick and simple PSD
P  = np.abs(scipy.fftpack.fft(y-np.mean(y),N))[:np.floor(N/2)]**2/N/freq; P[1:-1]=2*P[1:-1]
fs = np.linspace(0,freq/2,len(P))
In [12]:
fig, (ax1,ax2,ax3) = plt.subplots(3, 1, squeeze=True, figsize=(8, 5))
ax1.set_title('Temporal domain', fontsize=18)
ax1.plot(t, y, 'b', linewidth=2)
ax1.set_xlabel('Time [s]')
ax1.set_ylabel('y [V]')
ax1.locator_params(axis = 'both', nbins = 5)
ax2.set_title('Frequency domain', fontsize=18)
ax2.plot(fp, Pp,'r',  linewidth=2)
ax2.set_xlabel('Frequency [Hz]')
ax2.set_ylabel('PSD(y) $[V^2/Hz]$')
ax2.locator_params(axis = 'both', nbins = 5)
ax3.set_title('Frequency domain', fontsize=18)
ax3.plot(fw, Pw,'r',  linewidth=2)
ax3.set_xlabel('Frequency [Hz]')
ax3.set_ylabel('PSD(y) $[V^2/Hz]$')
ax3.locator_params(axis = 'both', nbins = 5)
plt.tight_layout()
In [13]:
F, P = signal.welch(y, fs=freq, window='hanning', nperseg=N/2, noverlap=N/4, nfft=N/2)
A = integrate.cumtrapz(P, F)
fm = np.trapz(F * P, F)/np.trapz(P, F)
f50 = F[np.nonzero(A >= 0.5*A[-1])[0][0]]
f95 = F[np.nonzero(A >= .95*A[-1])[0][0]]
fmax = F[np.argmax(P)]
$$ F_{mean} = \frac{ \sum_{i=1}^{N} F_i*P_i }{ \sum_{i=1}^{N} P_i } $$
In [14]:
from psd import psd
fp, mf, fmax, Ptot, F, P = psd(y, fs=freq, scales='linear', units='V')

Short-Time Fourier Transform

The Short-Time Fourier Transform is a procedure for time-frequency analysis where both frequency and temporal information are estimated. Here are examples of the Short-Time Fourier Transform in action.

In [18]:
fig, ax1 = plt.subplots(1, 1, figsize=(8, 4))
P, freqs, t, im = plt.specgram(y, NFFT=64, Fs=freq, noverlap = 32, cmap=plt.cm.gist_heat)
# P: array of shape (len(times), len(freqs)) of power,
# freqs: array of frequencies, 
# bins: time points the spectrogram is calculated over,
# im: matplotlib.image.AxesImage instance
ax1.set_title('Short-Time Fourier Transform ', fontsize=18)
ax1.set_xlabel('Time [s]')
ax1.set_ylabel('Frequency [Hz]')
plt.tight_layout()
In [19]:
import scipy.signal
t = np.linspace(0, (2**12-1)/1000, 2**12)
c = scipy.signal.chirp(t, f0=100, f1=300, t1=t[-1], method='linear')
fig, ax1 = plt.subplots(1, 1, figsize=(8, 4))
P, freqs, t, im = plt.specgram(c, NFFT=256, Fs=t.size/t[-1], noverlap = 128, cmap=plt.cm.gist_heat)
ax1.set_title('Short-Time Fourier Transform', fontsize=18)
ax1.set_xlabel('Time [s]')
ax1.set_ylabel('Frequency [Hz]')
plt.tight_layout()

Function psd.py

In [17]:
%loadpy ./../functions/psd.py
In [ ]:
#!/usr/bin/env python

"""Estimate power spectral density characteristcs using Welch's method."""

from __future__ import division, print_function
import numpy as np

__author__ = 'Marcos Duarte, https://github.com/demotu/BMC'
__version__ = 'tnorm.py v.1 2013/09/16'


def psd(x, fs=1.0, window='hanning', nperseg=None, noverlap=None, nfft=None,
        detrend='constant', show=True, ax=None, scales='linear', xlim=None,
        units='V'):
    """Estimate power spectral density characteristcs using Welch's method.

    This function is just a wrap of the scipy.signal.welch function with
    estimation of some frequency characteristcs and a plot. For completeness,
    most of the help from scipy.signal.welch function is pasted here.

    Welch's method [1]_ computes an estimate of the power spectral density
    by dividing the data into overlapping segments, computing a modified
    periodogram for each segment and averaging the periodograms.

    Parameters
    ----------
    x : array_like
        Time series of measurement values
    fs : float, optional
        Sampling frequency of the `x` time series in units of Hz. Defaults
        to 1.0.
    window : str or tuple or array_like, optional
        Desired window to use. See `get_window` for a list of windows and
        required parameters. If `window` is array_like it will be used
        directly as the window and its length will be used for nperseg.
        Defaults to 'hanning'.
    nperseg : int, optional
        Length of each segment.  Defaults to half of `x` length.
    noverlap: int, optional
        Number of points to overlap between segments. If None,
        ``noverlap = nperseg / 2``.  Defaults to None.
    nfft : int, optional
        Length of the FFT used, if a zero padded FFT is desired.  If None,
        the FFT length is `nperseg`. Defaults to None.
    detrend : str or function, optional
        Specifies how to detrend each segment. If `detrend` is a string,
        it is passed as the ``type`` argument to `detrend`. If it is a
        function, it takes a segment and returns a detrended segment.
        Defaults to 'constant'.
    show : bool, optional (default = False)
        True (1) plots data in a matplotlib figure.
        False (0) to not plot.
    ax : a matplotlib.axes.Axes instance (default = None)
    scales : str, optional
        Specifies the type of scale for the plot; default is 'linear' which
        makes a plot with linear scaling on both the x and y axis.
        Use 'semilogy' to plot with log scaling only on the y axis, 'semilogx'
        to plot with log scaling only on the x axis, and 'loglog' to plot with
        log scaling on both the x and y axis.
    xlim : float, optional
        Specifies the limit for the `x` axis; use as [xmin, xmax].
        The defaukt is `None` which sets xlim to [0, Fniquist].
    units : str, optional
        Specifies the units of `x`; default is 'V'.

    Returns
    -------
    Fpcntile : 1D array
        frequency percentiles of the power spectral density
        For example, Fpcntile[50] gives the median power frequency in Hz.
    mpf : float
        Mean power frequency in Hz.
    fmax : float
        Maximum power frequency in Hz.
    Ptotal : float
        Total power in `units` squared.
    f : 1D array
        Array of sample frequencies in Hz.
    P : 1D array
        Power spectral density or power spectrum of x.

    See Also
    --------
    scipy.signal.welch

    Notes
    -----
    An appropriate amount of overlap will depend on the choice of window
    and on your requirements.  For the default 'hanning' window an
    overlap of 50% is a reasonable trade off between accurately estimating
    the signal power, while not over counting any of the data.  Narrower
    windows may require a larger overlap.
    If `noverlap` is 0, this method is equivalent to Bartlett's method [2]_.

    References
    ----------
    .. [1] P. Welch, "The use of the fast Fourier transform for the
           estimation of power spectra: A method based on time averaging
           over short, modified periodograms", IEEE Trans. Audio
           Electroacoust. vol. 15, pp. 70-73, 1967.
    .. [2] M.S. Bartlett, "Periodogram Analysis and Continuous Spectra",
           Biometrika, vol. 37, pp. 1-16, 1950.

    Examples (also from scipy.signal.welch)
    --------
    >>> import numpy as np
    >>> from psd import psd
    #Generate a test signal, a 2 Vrms sine wave at 1234 Hz, corrupted by
    # 0.001 V**2/Hz of white noise sampled at 10 kHz and calculate the PSD:
    >>> fs = 10e3
    >>> N = 1e5
    >>> amp = 2*np.sqrt(2)
    >>> freq = 1234.0
    >>> noise_power = 0.001 * fs / 2
    >>> time = np.arange(N) / fs
    >>> x = amp*np.sin(2*np.pi*freq*time)
    >>> x += np.random.normal(scale=np.sqrt(noise_power), size=time.shape)
    >>> psd(x, fs=freq);
    """

    from scipy import signal, integrate

    if not nperseg:
        nperseg = np.ceil(len(x) / 2)
    f, P = signal.welch(x, fs, window, nperseg, noverlap, nfft, detrend)
    Area = integrate.cumtrapz(P, f, initial=0)
    Ptotal = Area[-1]
    mpf = integrate.trapz(f * P, f) / Ptotal  # mean power frequency
    fmax = f[np.argmax(P)]
    # frequency percentiles
    inds = [0]
    Area = 100 * Area / Ptotal  # + 10 * np.finfo(np.float).eps
    for i in range(1, 101):
        inds.append(np.argmax(Area[inds[-1]:] >= i) + inds[-1])
    fpcntile = f[inds]

    if show:
        _plot(x, fs, f, P, mpf, fmax, fpcntile, scales, xlim, units, ax)

    return fpcntile, mpf, fmax, Ptotal, f, P


def _plot(x, fs, f, P, mpf, fmax, fpcntile, scales, xlim, units, ax):
    """Plot results of the ellipse function, see its help."""
    try:
        import matplotlib.pyplot as plt
    except ImportError:
        print('matplotlib is not available.')
    else:
        if ax is None:
            fig, ax = plt.subplots(1, 1, figsize=(7, 5))
        if scales.lower() == 'semilogy' or scales.lower() == 'loglog':
            ax.set_yscale('log')
        if scales.lower() == 'semilogx' or scales.lower() == 'loglog':
            ax.set_xscale('log')
        plt.plot(f, P, linewidth=2)
        ylim = ax.get_ylim()
        plt.plot([fmax, fmax], [np.max(P), np.max(P)], 'ro',
                 label='Fpeak  = %.2f' % fmax)
        plt.plot([fpcntile[50], fpcntile[50]], ylim, 'r', lw=1.5,
                 label='F50%%   = %.2f' % fpcntile[50])
        plt.plot([mpf, mpf], ylim, 'r--', lw=1.5,
                 label='Fmean = %.2f' % mpf)
        plt.plot([fpcntile[95], fpcntile[95]], ylim, 'r-.', lw=2,
                 label='F95%%   = %.2f' % fpcntile[95])
        leg = ax.legend(loc='best', numpoints=1, framealpha=.5,
                        title='Frequencies [Hz]')
        plt.setp(leg.get_title(), fontsize=12)
        plt.xlabel('Frequency [$Hz$]', fontsize=12)
        plt.ylabel('Magnitude [%s$^2/Hz$]' % units, fontsize=12)
        plt.title('Power spectral density', fontsize=12)
        if xlim:
            ax.set_xlim(xlim)
        ax.set_ylim(ylim)
        plt.tight_layout()
        plt.grid()
        plt.show()