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

Interactive Fourier Analysis

\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 look at the Fourier Analysis of periodic functions. The exact function is compared to an N-term approximation. Then, the response to the exact and approximate inputs will be compared.\n", "\n", "We can write any periodic function as an infinite sum of sines and cosines:\n", "\n", "$ \\quad f(t) = \\sum_{n=0}^{\\infty}a_n\\cos(n\\omega_0t) + \\sum_{n=1}^{\\infty}b_n\\sin(n\\omega_0t) $\n", "\n", "where\n", "\n", "$ \\quad a_n = \\frac{\\omega_0}{\\pi}\\int_0^{\\frac{2\\pi}{\\omega_0}}f(t)\\cos(n\\omega_0t)dt $, \n", "\n", "$ \\quad b_n = \\frac{\\omega_0}{\\pi}\\int_0^{\\frac{2\\pi}{\\omega_0}}f(t)\\sin(n\\omega_0t)dt $, \n", "\n", "and\n", "\n", "$ \\quad a_0 = \\frac{\\omega_0}{2\\pi}\\int_0^{\\frac{2\\pi}{\\omega_0}}f(t)dt $\n", "\n", "For more information on this process, you can see the lectures at the [class website](http://www.ucs.louisiana.edu/~jev9637/MCHE485.html). In the code below, we'll solve the required integrals numerically.\n", "\n", "By changing the number of terms we use to approximate the original function we can approach its shape, as shown in Figure 1.\n", "

\n", "\t\"Fourier
\n", " Figure 1: Increasing the Number of Terms in the Approximation \n", "

\n", "\n", "### Interactive Plotting\n", "We'll set up interactive versions of the Fourier Approximation and the response resulting from using this approximation. Sliders will allow us to adjust the number of terms used in the approximation. There may be some flickering during the first time you change the sliders.\n", "\n", "The interactive portion will not run in the online notebook viewer, so you'll have to run it locally to play with the interactive part. The .gif in Figure 2 shows an example of what the interactive plot should look like after you run an interactive cell locally.\n", "

\n", "\t\"Interactive
\n", " Figure 2: Interactively Changing the Number of Terms in the Approximation \n", "



" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": true }, "outputs": [], "source": [ "import numpy as np # Grab all of the NumPy functions with nickname np" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# We want our plots to be displayed inline, not in a separate window\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": [ "#--------- Input your function to examine here --------\n", "\n", "# This is the square wave we did in class\n", "t = np.linspace(0, 4, 4000) # define the time to look at, easiest to just choose 1 period\n", "w0 = 2.0 * np.pi/4.0 # define the fundamental frequency (here, I know t(end)=tau)\n", "tau_0 = 2.0 * np.pi / w0 # define fundamental period based on w0\n", "\n", "y = (t > 0) - 2*(t>2)\n", "\n", "# This is the \"trapezoid\" wave example we worked in class\n", "# t = np.linspace(0,4,4000) # define the time to look at, easiest to just choose 1 period\n", "# w0 = 2*np.pi/t[-1] # define the fundamental frequency (here, I know t(end)=tau)\n", "# tau_0 = 2*np.pi/w0 # define fundamental period based on w0\n", "\n", "# F0 = 1\n", "# y = np.zeros((len(t),))\n", "\n", "# for ii in range(len(t)):\n", "# if t[ii] <= tau_0/3:\n", "# y[ii] = 3*F0/tau_0*t[ii]\n", "# elif t[ii] <= 2*tau_0/3:\n", "# y[ii] = F0\n", "# else:\n", "# y[ii] = -3*F0/tau_0*t[ii]+3*F0\n" ] }, { "cell_type": "code", "execution_count": 4, "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\n", "\n", "# Set up the function that plots the repsonse based on slider changes\n", "def plot_approx(num_terms = 11):\n", "\n", " # get the a0 term\n", " a0 = w0 / (2.0 * np.pi) * np.trapz(y, t) \n", "\n", " # fill arrays with zeros - Good practice, it speeds computation in large problems\n", " a = np.zeros((num_terms,))\n", " b = np.zeros((num_terms,))\n", " integral_cos = np.zeros((len(t),num_terms))\n", " integral_sin = np.zeros((len(t),num_terms))\n", " sin_term = np.zeros((num_terms,len(t)))\n", " cos_term = np.zeros((num_terms,len(t)))\n", "\n", " # cycle through the 1 to num_terms Fourier coefficients (a_n and b_n)\n", " for n in range(num_terms):\n", "\n", " # a_n calculations\n", " integral_cos[:,n] = y * np.cos((n+1) * w0 * t) # define the integral \"interior\"\n", " a[n] = w0/np.pi * np.trapz(integral_cos[:,n], t) # solve for a_n\n", "\n", " # b_n calculations\n", " integral_sin[:,n] = y * np.sin((n+1) * w0 * t) # define the integral \"interior\"\n", " b[n] = w0/np.pi * np.trapz(integral_sin[:,n], t) # solve for b_n\n", "\n", " sin_term[n,:] = np.sin((n+1) * w0 * t) # calculate the nth sine term\n", " cos_term[n,:] = np.cos((n+1) * w0 * t) # calculate the nth cosine_term\n", "\n", "\n", " # Generate the approximate input based on the Fourier coeff. calculated above\n", " approx = np.zeros_like(t) #First fill with zeros\n", "\n", " for ii in range(len(t)):\n", " approx[ii] = a0 + np.sum(a * cos_term[:,ii],0) + np.sum(b * sin_term[:,ii],0)\n", " \n", " # Now, let's plot the comparison\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 = plt.figure(figsize=(9,6))\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(),family='serif',fontsize=18)\n", " plt.setp(ax.get_xmajorticklabels(),family='serif',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(r'Time (s)', family='serif', fontsize=22, weight='bold', labelpad=5)\n", " plt.ylabel(r'$y(t)$', family='serif', fontsize=22, weight='bold', labelpad=10)\n", "\n", " plt.plot(t, y, '--', linewidth=2, label=r'Exact')\n", "\n", " f = str(num_terms) + '-Term Fourier Expansion'\n", " plt.plot(t, approx, linewidth=2, label=f)\n", "\n", " plt.ylim(-1.5,2.5)\n", "\n", " leg = plt.legend(loc='upper right', ncol = 1, fancybox=True)\n", " ltext = leg.get_texts() \n", " plt.setp(ltext,family='Serif',fontsize=18)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Interactive Plot\n", "Now, we'll create the interactive plot. The result will have a slider that allows us to change the number of terms used in the Fourier Approximation.\n", "\n", "***Reminder***: You will have the run the notebook locally to use the interactive features. They are *not* viewable at [http://nbviewer.jupyter.org](http://nbviewer.jupyter.org)." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "c3f3f5b77d3a493aa1606ef12658e001", "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=(IntSlider(value=11, description='num_terms', max=22, min=1), Output()), _dom_classes=('widget-interact',))" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Call the slider interaction\n", "# num_terms changes in the number of terms in the Fourier Approx, allowing between 1 and 21\n", "interact(plot_approx, num_terms=(1, 22, 1));" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Response Comparison\n", "Now, let's look at the response of a simple mass-spring-damper system like the one shown in Figure 3.\n", "

\n", "\t\"A
\n", " Figure 3: A Mass-Spring-Damper System \n", "

\n", "\n", "The equation of motion for the system is:\n", "\n", "$ \\quad m \\ddot{x} + c \\dot{x} + kx = c \\dot{y} + ky $\n", "\n", "We could also write this equation in terms of the damping ratio, $\\zeta$, and natural frequency, $\\omega_n$.\n", "\n", "$ \\quad \\ddot{x} + 2\\zeta\\omega_n\\dot{x} + \\omega_n^2x = 2\\zeta\\omega_n \\dot{y} + \\omega_n^2 y$\n", "\n", "For information on how to obtain this equation, you can see the lectures at the [class website](http://www.ucs.louisiana.edu/~jev9637/MCHE485.html).\n", "\n", "### Control System Library\n", "For the following simulations in this notebook, we will use the [Control Systems Library for Python](http://www.cds.caltech.edu/~murray/wiki/Control_Systems_Library_for_Python). Instructions on installation and use can be found at that link." ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "collapsed": true }, "outputs": [], "source": [ "import control # import the control system library" ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# Define the System Parameters\n", "m = 1.0 # kg\n", "k = (2.0 * np.pi)**2 # N/m (Selected to give an undamped natrual frequency of 1Hz)\n", "wn = np.sqrt(k / m) # Natural Frequency (rad/s)\n", "\n", "z = 0.25 # Define a desired damping ratio\n", "c = 2 * z * wn * m # calculate the damping coeff. to create it (N/(m/s))\n", "\n", "wd = wn * np.sqrt(1 - z**2) # Damped natural frequency (rad/s)" ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# Define the system to use in simulation - in transfer function form here\n", "num = [2.0 * z * wn, wn**2]\n", "den = [1, 2.0 * z * wn, wn**2]\n", "\n", "# Definte the transfer function\n", "sys = control.tf(num, den)\n", "\n", "# run the simulation - first with the exact input\n", "[T_out, yout_exact, xout_exact] = control.forced_response(sys, t, y)" ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# Set up the function that plots the repsonse based on slider changes\n", "def plot_response(num_terms = 11):\n", "\n", " # get the a0 term\n", " a0 = w0/(2.0*np.pi)*np.trapz(y,t) \n", "\n", " # fill arrays with zeros - Good practice, it speeds computation in large problems\n", " a = np.zeros((num_terms,))\n", " b = np.zeros((num_terms,))\n", " integral_cos = np.zeros((len(t),num_terms))\n", " integral_sin = np.zeros((len(t),num_terms))\n", " sin_term = np.zeros((num_terms,len(t)))\n", " cos_term = np.zeros((num_terms,len(t)))\n", "\n", " # cycle through the 1 to num_terms Fourier coefficients (a_n and b_n)\n", " for n in range(num_terms):\n", "\n", " # a_n calculations\n", " integral_cos[:,n] = y * np.cos((n+1) * w0 * t) # define the integral \"interior\"\n", " a[n] = w0/np.pi * np.trapz(integral_cos[:,n], t) # solve for a_n\n", "\n", " # b_n calculations\n", " integral_sin[:,n] = y * np.sin((n+1) * w0 * t) # define the integral \"interior\"\n", " b[n] = w0/np.pi * np.trapz(integral_sin[:,n], t) # solve for b_n\n", "\n", " sin_term[n,:] = np.sin((n+1) * w0 * t) # calculate the nth sine term\n", " cos_term[n,:] = np.cos((n+1) * w0 * t) # calculate the nth cosine_term\n", "\n", "\n", " # Generate the approximate input based on the Fourier coeff. calculated above\n", " approx = np.zeros_like(t) #First fill with zeros\n", "\n", " for ii in range(len(t)):\n", " approx[ii] = a0 + np.sum(a * cos_term[:,ii], 0) + np.sum(b * sin_term[:,ii], 0)\n", " \n", " \n", " # run the simulation - now with the approximate input\n", " [T_approx,yout_approx,xout_approx] = control.forced_response(sys,t,approx)\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 = plt.figure(figsize=(9,6))\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(),family='serif',fontsize=18)\n", " plt.setp(ax.get_xmajorticklabels(),family='serif',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('Time (s)', family='serif', fontsize=22, weight='bold', labelpad=5)\n", " plt.ylabel('Position (m)', family='serif', fontsize=22, weight='bold', labelpad=10)\n", "\n", " plt.plot(t,y,':', linewidth=2, label=r'Exact Command')\n", " plt.plot(t, yout_exact,'--', linewidth=2, label=r'Exact Response')\n", "\n", " f = str(num_terms) + '-Term Approx. Resp.'\n", " plt.plot(t, yout_approx, linewidth=2, label=f)\n", "\n", " # You may need to adjust these limits for best results\n", " plt.ylim(-2.5, 3.5)\n", "\n", " leg = plt.legend(loc='upper right', fancybox=True)\n", " ltext = leg.get_texts() \n", " plt.setp(ltext,family='serif',fontsize=18)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Interactive Response Plot\n", "Now, we'll create the interactive plot. The result will have a slider that allows us to change the number of terms used in the Fourier Approximation and resulting response approximation.\n", "\n", "***Reminder***: You will have the run the notebook locally to use the interactive features. They are *not* viewable at [http://nbviewer.ipython.org](http://nbviewer.ipython.org)." ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "77c3a9d398354d4c8df285f73c0344cf", "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=(IntSlider(value=11, description='num_terms', max=22, min=1), Output()), _dom_classes=('widget-interact',))" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Call the slider interaction\n", "# num_terms changes in the number of terms in the Fourier Approx, allowing between 1 and 22\n", "interact(plot_response, num_terms=(1,22,1));" ] }, { "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." ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\n" ], "text/plain": [ "" ] }, "execution_count": 11, "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 }