FourierTransform.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Discrete Fourier transform"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "import FFTW\n",
    "using BenchmarkTools\n",
    "using PyCall"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "const x = rand(2 ^ 13);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "naivedft (generic function with 1 method)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function dftmatrix(x::AbstractArray)\n",
    "    n = length(x)\n",
    "    W = ones(Complex{eltype(x)}, n, n)\n",
    "    ω = exp(-1im * 2π / n)\n",
    "    for i = 1:n, j = 1:n\n",
    "        W[j, i] = ω ^ ((i - 1) * (j - 1))\n",
    "    end\n",
    "    W\n",
    "end\n",
    "\n",
    "function naivedft(x::AbstractArray)\n",
    "    dftmatrix(x) * x\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "true"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "isapprox(FFTW.fft(x), naivedft(x))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "fft_dit (generic function with 1 method)"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function fft_dit(x::AbstractArray)\n",
    "    n = length(x)\n",
    "    \n",
    "    if n == 1\n",
    "        return x[1]\n",
    "    end\n",
    "    \n",
    "    ω = exp(-2π * im / n)\n",
    "    \n",
    "    x_even = fft_dit(x[1:2:end])\n",
    "    x_odd = fft_dit(x[2:2:end])\n",
    "\n",
    "    factor = ω .^ collect(0:n-1)\n",
    "    halfway = n ÷ 2\n",
    "    vcat(\n",
    "        x_even .+ factor[1:halfway] .* x_odd,\n",
    "        x_even .+ factor[halfway+1:end] .* x_odd\n",
    "    )\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "true"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "isapprox(FFTW.fft(x), fft_dit(x))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "fft (generic function with 1 method)"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function bit_reversal!(x::AbstractArray)\n",
    "    n = length(x)\n",
    "    j = 1\n",
    "    for i = 1:n-1\n",
    "        if i < j\n",
    "            x[i], x[j] = x[j], x[i]\n",
    "        end\n",
    "        k = n ÷ 2\n",
    "        while k < j\n",
    "            j -= k\n",
    "            k ÷= 2\n",
    "        end\n",
    "        j += k\n",
    "    end\n",
    "end\n",
    "\n",
    "function fft!(x::AbstractArray{<:Complex})\n",
    "    n = length(x)\n",
    "    bit_reversal!(x)\n",
    "    blocklen = 2\n",
    "    while blocklen <= n\n",
    "        middlepoint = blocklen ÷ 2\n",
    "        for innerindex = 1:middlepoint\n",
    "            m = -2 * im * (innerindex - 1) / blocklen\n",
    "            ω = exp(m * π)\n",
    "            for block = 1:n÷blocklen\n",
    "                i1 = (block - 1) * blocklen + innerindex\n",
    "                i2 = i1 + middlepoint\n",
    "                a = x[i1]\n",
    "                b = x[i2] * ω\n",
    "                x[i1] = a + b\n",
    "                x[i2] = a - b\n",
    "            end\n",
    "        end\n",
    "        blocklen *= 2\n",
    "    end\n",
    "end\n",
    "\n",
    "function fft(x::AbstractArray{<:Real})\n",
    "    y = complex(copy(x))\n",
    "    fft!(y)\n",
    "    y\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "true"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "isapprox(FFTW.fft(x), fft(x))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Julia Version 1.4.2\n",
      "Commit 44fa15b150* (2020-05-23 18:35 UTC)\n",
      "Platform Info:\n",
      "  OS: macOS (x86_64-apple-darwin17.7.0)\n",
      "  CPU: Intel(R) Core(TM) i7-7700K CPU @ 4.20GHz\n",
      "  WORD_SIZE: 64\n",
      "  LIBM: libopenlibm\n",
      "  LLVM: libLLVM-8.0.1 (ORCJIT, skylake)\n"
     ]
    }
   ],
   "source": [
    "versioninfo()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "BenchmarkTools.Trial: \n",
       "  memory estimate:  1.00 GiB\n",
       "  allocs estimate:  6\n",
       "  --------------\n",
       "  minimum time:     4.578 s (0.01% GC)\n",
       "  median time:      4.601 s (0.39% GC)\n",
       "  mean time:        4.601 s (0.39% GC)\n",
       "  maximum time:     4.624 s (0.77% GC)\n",
       "  --------------\n",
       "  samples:          2\n",
       "  evals/sample:     1"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "@benchmark naivedft(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "BenchmarkTools.Trial: \n",
       "  memory estimate:  256.16 KiB\n",
       "  allocs estimate:  4\n",
       "  --------------\n",
       "  minimum time:     49.349 μs (0.00% GC)\n",
       "  median time:      134.119 μs (0.00% GC)\n",
       "  mean time:        145.944 μs (10.17% GC)\n",
       "  maximum time:     3.462 ms (94.28% GC)\n",
       "  --------------\n",
       "  samples:          10000\n",
       "  evals/sample:     1"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "FFTW.set_num_threads(1)\n",
    "plan = FFTW.plan_fft(x)\n",
    "@benchmark (plan * x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "BenchmarkTools.Trial: \n",
       "  memory estimate:  15.54 MiB\n",
       "  allocs estimate:  147470\n",
       "  --------------\n",
       "  minimum time:     12.299 ms (0.00% GC)\n",
       "  median time:      12.520 ms (0.00% GC)\n",
       "  mean time:        13.259 ms (5.70% GC)\n",
       "  maximum time:     16.007 ms (14.61% GC)\n",
       "  --------------\n",
       "  samples:          377\n",
       "  evals/sample:     1"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "@benchmark fft_dit(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "BenchmarkTools.Trial: \n",
       "  memory estimate:  192.16 KiB\n",
       "  allocs estimate:  4\n",
       "  --------------\n",
       "  minimum time:     266.433 μs (0.00% GC)\n",
       "  median time:      309.516 μs (0.00% GC)\n",
       "  mean time:        317.771 μs (2.59% GC)\n",
       "  maximum time:     2.481 ms (87.26% GC)\n",
       "  --------------\n",
       "  samples:          10000\n",
       "  evals/sample:     1"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "@benchmark fft(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "┌ Warning: both Julia and NumPy are linked with MKL, which may cause conflicts and crashes (#433).\n",
      "└ @ PyCall /Users/hjkim/.julia/packages/PyCall/zqDXB/src/numpy.jl:73\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "BenchmarkTools.Trial: \n",
       "  memory estimate:  130.00 KiB\n",
       "  allocs estimate:  42\n",
       "  --------------\n",
       "  minimum time:     108.154 μs (0.00% GC)\n",
       "  median time:      195.517 μs (0.00% GC)\n",
       "  mean time:        214.026 μs (4.15% GC)\n",
       "  maximum time:     6.828 ms (44.58% GC)\n",
       "  --------------\n",
       "  samples:          10000\n",
       "  evals/sample:     1"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np = pyimport(\"numpy\")\n",
    "@benchmark np.fft.fft(x)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.4.2",
   "language": "julia",
   "name": "julia-1.4"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.4.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}