Fourier series, magnitude, phase and shift.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Fourier series, magnitude, phase and shift\n",
    "## Introduction\n",
    "In this notebook, the basics of Fourier analysis are discussed. With Fourier analysis, we can:\n",
    "* find out the relative importance of different frequencies by comparing their magnitudes\n",
    "* find out the phase of a frequency, thus when in the period a frequency has its maximum value\n",
    "* find out the phase shift between two time series\n",
    "\n",
    "In this tutorial, the basics are shown. We start by loading the necessary Python modules."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Populating the interactive namespace from numpy and matplotlib\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "from scipy import signal\n",
    "import matplotlib.pyplot as pl\n",
    "%pylab inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, we generate a time variable `t` that has 10 periods (`20.*np.pi`) distributed over 10,000 (`n`) points."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "n = 100000\n",
    "t = np.linspace(0., 20*np.pi, n)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then, we generate two time series. Series `y1` is a superposition of two sinus waves, series `y2` is `y1` with a linear trend added, and `y2d` is `y2`, but with a linear trend removed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/usr/local/lib/python2.7/site-packages/scipy/linalg/basic.py:884: RuntimeWarning: internal gelsd driver lwork query error, required iwork dimension not returned. This is likely the result of LAPACK bug 0038, fixed in LAPACK 3.2.2 (released July 21, 2010). Falling back to 'gelss' driver.\n",
      "  warnings.warn(mesg, RuntimeWarning)\n"
     ]
    }
   ],
   "source": [
    "y1 = np.sin(t) + 0.5*np.sin(4.*t)\n",
    "y2 = np.sin(t) + 0.5*np.sin(4.*t) + 3./(20.*np.pi)*t\n",
    "y2d = signal.detrend(y2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the code below we define a function to calculate the Fast Fourier Transform of a series. The Fast Fourier Transform is defined as:\n",
    "___\n",
    "$$ Y_k = \\sum_{n=0}^{N-1} y_n e^{-i2 \\pi k \\frac{n}{N}} \\qquad k=0,\\dots,N-1$$\n",
    "___"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def plot_spectrum(y, nfig):\n",
    "    # An FFT of real data.\n",
    "    yfft = np.fft.rfft(y)\n",
    "\n",
    "    # Division by the number of samples, as the FFT\n",
    "    # multiplies the values by that number.\n",
    "    yfft /= n\n",
    "\n",
    "    # Correct the data for the absense of the complex conjugate.\n",
    "    yfft[1:-1] *= 2.\n",
    "   \n",
    "    # Plot the series and its spectrum.\n",
    "    pl.figure(nfig)\n",
    "    pl.subplot(211)\n",
    "    pl.plot(t, y)\n",
    "    pl.subplot(212)\n",
    "    pl.semilogx(abs(yfft))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "plot_spectrum(y1 , 1)\n",
    "plot_spectrum(y2 , 1)\n",
    "plot_spectrum(y2d, 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Calculate the magnitude and phase of these two series.\n",
    "y3 = np.sin(t)\n",
    "y4 = 0.5*np.cos(t)\n",
    "\n",
    "plot_spectrum(y3, 2)\n",
    "plot_spectrum(y4, 2)\n",
    "\n",
    "def print_magnitude_phase(t, y):\n",
    "    yfft = np.fft.rfft(y) / n\n",
    "    yfft[1:-1] *= 2.\n",
    "\n",
    "    yfft10 = yfft[10]\n",
    "\n",
    "    # The magnitude is (Re(yfft)**2 + Im(yfft)**2)**(1/2)\n",
    "    yfft10_magnitude = abs(yfft10)\n",
    "\n",
    "    # The mathematical phase is the shift compared to the cosine.\n",
    "    yfft10_phase = np.arctan2(yfft10.imag, yfft10.real)\n",
    "    print 'Value = {0}, Magnitude = {1}, Phase = {2}'.format(yfft10,\\\n",
    "                                                             yfft10_magnitude,\\\n",
    "                                                             yfft10_phase)\n",
    "\n",
    "    return yfft10, yfft10_magnitude, yfft10_phase\n",
    "\n",
    "y3_10, y3_10_magnitude, y3_10_phase = print_magnitude_phase(t, y3)\n",
    "y4_10, y4_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Calculate the magnitude and phase of these two series.\n",
    "y3 = np.sin(t)\n",
    "y4 = 0.5*np.cos(t)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def print_magnitude_phase(t, y):\n",
    "    yfft = np.fft.rfft(y) / n\n",
    "    yfft[1:-1] *= 2.\n",
    "\n",
    "    yfft10 = yfft[10]\n",
    "\n",
    "    # The magnitude is (Re(yfft)**2 + Im(yfft)**2)**(1/2)\n",
    "    yfft10_magnitude = abs(yfft10)\n",
    "\n",
    "    # The mathematical phase is the shift compared to the cosine.\n",
    "    yfft10_phase = np.arctan2(yfft10.imag, yfft10.real)\n",
    "    print 'Value = {0}, Magnitude = {1}, Phase = {2}'.format(yfft10,\\\n",
    "                                                             yfft10_magnitude,\\\n",
    "                                                             yfft10_phase)\n",
    "\n",
    "    return yfft10, yfft10_magnitude, yfft10_phase"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "y3_10, y3_10_magnitude, y3_10_phase = print_magnitude_phase(t, y3)\n",
    "y4_10, y4_10_magnitude, y4_10_phase = print_magnitude_phase(t, y4)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}