Open in molab View on GitHub Download .py

Fourier transforms¶

From very naive approximations of the spectrum to more complex examples:¶

import matplotlib as mpl
import matplotlib.pyplot as plt
import numpy as np
from scipy import fft

Two basic ploting functions that repeat the actuall spectral analysis:¶

# redefine default figure size and fonts

# create signal
def create_signal(fs, N, a=[1., 1.5], f=[10.0, 35.7], DC=0):
    """ create a periodic signal with a DC and a Gaussian noise"""
    dt = 1./fs
    t = np.linspace(0, N*dt, N)
    y = np.random.normal(0, 1, N) + DC
    for aa, ff in zip(a, f):
        y += aa*np.sin(2*np.pi*ff*t)

    noise = np.random.normal(0, 1, N)
    y += noise
    return t, y

def spectrum(y, Fs):
    """
    Plots a Single-Sided Amplitude Spectrum of a sampled
    signal y(t), sampling frequency Fs (length of a signal
    provides the number of samples recorded)

    Following: http://goo.gl/wRoUn
    """
    n = len(y)  # length of the signal
    k = np.arange(n)
    T = n/Fs
    frq = k/T  # two sides frequency range
    frq = frq[range(int(n/2))]  # one side frequency range
    Y = 2*np.fft.fft(y)/n  # np.fft computing and normalization
    Y = Y[range(int(n/2))]
    return (frq, Y)

def plotSignal(A, ff, fs, N):

    T = N/fs  # sampling period
    t = np.arange(0.0, T, T/N)  # sampling time steps
    y = A*np.sin(2*np.pi*ff*t)  # sampled signal
    frq, Y = spectrum(y, fs)  # FFT(sampled signal)

    # Plot
    fig, ax = plt.subplots(2, 1, figsize=(10, 8))
    ax[0].plot(t, y, 'b.')
    ax[0].set_xlabel('$t$ [s]')
    ax[0].set_ylabel('Y [V]')
    ax[0].plot(t[0], y[0], 'ro')
    ax[0].plot(t[-1], y[-1], 'ro')

    ax[1].plot(frq, abs(Y), 'r')  # plotting the spectrum
    ax[1].set_xlabel('$f$ (Hz)')
    ax[1].set_ylabel('$|Y(f)|$')
    lgnd = str(r'N = %d, f = %d Hz' % (N, fs))
    ax[1].legend([lgnd], loc='best')

More elaborate example, demonstration of a leakage effect:¶

# We sample a signal at fs = 200 Hz and record 256 points"
A = 1.0
# true values
ff = 10.0  # Volt, amplitude
_fs = 200.0  # Hz, signal frequency, zero harmonics
_N = 256
# We will work with different sampling frequencies
# and different lengths of records:
plotSignal(A, ff, _fs, _N)  #Hz  # points

Note:¶

leakage at around 10 Hz
peak is below 1 Volt
we can get the spectrum up to half of sampling frequency
frequency resolution is \(\Delta f = 1/T = 0.781\) Hz

Let us try to minimize the leakage:¶

Let's first increase the sampling rate, stay with N = 256¶

_fs = 1000.0  #Hz
_N = 512  # points
plotSignal(A, ff, _fs, _N)

Note:¶

resolution in time is great BUT:
resolution in frequency is worse, \(\Delta f = 1/T = 1/0.25 = 3.91\) Hz
we see up to 500 Hz, but have here only 128 useful points
leakage is severe, the peak amplitude is about 0.66 Volt
peak location is now 11.7 Hz

MORE IS NOT NECESSARILY BETTER

Let us try to minimize the leakage problem sampling closer to the Nyquist frequency¶

_fs = 25.0  #Hz
_N = 256  # points
plotSignal(A, ff, _fs, _N)

Note:¶

sine does not seem to be a sine wave, time resolution is bad for Nyquist frequency
frequency resolution is better, about 0.1 Hz
peak is about 0.75 Volt, less leakage but still strong
peak location is close, 9.96 Hz

Let's get more points, keep same sampling frequency¶

_fs = 25.0  #Hz
_N = 512  # points
plotSignal(A, ff, _fs, _N)

Note:¶

time sampling is much longer, \(T = N/f_s = 512/(25 Hz) = 20.48\) s
peak is narrow, at 10.01 Hz
resolution is good, much lower leakage, value is close to 0.9 Volt

Still, how to get the perfect spectrum? Use tricks, for instance, sample at a specific frequency:¶

_fs = 25.6  #Hz
_N = 256  # points
plotSignal(A, ff, _fs, _N)

# After you know the real value, you can save a lot of time:
_fs = 25.6  #Hz
_N = 64  # points
plotSignal(A, ff, _fs, _N)