 Fourier Analysis I Problems (Problem Version).ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "0f2abc99",
   "metadata": {},
   "source": [
    "## Problem Set for Lecture I: Fourier Analysis (Bryan's Version) \n",
    "\n",
    "Bryan Scott, CIERA/Northwestern\n",
    "\n",
    "DSFP Session 20: Time Series Methods in Astronomy\n",
    "\n",
    "Adapted and Expanded from Exercises Originally given by L. Bugnet (Flatiron Institute) at Session 13 of the Data Science Fellowship Program \n",
    "\n",
    "version 0.1"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "86ab2b16",
   "metadata": {},
   "source": [
    "In this notebook, you will find four problems intended to build intuition about the development of the Fourier Transform from Fourier series. Because Fourier series are expressed in terms of coefficients computed from integral expressions, I have given you some expressions and asked you to compute others. \n",
    "\n",
    "This can be done either by hand or with a tool like sympy or Mathematica. Do not spend too much time on the integrals, but if you have not done a Fourier Series by hand before (or haven't in a long time) make sure you are able to work at least one of these out for yourself to see how it works. \n",
    "\n",
    "Many of the interesting or important features of Fourier analysis are asked in the form of questions in this notebook. Make sure you take a few minutes to think about them and discuss with your neighbor. The last problem will model applications of Fourier analysis on real data from the $\\textit{Kepler}$ mission. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0855f3cf",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8e91b27e",
   "metadata": {},
   "source": [
    "# Problem 1: Fourier Series for discontinuous functions "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0d7cd8bb",
   "metadata": {},
   "source": [
    "## Problem 1: Use the formula defining the Fourier Series to write down a series approximation to the square wave defined by:\n",
    "\n",
    "$$\n",
    "f(t) \\begin{cases} \n",
    "      +1 & 0 \\leq x\\leq \\pi \\\\\n",
    "      0 & \\pi \\leq x\\leq 2\\pi \\\\\n",
    "   \\end{cases}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "448b4875",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cf52d482",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot f over 3T-wide intervals:\n",
    "# you can use this template or np.piecewise/lambda function here for a slightly fancier solution\n",
    "\n",
    "x1=\n",
    "x2=\n",
    "\n",
    "plt.figure()\n",
    "for ii in range(0,4):\n",
    "    plt.plot()\n",
    "    plt.plot()\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "020e007f",
   "metadata": {},
   "source": [
    "### Part A) Show that the Fourier Series expansion for this function is given by:\n",
    "\n",
    "$$\n",
    "F(t) = 1/2 + \\Sigma_{n=1}^{\\infty} \\left( \\frac{sin n\\pi}{n \\pi} cos(n t) + \\frac{1 - cos(n\\pi)}{n \\pi} sin(n t) \\right)\n",
    "$$\n",
    "\n",
    "by computing the integrals for $a_n$, $b_n$, and $a_0$."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b201f6e2",
   "metadata": {},
   "source": [
    "### Part B) Now write a function to plot the Fourier Series expressions for different numbers of sine and cosine functions.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6c96f0de",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Overplot Fourier Serie F:\n",
    "x = \n",
    "a0= 1./2.\n",
    "Na =  #number of trig functions to consider in the Fourier Serie\n",
    "\n",
    "#Fucntion that calculates F from a0, an, bn:\n",
    "def f_N(x,N, a0):\n",
    "    \n",
    "    for n in na:\n",
    "        an = \n",
    "        bn = \n",
    "        f_N \n",
    "    return f_N\n",
    "    \n",
    "    \n",
    "for N in Na:\n",
    "    plt.plot()\n",
    "plt.xlabel('x')\n",
    "plt.legend(loc='best')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "db937e3c",
   "metadata": {},
   "source": [
    "## Part 2: Repeat the above for the triangle wave that I showed in my lecture and replicate my plots\n",
    "\n",
    "You will need the fact that the Fourier series for the triangle wave is\n",
    "\n",
    "$$\n",
    "\\tilde{f}(x) = \\frac{2}{\\pi(2n -1))^2} cos( 2 \\pi (2n -1)x)\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "421c1a62",
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot the triangle wave function here\n",
    "\n",
    "x = \n",
    "fN = \n",
    "\n",
    "plt.plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e8927d98",
   "metadata": {},
   "outputs": [],
   "source": [
    "  \n",
    "#Overplot Fourier Serie F:\n",
    "x = np.linspace(-1, 1, 1000)\n",
    "a0=1./4\n",
    "Na = [ 1, 3, 5, 11, 101] #number of sine function to consider in the Fourier Serie\n",
    "\n",
    "#Fucntion that calculates F from a0, an, bn:\n",
    "def f_N(x,N, a0):\n",
    "\n",
    "    return f_N\n",
    "    \n",
    "    \n",
    "for N in Na:\n",
    "    \n",
    "\n",
    "plt.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "387b6a34",
   "metadata": {},
   "source": [
    "## Part 3: What do you notice about the behavior of the approximations you just derived? Why do you think this happens?\n",
    "\n",
    "Hint: Take the derivatives of the two functions. What differs between them? We'll come back to this in the afternoon session. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c7484fe8",
   "metadata": {},
   "source": [
    "# Problem 2: Fourier Transforms of Periodic and Non-periodic Functions\n",
    "\n",
    "Plot the periodic function\n",
    "\n",
    "$$\n",
    "f(x) = \\sin(6\\pi x)\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3072e097",
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.linspace(0,3, )\n",
    "y = \n",
    "plt.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6e5076d9",
   "metadata": {},
   "source": [
    "## Part 1: Visualizing the Fourier Transforms of a Periodic and Non-periodic Function \n",
    "\n",
    "Now construct a counterpart function that decays exponentially. Introducing an exponential decay is a common way of making a non-periodic function periodic. We will study this during the afternoon lecture, but can you see why this works from the definition of the Fourier transform? Why this works physically? And why might it mean that damped physical systems have EOM that can be solved by Fourier transform methods?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ae98ea82",
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.linspace(0,3, 10000)\n",
    "fig, ax = plt.subplots(1,2, figsize=(9,4))\n",
    "\n",
    "#plot function f1:\n",
    "ax[0].plot()\n",
    "\n",
    "ax[0].set_xlabel('x')\n",
    "ax[0].set_ylabel(r'$f_1(x)$')\n",
    "\n",
    "#plot function f2:\n",
    "ax[1].plot()\n",
    "\n",
    "ax[1].set_xlabel('x')\n",
    "ax[1].set_ylabel(r'$f_1(x)$')\n",
    "\n",
    "plt.tight_layout();"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7c154ac9",
   "metadata": {},
   "source": [
    "Now plot the fourier transforms of each function. The analytic fourier transform for my choice of exponential decay is:\n",
    "\n",
    "$$\n",
    "\\frac{6\\pi}{(2\\pi+2\\pi i \\nu)^2+(6 \\pi)^2}\n",
    "$$.\n",
    "\n",
    "How do you deal with the imaginary term in this equation? What will the Fourier transform of the periodic function be? After plotting the fourier transforms of both functions, discuss what you notice. What does this tell you about the Fourier transforms of periodic and non-periodic functions, respectively? "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b2dc1b36",
   "metadata": {},
   "outputs": [],
   "source": [
    "nu = np.linspace(0,50, 10000)\n",
    "\n",
    "#we approximate FT1 with:\n",
    "FT1 = np.zeros(10000)\n",
    "FT1[]=\n",
    "\n",
    "#plot FT1:\n",
    "plt.plot(nu, )\n",
    "\n",
    "\n",
    "#plot FT2\n",
    "FT2 = \n",
    "plt.plot(, color='crimson')\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7cf2eab1",
   "metadata": {},
   "source": [
    "## Part 2: Recreate My Fourier Transform as the Limiting Function Figure \n",
    "\n",
    "In the lecture portion, I showed a figure where I increased the number of c_n terms plotted and derived a \"bounding\" function for the Fourier Series of the RECT function. Recreate this figure by writing down an expression for the Fourier coefficients and then plotting them for different values of n. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "689a9fe0",
   "metadata": {},
   "outputs": [],
   "source": [
    "def c_n(n, T):\n",
    "    if n == 0:\n",
    "        return \n",
    "    else:\n",
    "        return \n",
    "    \n",
    "c_n_list = []\n",
    "nlist = \n",
    "for i in nlist:\n",
    "    c_n_list.append()\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a1cd5661",
   "metadata": {},
   "outputs": [],
   "source": [
    "for i,c_n in enumerate(c_n_list):\n",
    "    plt.vlines()\n",
    "    \n",
    "plt.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "07f9d5b4",
   "metadata": {},
   "source": [
    "### Now repeat this exercise by taking the fft of a piecewise defined rect function. \n",
    "\n",
    "Hint: You may find it helpful to make use of the np.fft.fft and np.fft.fftfreq helper function, as well as using np.piecewise with lambda function syntax. See the example in the np.piecewise documentation: https://numpy.org/doc/stable/reference/generated/numpy.piecewise.html "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3c7c6c4a",
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.linspace(-1.5, 1.5, )\n",
    "\n",
    "t = np.piecewise()\n",
    "\n",
    "freq = \n",
    "sp = \n",
    "plt.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1ddb21b8",
   "metadata": {},
   "source": [
    "# Problem 3:  Power Spectral Densities "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "830b2d3f",
   "metadata": {},
   "source": [
    "The power spectral density is defined as:\n",
    "\n",
    "$$\n",
    "PSD(\\nu) = \\lim_{T \\rightarrow \\infty} \\frac{1}{T} \\|{\\tilde{f}(\\nu)}\\|^2\n",
    "$$\n",
    "\n",
    "In this problem, you will investigate the origin of the 1/T term and the utility of computing this quantity as opposed to simply computing the Fourier transform itself. Luckily for us, the above term can be computed using the scipy.welch function. \n",
    "\n",
    "Begin by plotting a 1D periodic signal (I've provided some reproducible parameter choices for you), then compute the FFT and PSD as you increase the number of points in time that you 'sample the signal at'. We will discuss what sampling means more formally in the afternoon.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4a39f9c9",
   "metadata": {},
   "outputs": [],
   "source": [
    "fs = 10e3\n",
    "N = 1e3\n",
    "amp = 2*np.sqrt(2)\n",
    "freq = 1234.0\n",
    "time = np.arange(N) / fs\n",
    "x = \n",
    "f, Pxx_den = signal.welch()\n",
    "\n",
    "fig, ax = plt.subplots(2,3, figsize=(9,6))\n",
    "\n",
    "ax[0,0].plot()\n",
    "\n",
    "# fill in the appropriate plots \n",
    "\n",
    "\n",
    "ax[0,1].set_xlabel('frequency [Hz]')\n",
    "ax[0,1].set_ylabel('PSD')\n",
    "ax[0,1].set_xlabel('frequency [Hz]')\n",
    "ax[0,2].set_ylabel('FFT')\n",
    "\n",
    "\n",
    "# now repeat the exercise with a higher sampling rate N, \n",
    "# what do you notice about the behavior of the PSD and FFT respectively?\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7688a355",
   "metadata": {},
   "source": [
    "# Problem 4: Application to Real Data "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f450213a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# pip install lightkurve \n",
    "from lightkurve import search_targetpixelfile\n",
    "import lightkurve as lk"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "53d709e9",
   "metadata": {},
   "source": [
    "### Part a) Load a Kepler Lightcurve \n",
    "\n",
    "The lightkurve package contains a number of convenience methods for loading and working with light curve data. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6e4860d5",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Download KIC 3541346 lightcurve (flux VS time):\n",
    "search_result = lk.search_lightcurve('KIC 3541346', author='Kepler', cadence='long')\n",
    "print(search_result)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "83853590",
   "metadata": {},
   "outputs": [],
   "source": [
    "lc_collection = search_result.download_all()\n",
    "#we stitch all the quarters together:\n",
    "lc_stitched = lc_collection.stitch()\n",
    "lc_stitched.plot();"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "832317f9",
   "metadata": {},
   "source": [
    "### Part b) How to handle gaps in time series\n",
    "\n",
    "The results we saw in the lecture all depend on even sampling of the time series and we will see more formally well posed methods for dealing with gaps in data in tomorrow's lectures. For now, we need to (evenly) fill in the gaps in the light curve that we downloaded so that we can apply our Fourier analysis tools. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9ea67b8e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Let's start by seeing where the largest gaps are\n",
    "time = lc_stitched.time.to_value(format='jd')\n",
    "flux = lc_stitched.flux.to_value()\n",
    "\n",
    "plt.figure()\n",
    "plt.plot(np.diff(time), color='dodgerblue')\n",
    "\n",
    "print(np.median(np.diff(time)))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "29a8beaa",
   "metadata": {},
   "source": [
    "It looks like time differences > 0.02 should be considered gaps. We also need to resample the time series onto an *evenly* spaced grid. Two options: \n",
    "\n",
    "1) just resample to the nearest point and append zero \n",
    "2) generate an interpolating function and resample onto an even time grid\n",
    "\n",
    "Give both a try and compare the resulting time series"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "750ae147",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "plt.plot(new_time, new_flux, color='dodgerblue')\n",
    "plt.plot(time, flux, color = 'orange')\n",
    "plt.xlabel('Time [days]')\n",
    "plt.ylabel('Flux')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "28ac022b",
   "metadata": {},
   "source": [
    "### Part c) Compute the FFT and estimate the period. As we saw before, the FFT depends on the number of sample points and the bandwidth. We'll understand the theory of this better in the afternoon, but compute the FFT for different time windows and compare the results. What happens?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a8b41338",
   "metadata": {},
   "outputs": [],
   "source": [
    "# ---Fill and excecute this cell---\n",
    "time_res=(new_time[1]-new_time[0])#/len(time)\n",
    "# Number of sample points\n",
    "N = len(new_flux)\n",
    "\n",
    "print(N)\n",
    "# sample spacing\n",
    "T = time_res\n",
    "\n",
    "yf = np.fft.fft(new_flux)\n",
    "xf = np.fft.fftfreq(N, T)[1:N//2]\n",
    "\n",
    "print(xf)\n",
    "print(yf)\n",
    "\n",
    "fig, ax = plt.subplots(1,2, figsize=(9,4))\n",
    "ax[0].set_yscale('log')\n",
    "ax[0].plot(xf, np.abs(yf[1:N//2]), color='dodgerblue')\n",
    "\n",
    "wh = np.argwhere((np.abs(yf[10:N//2])==np.max(np.abs(yf[10:N//2]))) )\n",
    "frequency_max = xf[9:][wh]\n",
    "\n",
    "print(frequency_max)\n",
    "\n",
    "Period = 1/frequency_max[0][0]\n",
    "ax[0].plot(frequency_max[0][0],np.abs(yf[10:N//2])[wh], '+', color='dodgerblue', markersize=5 )\n",
    "ax[0].set_xlabel(r'Frequency [day$^{-1}$]')\n",
    "ax[0].set_ylabel(r'FFT Amplitude')\n",
    "\n",
    "print('Period from FFT1 = '+str(Period)+' days')\n",
    "\n",
    "#compare with shorter time serie:\n",
    "time2 = time[:130]\n",
    "flux2 = flux[:130]\n",
    "# Number of sample points\n",
    "N = len(flux2)\n",
    "# sample spacing\n",
    "T = time_res\n",
    "\n",
    "yf = np.fft.fft(flux2)\n",
    "xf = np.fft.fftfreq(N, T)[1:N//2]\n",
    "\n",
    "ax[1].set_yscale('log')\n",
    "ax[1].plot(xf, np.abs(yf[1:N//2]), color='dodgerblue')\n",
    "\n",
    "wh = np.argwhere((np.abs(yf[1:N//2])==np.max(np.abs(yf[1:N//2]))) )\n",
    "frequency_max2 = xf[:][wh]\n",
    "Period2 = 1/frequency_max2[0][0]\n",
    "ax[1].plot(frequency_max2[0][0],np.abs(yf[1:N//2])[wh], '+', color='dodgerblue', markersize=5 )\n",
    "ax[1].set_xlabel(r'Frequency [day$^{-1}$]')\n",
    "ax[1].set_ylabel(r'FFT Amplitude')\n",
    "\n",
    "print('Period from FFT2 = '+str(Period2)+' days')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3c91668b",
   "metadata": {},
   "outputs": [],
   "source": [
    "newest_flux[np.where(np.isnan(newest_flux))] = 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d3ab99eb",
   "metadata": {},
   "outputs": [],
   "source": [
    "np.sum(np.isnan(new_flux))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "929e78b3",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.9.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}