Fourier vs Goertzel

fourier.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def fourier(data, omega):\n",
    "    eiw = complex(np.cos(omega), np.sin(omega))\n",
    "    f = 0\n",
    "    for s in data:\n",
    "        f = eiw*f + s\n",
    "\n",
    "    return f"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def goertzel(data, omega):\n",
    "    cosine = np.cos(omega)\n",
    "    sine = np.sin(omega)\n",
    "    coeff = 2*cosine\n",
    "\n",
    "    g0 = g1 = 0\n",
    "    for s in data:\n",
    "        g0, g1 = g1, g1*coeff - g0 + s\n",
    "\n",
    "    return complex(g1 - cosine*g0, sine*g0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def next_fourier(F, sample_out, sample_in, omega, n):\n",
    "    eiw = complex(np.cos(omega), np.sin(omega))\n",
    "    eiwn = complex(np.cos(omega*(n-1)), np.sin(omega*(n-1)))\n",
    "    return eiw*(F - sample_out*eiwn) + sample_in"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "((-15.08813423269332-19.897402665645426j),\n",
       " (-15.088134232693687-19.897402665644755j),\n",
       " (3.659295089164516e-13-6.714628852932947e-13j),\n",
       " 3.161915174132446e-13)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "N = 4410\n",
    "omega = 441\n",
    "data = 2*np.random.random_sample(2*N) - 1\n",
    "\n",
    "F = fourier(data[:N], omega)\n",
    "G = goertzel(data[:N], omega)\n",
    "\n",
    "F, G, F-G, abs(F)-abs(G)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "((12.24881925474038-21.688659683576557j),\n",
       " (12.248819254740404-21.688659683576585j),\n",
       " (-2.4868995751603507e-14+2.8421709430404007e-14j),\n",
       " -3.907985046680551e-14)"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "F2 = fourier(data[1:N+1], omega)\n",
    "FF2 = next_fourier(F, data[0], data[N], omega, N)\n",
    "F2, FF2, F2-FF2, abs(F2)-abs(FF2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1000 loops, best of 3: 539 µs per loop\n",
      "1000 loops, best of 3: 901 µs per loop\n",
      "100000 loops, best of 3: 6.49 µs per loop\n"
     ]
    }
   ],
   "source": [
    "%timeit fourier(data[1:N+1], omega)\n",
    "%timeit goertzel(data[1:N+1], omega)\n",
    "%timeit next_fourier(F, data[0], data[N], omega, N)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "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.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}