# 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.

--------
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.

--------
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()