{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "

An Interactive Fast Fourier Transform (FFT)

\n", "

MCHE 485: Mechanical Vibrations

\n", "

Dr. Joshua Vaughan
\n", "joshua.vaughan@louisiana.edu
\n", "http://www.ucs.louisiana.edu/~jev9637/

" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This notebook will present an interactive [Fast Fourier Tranform](https://en.wikipedia.org/wiki/Fast_Fourier_transform), or FFT for a simple, term-term harmonic signal. For more infomation on the FFT, see the [Introduction to the FFT Notebook](http://nbviewer.jupyter.org/github/DocVaughan/MCHE485---Mechanical-Vibrations/blob/Spring2016/Jupyter%20Notebooks/FFT%20-%20Introduction.ipynb).\n", "\n", "The signal that this notebook analyzing is a simple sum of three sine waves:\n", "\n", "$ \\quad f(t) = A_1 \\sin{(\\omega_1 t)} + A_2 \\sin{(\\omega_2 t)} + A_3 \\sin{(\\omega_3 t)} $ \n", "\n", "The sliders presented in the interactive portion of the notebook allow amplitudes, $A_1$, $A_2$, and $A_3$, and frequencies, $\\omega_1$, $\\omega_2$, and $\\omega_3$, to be changed. Any time any of these parameters is changed, the plots will automatically regenerate.\n", "\n", "The \"best\" way to use this notebook is to use the Cell... Run All command from the menu bar, then scroll down to the bottom of the notebook and play with the sliders.\n", "\n", "*NOTE*: The interactive portion of this notebook will not run in the online notebook viewer, so you'll have to run it locally to play with the interactive part." ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# Grab all of the NumPy functions with namespace (nickname) np\n", "import numpy as np " ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# Set up the notebook to display plots inline\n", "%matplotlib inline\n", "\n", "# Import the plotting functions \n", "import matplotlib.pyplot as plt" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# import the IPython widgets\n", "from ipywidgets.widgets import interact\n", "from ipywidgets import widgets # Widget definitions\n", "from IPython.display import display # Used to display widgets in the notebook" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "collapsed": true }, "outputs": [], "source": [ "t = np.linspace(0, 25, 5001) # Time, 0-25s with 5001 samples in the range" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "collapsed": true }, "outputs": [], "source": [ "def plot_FFT(A1 = 1.0, f1 = 1.0, A2 = 0.5, f2 = 2.75, A3 = 1.5, f3 = 7.5):\n", " \n", " w1 = f1 * (2 * np.pi) # frequency of first sine (rad/s)\n", " w2 = f2 * (2 * np.pi) # frequency of second sine (rad/s)\n", " w3 = f3 * (2 * np.pi) # frequency of third sine (rad/s)\n", "\n", " orig_signal = A1 * np.sin(w1 * t) + A2 * np.sin(w2 * t) + A3 * np.sin(w3 * t)\n", " \n", " # We can also use the FFT to get the natrual frequency\n", " freq, mag = CRAWLAB_fft(orig_signal, t, False)\n", " \n", " # Make the figure pretty, then plot the results\n", " # \"pretty\" parameters selected based on pdf output, not screen output\n", " # Many of these setting could also be made default by the .matplotlibrc file\n", " fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12,4))\n", "\n", " plt.subplots_adjust(bottom=0.12,left=0.17,top=0.96,right=0.96)\n", " plt.setp(ax1.get_ymajorticklabels(),family='serif',fontsize=18)\n", " plt.setp(ax1.get_xmajorticklabels(),family='serif',fontsize=18)\n", "\n", " ax1.spines['right'].set_color('none')\n", " ax1.spines['top'].set_color('none')\n", " ax1.xaxis.set_ticks_position('bottom')\n", " ax1.yaxis.set_ticks_position('left')\n", " ax1.grid(True,linestyle=':',color='0.75')\n", " ax1.set_axisbelow(True)\n", "\n", " ax2.spines['right'].set_color('none')\n", " ax2.spines['top'].set_color('none')\n", " ax2.xaxis.set_ticks_position('bottom')\n", " ax2.yaxis.set_ticks_position('left')\n", " ax2.grid(True,linestyle=':',color='0.75')\n", " ax2.set_axisbelow(True)\n", "\n", " # Original Signal\n", " ax1.set_xlabel('Time (s)', family='serif', fontsize=22, weight='bold', labelpad=5)\n", " ax1.set_ylabel('Signal Amplitude', family='serif', fontsize=22, weight='bold', labelpad=10)\n", "\n", " ax1.plot(t, orig_signal, linewidth=2, linestyle='-', label=r'Original Signal')\n", " ax1.set_xlim(0, 5)\n", "\n", " # FFT\n", " ax2.set_xlabel('Frequency (Hz)',fontsize=22, labelpad=5)\n", " ax2.set_ylabel('FFT magnitude',fontsize=22, labelpad=10)\n", "\n", " ax2.plot(freq, mag, linewidth=2, linestyle='-', label = 'FFT')\n", "\n", " # Let's also annotate the plot with vertical lines at each of the frequencies in our original signal\n", " ax2.plot([w1 / (2 * np.pi), w1 / (2 * np.pi)], [0, 1], linewidth = 2, linestyle = '--', label = r'$\\omega_1$')\n", " ax2.plot([w2 / (2 * np.pi), w2 / (2 * np.pi)], [0, 1], linewidth = 2, linestyle = '-.', label = r'$\\omega_2$')\n", " ax2.plot([w3 / (2 * np.pi), w3 / (2 * np.pi)], [0, 1], linewidth = 2, linestyle = ':', label = r'$\\omega_3$')\n", "\n", " ax2.set_xlim(0, 10)\n", " ax2.set_ylim(0, 1.0)\n", "\n", " ax2.leg = ax2.legend(loc='upper right', ncol = 2, fancybox=True)\n", " ltext = ax2.leg.get_texts()\n", " plt.setp(ltext,family='Serif',fontsize=16)\n", "\n", " # Adjust the page layout filling the page using the new tight_layout command\n", " plt.tight_layout(pad=0.5, w_pad=5.0)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "I've written a function for my lab that takes care of a lot of the boilerplate code necessary to complete and plot the FFT. The function, defined in the next cell, is named ```CRAWLAB_fft```. It relies on the ```fft``` function from the ```fftpack``` module of SciPy.\n", "\n", "The function recenters data about zero, applies a Hanning window to the data, and selects the real-valued components returned from the SciPy ```fft``` function. If ```plotflag``` is ```True```, then it will also automatically generate the plot of the FFT magnitude. Reasons for doing these things are a bit beyond this tutorial, so, for now, you can just think of ```CRAWLAB_fft``` as implementing the FFT algorithm for you." ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "collapsed": true }, "outputs": [], "source": [ "def CRAWLAB_fft(data, time, plotflag):\n", " ''' Function to get the FFT for a response\n", " #\n", " # Inputs:\n", " # time = time array corresponding to the data\n", " # data = the response data array (only pass a single dimension/state at at time)\n", " # plotflag = will plot the FFT if nonzero\n", " # \n", " # Output:\n", " # fft_freq = an array of the freqs used in the FFT\n", " # fft_mag = an array of the amplitude of the FFT at each freq in fft_freq\n", " #\n", " # Created: 03/28/14\n", " # - Joshua Vaughan\n", " # - joshua.vaughan@louisiana.edu\n", " # - http://www.ucs.louisiana.edu/~jev9637\n", " ######################################################################################\n", " '''\n", " \n", " from scipy.fftpack import fft\n", " \n", " # correct for any DC offset\n", " offset = np.mean(data) \n", "\n", " # Get the sampling time\n", " sample_time = time[1] - time[0]\n", " \n", " # Get the length of the dataset\n", " n = len(data)\n", "\n", " # Calculate the FFT of the data, removing the offset and using a Hanning Window\n", " fft_mag = fft((data - offset) * np.hanning(len(data)))\n", " \n", " # Define the frequency range of the output\n", " fft_freq = np.linspace(0.0, 1.0 / (2.0*sample_time), int(np.ceil(n/2)))\n", " \n", " # Only return the \"useful\" part of the fft\n", " fft_mag = 2.0/n * np.abs(fft_mag[0:int(np.ceil(n/2))])\n", " \n", " # If plotflag is nonzero (True), plot the FFT before returning the magnitude and phase\n", " if plotflag:\n", " # Plot the relationshiop\n", " # Many of these setting could also be made default by the .matplotlibrc file\n", " fig = plt.figure(figsize=(6,4))\n", " ax = plt.gca()\n", " plt.subplots_adjust(bottom=0.17,left=0.17,top=0.96,right=0.96)\n", " plt.setp(ax.get_ymajorticklabels(),fontsize=18)\n", " plt.setp(ax.get_xmajorticklabels(),fontsize=18)\n", " ax.spines['right'].set_color('none')\n", " ax.spines['top'].set_color('none')\n", " ax.xaxis.set_ticks_position('bottom')\n", " ax.yaxis.set_ticks_position('left')\n", " ax.grid(True, linestyle=':', color='0.75')\n", " ax.set_axisbelow(True)\n", "\n", " plt.xlabel('Frequency (Hz)', fontsize=22, labelpad=8)\n", " plt.ylabel('FFT magnitude', fontsize=22, labelpad=10)\n", " \n", " plt.plot(fft_freq, fft_mag, linewidth=2, linestyle='-')\n", " \n", " # Adjust the page layout filling the page using the new tight_layout command\n", " plt.tight_layout(pad=0.5)\n", " plt.show()\n", " \n", " return fft_freq, fft_mag" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now will make the call to the interact function, setting up the ranges for each parameter. The result will be six sliders, representing the amplitude and frequency of each component of the signal. Each time one of the sliders is changed, the plot of the signal and the FFT should update." ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "0c99abc482414d8c90ecaf5d980acb5d", "version_major": 2, "version_minor": 0 }, "text/html": [ "

Failed to display Jupyter Widget of type interactive.

\n", "

\n", " If you're reading this message in the Jupyter Notebook or JupyterLab Notebook, it may mean\n", " that the widgets JavaScript is still loading. If this message persists, it\n", " likely means that the widgets JavaScript library is either not installed or\n", " not enabled. See the Jupyter\n", " Widgets Documentation for setup instructions.\n", "

\n", "

\n", " If you're reading this message in another frontend (for example, a static\n", " rendering on GitHub or NBViewer),\n", " it may mean that your frontend doesn't currently support widgets.\n", "

\n" ], "text/plain": [ "interactive(children=(FloatSlider(value=1.0, description='A1', max=1.0), FloatSlider(value=1.0, description='f1', max=10.0, step=0.25), FloatSlider(value=0.5, description='A2', max=1.0), FloatSlider(value=2.75, description='f2', max=10.0, step=0.25), FloatSlider(value=1.0, description='A3', max=1.0), FloatSlider(value=7.5, description='f3', max=10.0, step=0.25), Output()), _dom_classes=('widget-interact',))" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Call the slider interaction\n", "interact(plot_FFT, A1 = (0, 1, 0.1), \n", " f1 = (0, 10, 0.25), \n", " A2 = (0, 1, 0.1), \n", " f2 = (0, 10, 0.25), \n", " A3 = (0, 1, 0.1), \n", " f3 = (0, 10, 0.25));" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "If you'd like to explore the FFT further, the tutorials below and the links referenced within them have much more information:\n", "\n", "* [Understanding the FFT Algorithm][0-3]\n", "* [The Math Trick Behind MP3s, JPEGs, and Homer Simpson’s Face][0-4]\n", "* [NI FFT Fundamentals][0-1]\n", "* [SciPy Lectures Notes FFT][0-5]\n", "* [Fourier transform for dummies][0-2]\n", "\n", "[0-1]: http://zone.ni.com/reference/en-XX/help/372416B-01/svtconcepts/fft_funda/ \"FFT Fundamentals (Sound and Vibration Measurement Suite) - Sound and Vibration Measurement Suite 7.0 Help - National Instruments\"\n", "[0-2]: http://nipunbatra.github.io/2016/01/fft/ \"Fourier transform for dummies | Nipun Batra |\"\n", "[0-3]: http://jakevdp.github.io/blog/2013/08/28/understanding-the-fft/ \"Understanding the FFT Algorithm\"\n", "[0-4]: http://nautil.us/blog/the-math-trick-behind-mp3s-jpegs-and-homer-simpsons-face \"The Math Trick Behind MP3s, JPEGs, and Homer Simpson’s Face - Facts So Romantic - Nautilus\"\n", "[0-5]: http://www.scipy-lectures.org/intro/scipy.html#fast-fourier-transforms-scipy-fftpack \"1.5. Scipy : high-level scientific computing — Scipy lecture notes\"" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Licenses\n", "Code is licensed under a 3-clause BSD style license. See the licenses/LICENSE.md file.\n", "\n", "Other content is provided under a [Creative Commons Attribution-NonCommercial 4.0 International License](http://creativecommons.org/licenses/by-nc/4.0/), CC-BY-NC 4.0.\n" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\n" ], "text/plain": [ "" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# This cell will just improve the styling of the notebook\n", "from IPython.core.display import HTML\n", "import urllib.request\n", "response = urllib.request.urlopen(\"https://cl.ly/1B1y452Z1d35\")\n", "HTML(response.read().decode(\"utf-8\"))" ] } ], "metadata": { "anaconda-cloud": {}, "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.6.2" } }, "nbformat": 4, "nbformat_minor": 1 }