.. include:: /project-links.txt
.. include:: /abbreviation.txt



.. getthecode:: fft.py
    :language: python3
    :hidden:


========================
 Fast Fourier Transform
========================

This example shows how to compute a FFT of a signal using the scipy Scientific Python package.


.. code-block:: py3

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


We will first compute the spectrum of the sum of two sinusoidal waveforms.


.. code-block:: py3

    
    N = 1000 # number of sample points
    dt = 1. / 500 # sample spacing
    
    frequency1 = 50.
    frequency2 = 80.
    
    t = np.linspace(0.0, N*dt, N)
    y = np.sin(2*np.pi * frequency1 * t) + .5 * np.sin(2*np.pi * frequency2 * t)
    
    yf = fft(y)
    tf = np.linspace(.0, 1./(2.*dt), N/2)
    spectrum = 2./N * np.abs(yf[0:N//2])
    
    figure1 = plt.figure(1, (20, 10))
    plt.plot(tf, spectrum, 'o-')
    plt.grid()
    for frequency in frequency1, frequency2:
        plt.axvline(x=frequency, color='red')
    plt.title('Spectrum')
    plt.xlabel('Frequency [Hz]')
    plt.ylabel('Amplitude')
    


.. image:: fft-sum-of-sin.png
    :align: center


Now we will compute the spectrum of a square waveform.

The Fourier series is given by:

.. math::

  \frac{4}{\pi} \sum_{n=1, 3, 5, \ldots}^{\inf} \frac{1}{n} \sin(n 2\pi f t)


.. code-block:: py3

    
    N = 1000 # number of sample points
    dt = 1. / 1000 # sample spacing
    
    frequency = 5.
    
    t = np.linspace(.0, N*dt, N)
    y = signal.square(2*np.pi*frequency*t)
    
    figure2 = plt.figure(2, (20, 10))
    
    plt.subplot(211)
    plt.plot(t, y)
    y_sum = None
    for n in range(1, 20, 2):
        yn = 4/(np.pi*n)*np.sin((2*np.pi*n*frequency*t))
        if y_sum is None:
            y_sum = yn
        else:
            y_sum += yn
        if n in (1, 3, 5):
            plt.plot(t, y_sum)
    plt.plot(t, y_sum)
    plt.xlim(0, 2/frequency)
    plt.ylim(-1.5, 1.5)
    
    yf = fft(y)
    tf = np.linspace(.0, 1./(2.*dt), N/2)
    spectrum = 2./N * np.abs(yf[0:N//2])
    
    plt.subplot(212)
    plt.plot(tf, spectrum)
    n = np.arange(1, 20, 2)
    plt.plot(n*frequency, 4/(np.pi*n), 'o', color='red')
    plt.grid()
    plt.title('Spectrum')
    plt.xlabel('Frequency [Hz]')
    plt.ylabel('Amplitude')
    


.. image:: fft-square-waveform.png
    :align: center


.. code-block:: py3

    
    
    plt.show()