{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c271abd4-6464-4d0d-a598-18623c080570",
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib widget\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from setupFigure import SetupFigure"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "33012d14-c18a-4e61-89c5-10ac079f4121",
   "metadata": {},
   "outputs": [],
   "source": [
    "def dft_coeff(x):\n",
    "    \"\"\"\n",
    "    Evaluate c_n = 1/N*sum_{k=0}^{N-1}x_k exp(-i*2*pi*nk/N)\n",
    "    :param x: samples of function with length N\n",
    "    \"\"\"\n",
    "    ns = np.size(x)\n",
    "    c = np.zeros(ns, dtype=np.complex64)  # zero coefficients before summing\n",
    "    for n in range(ns):\n",
    "        for k in range(ns):\n",
    "            arg = -2*np.pi*1j*n*k/ns\n",
    "            c[n] = c[n]+x[k]*np.exp(arg)\n",
    "        c[n] = c[n]/ns\n",
    "    return c"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "80d7b404-00c7-49e8-9061-a262e7b84bb7",
   "metadata": {},
   "outputs": [],
   "source": [
    "def dft_inverse(c):\n",
    "    \"\"\"\n",
    "    Evaluate x_k = sum_{n=0}^{N-1}c_n exp(i*2*pi*nk/N)\n",
    "    :param c: complex DFT coefficients\n",
    "    \"\"\"\n",
    "    ns = np.size(c)\n",
    "    x = np.zeros(ns)   # zero values before summing\n",
    "    for k in range(ns):\n",
    "        for n in range(ns):\n",
    "            arg = +2*np.pi*1j*n*k/ns\n",
    "            x[k] = x[k]+np.real(c[n]*np.exp(arg))   # since x[k] is real, we only sum the real parts\n",
    "    return x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b24ed02f-2388-444a-943c-3d6341f9490a",
   "metadata": {},
   "outputs": [],
   "source": [
    "def dft_inverse_any_time(c, period, dt, nt, mc=1):\n",
    "    \"\"\"\n",
    "    Evaluate x(k*dt) = sum_{n=0}^{N-1}c_n exp(i*2*pi*n*(k*dt)/T)\n",
    "    :param c: complex DFT coefficients\n",
    "    :param period: period of time series\n",
    "    :param dt: new time sampling interval\n",
    "    :param nt: number of time samples\n",
    "    :param mc: take mc*N Fourier coefficients\n",
    "    \"\"\"\n",
    "    ns = np.size(c)\n",
    "    x = np.zeros(nt)\n",
    "    for k in range(nt):\n",
    "        for n in range(mc*ns):\n",
    "            arg = +2*np.pi*1j*n*(k*dt)/period\n",
    "            dum, n2 = divmod(n,ns)                   # use periodicity c[n+N] = c[n]\n",
    "            x[k] = x[k]+np.real(c[n2]*np.exp(arg))   # since x[k] result is real, we only sum the real parts\n",
    "    return x"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4112da9b-5ff0-4306-a4a4-d7cf3924b274",
   "metadata": {},
   "source": [
    "### The sampled Gaussian function"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "941b218a-4bbf-4588-962d-177a300a1051",
   "metadata": {},
   "source": [
    "Setup of Gaussian function $x(t)=\\exp(-(t-3)^2)$ in range $0 < t < 6$. We take $\\Delta t=0.1$ producing 60 intervals and N=61 samples.\n",
    "Note that the period of the function is then $T=N\\Delta t=6.1$ and not $6.0$; hence $x(T=N\\Delta t) = x(0)$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bd7afa4c-1091-4b0d-a211-38933215d782",
   "metadata": {},
   "outputs": [],
   "source": [
    "tmax = 6.0\n",
    "dt = 0.1\n",
    "ns = int(tmax/dt)+1\n",
    "period = ns*dt\n",
    "x = np.zeros(ns)\n",
    "t = np.linspace(0, tmax, ns)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "03f1b9cb-be96-4e54-a479-1a3ed2da3d0a",
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.exp(-(t-3)**2)   # compute ns samples of Gaussian"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fe9d7345-81b0-4c49-afab-8d15acf62bbe",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig1, ax1 = SetupFigure(10, 5, \"Time [s]\", \"x(t)\", \"Gaussian, DT=0.1, TMAX=6\", 14)\n",
    "ax1.plot(t, x, color='blue', ls=':', marker='o', markersize=4.0)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ab863cff-159d-4586-925e-580d5aedcaaf",
   "metadata": {},
   "source": [
    "### The Fourier coefficients of the sampled Gaussian"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1e2933b2-adbf-4c20-96a0-34d5c0bb0077",
   "metadata": {},
   "source": [
    "Compute the Fourier coefficients using $c_n = \\frac{1}{N}\\sum\\limits_{k=0}^{N-1} x_k \\exp(-2\\pi i nk/N)$. The frequencies to $c_n$ are $f_n = n/T$ or $\\omega_n = 2\\pi n/T$. The frequency to $c_{N-1}$ is $\\frac{N-1}{T}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "015ba673-f453-46c8-b034-03bf72ffe965",
   "metadata": {},
   "outputs": [],
   "source": [
    "c = dft_coeff(x)\n",
    "f = np.linspace(0, (ns-1)/period, ns)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1ecdfb4c-28c0-447c-8110-52b3c6148992",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig2, ax2 = SetupFigure(10, 5, \"Frequency [Hz]\", \"c_n\", \"Fourier coefficients of Gaussian, DT=0.1, TMAX=6\", 14)\n",
    "ax2.plot(f, np.absolute(c), color='black', ls='--', marker='d', label='abs')\n",
    "ax2.plot(f, np.real(c), color='red', ls=':', marker='*', markersize=4.0, label='real')\n",
    "ax2.plot(f, np.imag(c), color='blue', ls=':', marker='o', markersize=4.0, label='imag')\n",
    "ax2.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d8e7ff1d-c43b-418b-9d04-e104fd127a3f",
   "metadata": {},
   "source": [
    "### The recovered Gaussian at the original samples"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d9fbd31d-af26-4e7c-8db9-eba38cd2d9b0",
   "metadata": {},
   "source": [
    "Recover the samples by inverse discrete Fourier transform: $x_k = \\sum\\limits_{n=0}^{N-1} c_n \\exp(2\\pi i nk/N)$ belonging to times $t_k=k\\Delta t$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ad695630-ed05-4508-abb7-802b1ef0556f",
   "metadata": {},
   "outputs": [],
   "source": [
    "xr = dft_inverse(c)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cf73ed2b-5571-4d8e-bfb1-aa917dc351bc",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig3, ax3 = SetupFigure(10, 5, \"Time [s]\", \"xr(t)\", \"Recovered Gaussian, DT=0.1, TMAX=6\", 14)\n",
    "ax3.plot(t, xr , color='blue', ls='--',  marker='o', markersize=4.0)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d741cc78-d48e-4ffd-9029-811367b4b972",
   "metadata": {},
   "source": [
    "### The recovered Gaussian at original sampling interval up to TMAX=12"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5008030a-6f9f-4613-9937-ebbb6e0d9298",
   "metadata": {},
   "source": [
    "Evaluate the inverse transform at times beyond $T$ using $x_k = \\sum\\limits_{n=0}^{N-1} c_n \\exp(2\\pi in t/T)$. Here, we take the same sampling interval but extend the time range to $12$. We get a periodic continuation of the Gaussian."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fcef58bc-051d-49e9-b37f-40b6bc91a87e",
   "metadata": {},
   "outputs": [],
   "source": [
    "tmax = 12\n",
    "nt = int(tmax/dt)+1\n",
    "tlong = np.linspace(0, tmax, nt)\n",
    "xrlong = dft_inverse_any_time(c, period, dt, nt)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8e527baa-c2ea-43a0-b10c-8366048d6eac",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig4, ax4 = SetupFigure(10, 5, \"Time [s]\", \"xrlong(t)\", \"Recovered Gaussian, DT=0.1, TMAX=12\", 14)\n",
    "ax4.plot(tlong, xrlong , color='blue', ls='--',  marker='o', markersize=4.0)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aa22d27c-a959-4e3f-b7bf-5fc1561d9956",
   "metadata": {},
   "source": [
    "### The recovered Gausian with dense sampling (DT=0.01)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "08123a90-df5b-4bf8-aade-4819b55d73db",
   "metadata": {},
   "source": [
    "Here, we also evaluate the Fourier series for times in between the original samples by choosing times $t=0.1\\,k\\Delta t$. Evidently, this is not a suitable interpolation method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "10458503-2f03-4a16-b20e-8fb4aa0c99b0",
   "metadata": {},
   "outputs": [],
   "source": [
    "tmax = 6\n",
    "dtnew = 0.1*dt\n",
    "nt = int(tmax/dtnew)+1\n",
    "tdense = np.linspace(0, tmax, nt)\n",
    "xdense = dft_inverse_any_time(c, period, dtnew, nt, mc=1)        # dense sampling"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2cc0d1c9-f2a6-4bfd-84d9-bd8af31818d0",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig5, ax5 = SetupFigure(10, 5, \"Time [s]\", \"xrdense(t)\", \"Recovered Gaussian, DT=0.01, TMAX=6\", 14)\n",
    "ax5.plot(tdense, xdense , color='blue', ls=':',  marker='o', markersize=4.0)\n",
    "ax5.plot(t, x , color='red', ls='--',  marker='d', markersize=4.0)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "061d0254-a311-4a70-b26f-f51f1ae523ac",
   "metadata": {},
   "source": [
    "### The recovered Gausian with dense sampling and high order Fourier coefficients (DT=0.01, mc*N-1)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "20d7b156-6a75-43e7-b5d1-04d9663b15be",
   "metadata": {},
   "source": [
    " Here, we again use dense sampling and, in addition, sum up the coefficients to $mc*N-1$ thus introducing higher frequencies. By doing this, we aproximate the time function $x(t)=\\sum_{k=0}^{N-1} x_k\\delta(t-k\\Delta t)$ with high peaks at the samples and zero values in between. Play with increasingly higher values of $mc$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "902ecafc-1f7e-4b16-a518-19f1969e2990",
   "metadata": {},
   "outputs": [],
   "source": [
    "xhigh = dft_inverse_any_time(c, period, dtnew, nt, mc=60)         # dense sampling and higher freqeuncies"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0d3e681e-49e5-4549-b864-bbf6df6f61de",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig6, ax6 = SetupFigure(10, 5, \"Time [s]\", \"xrhigh(t)\", \"Recovered Gaussian, DT=0.01, TMAX=6, high frequencies\", 14)\n",
    "ax6.plot(tdense, xhigh , color='blue', ls=':',  marker='o', markersize=4.0)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.12.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}