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

Fourier Analysis using SymPy

\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. In this notebook, we'll use [SymPy](http://sympy.org/en/index.html), \"a Python library for symbolic mathematics.\"\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).\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", "The remainder of this notebook will focus on determining the $a_n$ and $b_n$ terms in the Fourier Expansion.\n", "
\n", "
" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# import Sympy and start \"pretty printing\"\n", "import sympy\n", "sympy.init_printing()" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# Define the sympy symbolic variables we'll need\n", "t, w0, tau_0 = sympy.symbols(['t', 'omega_0', 'tau_0'], real=True, positive=True)" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": true }, "outputs": [], "source": [ "#--------- Input your function to examine here --------\n", "\n", "# Use the sympy Piecewise function to define the square wave - This matches the one in the Figure 1 above.\n", "y = 2 + sympy.Piecewise((1, t < sympy.pi/w0), (-1, t > sympy.pi/w0))\n", "\n", "\n", "# Use the sympy Piecewise function to define the triangle wave\n", "# First define F0\n", "# F0 = sympy.symbols('F0')\n", "# y = sympy.Piecewise((F0/2*t, t < sympy.pi/w0), (-(F0/2)*t + 2*F0, t >= sympy.pi/w0))\n", "\n", "\n", "# Use the sympy Piecewise function to define a trapezoid function\n", "# y = sympy.Piecewise((3*F0*w0/(2*sympy.pi)*t, t < (2*sympy.pi/(3*w0))), (F0, t < (4*sympy.pi/(3*w0))),\n", "# (-3*F0*w0/(2*sympy.pi)*t + 3*F0, t > (4*sympy.pi/(3*w0))))" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# define the number of terms to use in the approximation\n", "num_terms = 7" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# get the a0 term\n", "a0 = w0 / (2*sympy.pi) * sympy.integrate(y, (t, 0, 2*sympy.pi/w0))" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# Define matrices of 0s to fill the the an and bn terms\n", "a = sympy.zeros(1, num_terms)\n", "b = sympy.zeros(1, num_terms)\n", "\n", "# cycle through the 1 to num_terms Fourier coefficients (a_n and b_n)\n", "for n in range(num_terms):\n", " integral_cos = y * sympy.cos((n+1)*w0*t) # define the integral \"interior\"\n", " a[n] = w0 / sympy.pi * sympy.integrate(integral_cos, (t, 0, 2*sympy.pi/w0)) # solve for a_n\n", "\n", " integral_sin = y * sympy.sin((n+1)*w0*t) # define the integral \"interior\"\n", " b[n] = w0 / sympy.pi * sympy.integrate(integral_sin, (t, 0, 2*sympy.pi/w0)) # solve for b_n" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAAkAAAAOBAMAAAAPuiubAAAALVBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAADAOrOgAAAADnRSTlMAIpm7MhCriUTv3c12VLge\nopIAAAAJcEhZcwAADsQAAA7EAZUrDhsAAABOSURBVAgdY2BUMnZgYAhjYH/BwJDKwDCTgWEWA0Oe\nA8O+ABAJBOsCgATHcxCTKwFEKoEIHgUQeYmBUYCBRYGBR4BBqrwoi4Fh37t3rxgAK5QOlzv7snYA\nAAAASUVORK5CYII=\n", "text/latex": [ "$$2$$" ], "text/plain": [ "2" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Simplify and display a0\n", "sympy.simplify(a0)" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAM8AAAAZBAMAAABkyGdeAAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAu90iEJmJZjLNVO+r\nRHY8nXFNAAAACXBIWXMAAA7EAAAOxAGVKw4bAAABVUlEQVRIDe2WrU4DQRSFT9mBliYNPxpBFofi\nDdqQgEFQRVK3DjDImor2DUCj9g2KRK5E8wLwAAQFCEIIs3PTzO0R3BSyq7qqZ7+T+2VmdpJiOz1A\n5c9jmmG3ckspSIpaRe5sONHrOjq/0JExZ6MesKxoPXd7evItDvW5MeZs1AMW0Qi4VKLmHdZuVCYM\nykZdsIh2gGkvTm4UaL/FCMKcjbpgEX0BT3mcvFGg8xkjCHM26oKDyL170XOcfNpH5ztGxpyNuuAg\nSj6A6yxO7mZo+VezhzFnoy5YRH5Fc6L+vIhwQrn7e11wvVtXnu40n+0U4I+vTR+DxgvWZVpYEbaA\ncS+KGpto6s+b8IJ1mSYifwMH0QN/xVbpwmpcXlidjbpgEa3k7l6JsI+TicqMORv1gEXkrl70YBy/\nPigPGHM26gGLSE+t6PdS9OeNXW7df7Yured/XfEDnc+x6uTMvBQAAAAASUVORK5CYII=\n", "text/latex": [ "$$\\left[\\begin{matrix}0 & 0 & 0 & 0 & 0 & 0 & 0\\end{matrix}\\right]$$" ], "text/plain": [ "[0 0 0 0 0 0 0]" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Simplify and display the an terms\n", "sympy.simplify(a)" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPsAAAAbBAMAAABIGOQXAAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAu90iVBCJdjJmmavN\nRO+PBpafAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAC10lEQVRIDY2WTWgTQRzFX7rZNEldG716UBI9\nCvGoCMlFEEG6CK3Qg81B/EIwIFWoPUQQ9WQjgsc2PRhNEAyo0IvYg3i19ipij3oT/KRF1tlkZudz\nd2cOmf//zcvvdWdnNwVGo0TnhCmTsMaXLECAisr+4oC4asqPWxF0GxDAUXvLVfLtmfR495pNvAUI\noKjNcgMHhn/6fHr8+JxNvAUIYChnMIrPn0iPv2MTbwMCGIrF59z0+IZNvA0IYCgW/0iJd7uPW8Ob\nwj8KVS1+9vIlvj6qVBBRdVeEovFuU4nfVXcrCvkstPh3OBMeW2FoILKmuyIUjc/2+9ttgYJl4KrY\nk/pU/0VPlgpbyK3LkgYCDK4IxTYfjnzv94mPJkuoKSc/M4D+lCsg8oYxuSgqin+4U2cp4bwDfJEE\nohX2vwqX+JgcwPvJ21GlgACTi6GieJnibpP4pqzp3VwH3h9dVpQkV0y88xtYaigcra01kCfGlJHk\niosnV28R37GKT3DR+CAc3/llaJs/OXSEH9wEbVsjk4DSXCIq5urDozdVF5KMJTlUWe3oac4kV1z8\nHmBFesxyC681cKaEgvS85teOttNd56enmSkunrx2njPPcO5hVerDhrxQiuui6gVk19Shud4GQYma\n4uLH6u5XidPFitQPm4M41xJV76S0Y3RJdb35/JHZ4uLdi7ckMCHpV4+Zm2tiOrwNqaWN6qrjNLPR\neOdDEPxlmnkmu3GcnOyWeXWkek8XfBtUg6FofPH+/BO2ITH8xSae3ZuNDo3R5WwUShaoCZ+haLw/\nXu0YgYK46uJ6UeiNZf6HBeo2GIrd+4zTMeIEcamO0jGhN5QTvvsP6agr5KsjFIs/4g0MOC6RO19r\nO6VPXDFVu6vk5zYNBbwHKIrF9/JbJh7XKlj2i+0HXDBVOT/XQSoK3wCKIvHl8P/8Cg6ZeFy7e+Ml\nxqpjKefzQtcChcOgqM3y4D/bqPV3AxBS6AAAAABJRU5ErkJggg==\n", "text/latex": [ "$$\\left[\\begin{matrix}\\frac{4}{\\pi} & 0 & \\frac{4}{3 \\pi} & 0 & \\frac{4}{5 \\pi} & 0 & \\frac{4}{7 \\pi}\\end{matrix}\\right]$$" ], "text/plain": [ "⎡4 4 4 4 ⎤\n", "⎢─ 0 ─── 0 ─── 0 ───⎥\n", "⎣π 3⋅π 5⋅π 7⋅π⎦" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Simplify and diplay the bn terms\n", "sympy.simplify(b)" ] }, { "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": 10, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\n" ], "text/plain": [ "" ] }, "execution_count": 10, "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 }