Lesson_4.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Fourier Transforms\n",
    "\n",
    "A Fourier transform is a mathematical transformation that takes an existing function and decomposes it into an infinite series of sine or cosine waves of different periods and phases.  It has many uses in science (and in everyday life, as we'll see).  This tutorial will not discuss much of the mathematics behind the transform, but instead today we will concentrate on getting a better intuitive understanding of what Fourier transforms represent, and some ways they can be used. \n",
    "\n",
    "## Definitions\n",
    "\n",
    "First, it's useful to define a few terms.  \n",
    "\n",
    "### Continuous functions (e.g. mathematical functions)\n",
    "\n",
    "Any function can be decomposed into a *Fourier series*, an infinite sum of sines and cosines of different frequencies.  (Another way to say this is that the set of sines and cosines of all frequencies is a complete basis set.) The coefficients of this infinite series for a given function are called the *Fourier coefficients* of that function.   The mathematical operation of determining these coefficients for a particular function is called a *Fourier transform*, which we'll refer to as FT for short. The *inverse Fourier transform* is (as the name implies) the inverse process, which when applied to the Fourier series, recovers the original function.  \n",
    "\n",
    "###  Discrete data points (samples of functions at particular points)\n",
    "\n",
    "When we're dealing with experimental data, we generally do not have continuous sampling of our function of interest.  For example, we might have time samples from observing a particular phenomenon at regular intervals.  This is also the case with digital representations of continuous waveforms, for example a digital audio file containing a certain number of samples of the audio waveform each second.  In this case, we need a slightly different mathematical operation, one that can operate on a discrete set of datapoints rather than on a function.  This operation is called a *discrete Fourier transform*, or DFT.  This operation can be computationally intensive, especially for large datasets, so people have settled upon a particularly efficient implementation of the discrete Fourier transform, called the *fast Fourier transform*, or FFT.  This is so common that \"FFT\" is sometimes used as being synonymous with \"discrete Fourier transform\". \n",
    "\n",
    "###  Time domain vs. Frequency domain\n",
    "\n",
    "One important aspect of Fourier transforms is that the transform and inverse transform take us back and forth between different \"spaces\" or \"domains\".  This is easiest to see by considering an example.  If we have a dataset that is a function of time, then each point in our dataset represents the data at one particular time, e.g. 1 second, 2 seconds, etc.  When we Fourier transform that dataset, then each point in our new dataset represents the behavior of the data at one particular *frequency*, e.g. 1 Hz, 2 Hz, etc.  Note here that frequency is the inverse of time.  And in general it is true that each point in the FT has the inverse units of each point in the original data.  If our dataset represents spatial samples (e.g. elevation as a function of position, in units of meters), then the FT points will have units of 1/meters, i.e. they are *spatial* frequency.   The idea of spatial frequency might seem odd, but it is exactly analogous to temporal frequency.  Each point in that spatial FT represents the coefficient (amplitude and phase) of a sinusoidal wave with a particular wavelength, and how big or small it would have to be to add up to some particular rolling landscape.  Quickly varying terrain would need high (large) spatial frequencies, while slowly varying terrain would need low (small) spatial frequencies. \n",
    "\n",
    "Finally, implicit in this domain-switching idea is a key concept of Fourier analysis, and something that gives it much of its power:  any *one* point in one domain affects the *whole dataset* in the other domain.  This is perhaps easiest to see in going from time to frequency.  If I have a signal that is a pure sine wave, and I Fourier transform it, I only get one non-zero Fourier component, at whatever frequency point is need to represent that signal.  If I now change one other point in the Fourier transform(i.e. in the frequency domain) to be non-zero, and then transform back to the time domain, I have changed the whole signal, because I have added a second sinusoidal component, which will show up at *all* time points.   We'll experiment with this below. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Exercise 1\n",
    "\n",
    "In order to test out what Fourier transforms (FT) can do, we'll want to generate some fake datasets.  Sine waves are easy to deal with since we know what their FTs should look like. \n",
    "\n",
    "Use `numpy`'s `linspace` function to generate a set of time values that are evenly sampled between 0 and 10 seconds.  To start, sample this interval with 100 points.  Set the keyword \"endpoint=False\" to avoid duplicating the zero point at the end.  Calculate the sine of these time points, using a sine wave with a frequency of 0.5 Hz.  You may find it useful to recall that a time-dependent sine wave in terms of frequency is $\\sin(2\\pi t f)$.  Plot the resulting sine wave. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Exercise 2a\n",
    "\n",
    "Now let's Fourier transform our data. `Numpy` has a package called `fft`, which implements various version of the FFT (fast Fourier transform).  Use `np.fft.rfft` to compute the Fourier transform of your sine curve.  (You can look up the exact calling sequence on-line.)   The `r` before `fft` there refers to the fact that we are only dealing with real data here.\n",
    "\n",
    "If you examine the result, you will find that each component of the Fourier transform is a complex number, not a real number.  You can think of the real and imaginary parts of each component as *amplitude* and the *phase* of the sine wave at that frequency, or alternatively as the amplitudes of the sine and cosine terms at that frequency.  (You should be able to convince yourself that these amount to the same thing.) \n",
    "\n",
    "To get the frequencies corresponding to each point in the FT, use `np.fft.rfftfreq`.  You need to pass it the number of points you used ($N$), and the step size ($\\Delta t$). \n",
    "\n",
    "$N$ and $\\Delta t$ are the only parameters you need to figure out which frequencies correspond to the amplitude information in the FFT.  To see this, calculate the frequencies corresponding to your original signal, then calculate the frequencies corresponding to:\n",
    "* a sample of 50 points covering the same 0 to 10 seconds\n",
    "* a sample of 50 points covering only 0 to 5 seconds\n",
    "\n",
    "Remember to set \"endpoint=False.\"\n",
    "\n",
    "Calculate the maximum frequency and the frequency resolution $\\Delta f$ in each case, and verify the following:\n",
    "* the maximum frequency is also $f_{samp}/2$, where $f_{samp} = 1/\\Delta t$\n",
    "* $\\Delta f$ is given by $1/T$ where $T$ is the total time of the record"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "#### Exercise 2b\n",
    "\n",
    "Now look at the amplitude data given by `np.fft.rfft`:\n",
    "\n",
    "* On one plot, show the real and imaginary parts of the FT, plotted as a function of frequency.  Label each axis appropriately.\n",
    "\n",
    "* On a second plot, show the *power spectrum* of the FT, as a function of frequency.  The power spectrum is the squared amplitudes of the FT components, and is a useful way to see which Fourier components (i.e. which frequencies) are strongest.  Recall that for a complex number $a + bi$, the modulus is $\\sqrt{a^2 + b^2}$.  We want the square of that, which you can calculate for a complex number $x$ by using `np.abs(x)**2`.\n",
    "\n",
    "Make sure that the plots you get (especially the power spectrum) make sense to you, given the input signal you used. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "In preparation for the next step, we need to note one coding tip.  We're going to be messing with the FT array, so it will be a good idea to work on a *copy* of the array, so that we don't inadvertently change multiple different points and thus get unexpected results.  You can get an independent copy of a `numpy` array using `.copy()`, like this:\n",
    "\n",
    "`x = y.copy()`\n",
    "\n",
    "Below is a bit more on how this works:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "a = np.array([1,2,3])\n",
    "b = a.copy()  # Now b can be changed without affecting a. \n",
    "b[1] = 100\n",
    "print(\"Working with a copy, b = a.copy()\")\n",
    "print(\"a is\", a)\n",
    "print(\"b is\", b)\n",
    "\n",
    "# Note that this does *not* make b independent:\n",
    "a = np.array([1,2,3])\n",
    "b = a      \n",
    "#  Here b is just a \"view\" of a, so changing b changes a:\n",
    "b[1] = 100\n",
    "print(\"Working with view, b = a\")\n",
    "print(\"a is\", a)\n",
    "print(\"b is\", b)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Exercise 3\n",
    "\n",
    "Now let's see how a tweak in the frequency domain transforms back to a different signal in the time domain. \n",
    "\n",
    "* Make a copy of the FT array to work with. \n",
    "* Pick one point in this FT array, and change its amplitude to be comparable to that of the largest peak.  \n",
    "* Re-plot the power spectrum to make sure it looks like you expect. \n",
    "* Transform that FT array back into the time domain using `np.fft.irfft`, and plot the resulting time series. \n",
    "* Explain why the result you get is what you should expect. \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###  Uses of Fourier transforms: finding periodicities in time series data\n",
    "\n",
    "As we've seen, the power spectrum contains very clear information about what periodic signals are present in a dataset.  So far, we've only used very clean sine waves to see that, but it turns out that the Fourier transform is quite robust at picking out periodic signals even if they are heavily obscured by noise in a dataset (as long as the noise is random, and not periodic itself).  \n",
    "\n",
    "#### Exercise 4\n",
    "\n",
    "Using the same array of time values as above, generate a new signal that is the sum of two different sine waves.  Also generate an array of \"noise\" values, with comparable amplitude to the signal.  Since real experimental noise often follows a Gaussian (or normal) distribution, a good way to do this is to generate a set of random numbers using `np.random.normal` (look up the usage online, and note in particular the `scale` parameter.)  Note that adding artificial noise to a signal is another use of random numbers which we studied last week.  This is often a useful way to see how noise might affect an experiment you're planning (or have already carried out.) \n",
    "\n",
    "* Plot your original signal. \n",
    "* Plot the signal with noise added. \n",
    "* Plot the power spectrum of the noisy data - can you find the input frequencies? \n",
    "* Play around with how big you make the noise component, making it big enough that it is hard to see clearly the signal in the original data.  Does the FT power spectrum still pull out the right frequencies? "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Audio filtering and compression\n",
    "\n",
    "Another use of Fourier transforms is in audio signal processing.  They have wide applicability in filtering or modifying signals, and in compressing the files.  In its simplest form, a digital audio file is just a series of time samples of the audio waveform, very analogous to the simple time series we have been working with above.  Here we'll read one of these in from a file, storing it as a `numpy` array, and we'll be able to work with just as we have above. \n",
    "\n",
    "For the exercises below, you may find it useful to know that the frequency range of human hearing is approximately 20 Hz to 20 kHz.  Use this as a sanity check on your plots - you may see frequencies somewhat outside this range, but if the frequencies you are plotting are wildly different, you have probably made a mistake somewhere. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Here's some sample code to show you how to read audio into a numpy array:\n",
    "import scipy.io.wavfile as wv\n",
    "\n",
    "audiofile = 'Cavalry_short.wav'\n",
    "# This returns both the data and the sampling rate, which is the \n",
    "# inverse of the time step.  You'll need to pass this to the \n",
    "# np.fft.fftfreq routine to get the proper frequencies back. \n",
    "bitrate, audiodata = wv.read(audiofile)\n",
    "\n",
    "# Since this example clip is only mono, the numpy array only \n",
    "# has one dimension.  If it were a stereo file, there would \n",
    "# be a second dimension of size 2, i.e. there would be two\n",
    "# equal-length rows in the array.\n",
    "print(\"The shape of the audio data array is\", audiodata.shape)\n",
    "\n",
    "\n",
    "print(\"The audio is sampled at a rate of %0.1f kHz.\" % (bitrate/1000.))\n",
    "print(\"The clip is %0.1f seconds long.\" % (audiodata.size/bitrate))\n",
    "\n",
    "# We can also embed an audio player directly in the notebook:\n",
    "import IPython.display as disp\n",
    "disp.display(disp.Audio(data=audiodata, rate=bitrate))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "#### Exercise 5\n",
    "\n",
    "Using the above audio data (or a different .WAV file if you have one handy to play around with), do the following: \n",
    "\n",
    "* Plot the waveform vs. time.\n",
    "* Plot the Fourier transform power spectrum vs. frequency.  Since human hearing spans three decades of frequency, it may be useful to make this plot have a logarithmic x axis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Exercise 6\n",
    "\n",
    "Now let's do something with the audio! \n",
    "\n",
    "* Split the frequency spectrum in half, i.e. figure out a frequency in the middle of the range.  (You might want to think about what frequency makes the most sense here - arithmetic mean?  geometric mean?  something else?)   \n",
    "* Generate a 'mask' to only pass certain frequencies.  Remember, you can create a numpy boolean array with the same dimensions as array `x` by using a statement like `big = x > 1000`; your array `big` will then be True whereever the array `x` had values greater than 1000. \n",
    "* Using your mask array, create a new copy of the FT with all of the high frequency components set to zero.  Inverse transform the FT, and play the resulting audio.  How does it sound? \n",
    "* Do the same thing, but this time zero out the lows, and play the result. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Exercise 7\n",
    "\n",
    "Often it is useful to be able to compress a signal, so that, e.g., all those music files don't take up a lot of space on your computer or your phone. The above techniques are one form of compression - they reduce the amount of data we're dealing by a factor of two!  (At least they would if we only stored the Fourier transform in the file - you can check for yourself that the inverse-transformed audio file is the same size as before.)  But the resulting sound is pretty bad. Obviously, if possible we'd prefer to do something that doesn't distort the sound so much, but still reduces the file size. \n",
    "\n",
    "One very simple way to do this is to look at the FT power spectrum, and throw away components with very low power.  By definition, these contribute less to the final overall signal, but it costs us just as much storage space to store a floating-point number with a low value as to store one with a high value.  Instead, we just set the low-power components to zero, and store the FT amplitudes and frequencies.  The software that reconstructs the sound for playback can be written to assume that any frequency not present in the FT has an amplitude of zero, and to fill in the grid before doing the inverse transform.\n",
    "\n",
    "Try it:\n",
    "\n",
    "* Find the entries in the FT power spectrum that have the lowest amplitudes, and set them all to zero.  Set your threshold so that you are discarding half the datapoints. \n",
    "* Inverse transform, and play back the resulting audio.  Can you hear the difference? \n",
    "* Try it again, for other cutoff values (i.e. throwing out more or less data).  What percentage of datapoints do you need to discard before you start to hear a difference compared to the original audio? "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Higher-dimensional Fourier transforms\n",
    "\n",
    "Finally, we can generalize the Fourier transform to higher dimensions.  In particular, a 2-D Fourier transform of image data can be a useful way of analyzing the image for spatial information, especially information about repeating patterns.  For example, the plots you often see of the cosmic microwave background (CMB) anisotropies show a power spectrum (bottom) that results from Fourier transforming the sky brightness distribution of hot and cool spots (top): <p><img src=\"http://wmap.gsfc.nasa.gov/media/121238/ilc_9yr_moll720.jpg\" style=\"width: 60%\"> <br /><img src=\"http://wmap.gsfc.nasa.gov/media/060911/PowerSpectrum150.png\" style=\"width: 52%\"></p>  (Image credit: NASA/WMAP)\n",
    "\n",
    "(Although the above power spectrum is only one-dimensional, that is because it is essentially the azimuthal average of symmetric, two-dimensional power spectrum.) \n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Exercise 8\n",
    "\n",
    "Generate an array of zeros that is 100x100 in size. (`np.zeros` will do this for you.)  If we think of this as a Fourier transform, it is an FT of a completely empty image, i.e. every Fourier component has a value of zero. \n",
    "\n",
    "Now change *one* pixel of that array to 1 instead of 0, and inverse-transform it.  (You can use the tasks `np.fft.fft2` and `np.fft.irfft2` to do the forward and inverse FTs of 2D data.)  \n",
    "\n",
    "Display the resulting image, but before you do, try to predict what it will look like.  What should changing one Fourier component do?  How might you expect it to look different if you change which pixel you set to 1? How does it change if you change the $x$ coordinate of your key pixel, vs. changing the $y$ coordinate? \n",
    "\n",
    "Make a few representative plots, and write a few sentences explaining what you see. \n",
    "\n",
    "Some useful tasks:  you can display images with `plt.imshow(myimage.real, cmap='gray')`.  You'll find that the returned image data after the inverse transform has some imaginary components.  You can avoid these in plotting your image by just referring to `myimage.real` as above.  (In general you can do this with any Python complex number.) \n",
    "\n",
    "You can get the frequencies associated with each axis using `np.fft.fftfreq(n)`.  Note that the frequency array is one-dimensional, but gives the frequencies along a given axis.  These are the same for both axes as long as the image is square.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Exercise 9\n",
    "\n",
    "One principle of optics that involves Fourier transforms is that the diffraction pattern of light passing through an aperture of a given shape is the Fourier transform of that shape.  For example, when light passes through a circular aperture (e.g. in a telescope), the image of a point source becomes an Airy pattern:\n",
    "\n",
    "<a title=\"Anaqreon at English Wikipedia [CC0], via Wikimedia Commons\" href=\"https://commons.wikimedia.org/wiki/File%3AAiry_disk_created_by_laser_beam_through_pinhole.jpg\"><img  alt=\"Airy disk created by laser beam through pinhole\" src=\"https://upload.wikimedia.org/wikipedia/commons/thumb/4/4b/Airy_disk_created_by_laser_beam_through_pinhole.jpg/512px-Airy_disk_created_by_laser_beam_through_pinhole.jpg\"/></a>\n",
    "\n",
    "Here you will show that this Fourier transform relation is correct, and explore how it depends on other factors, such as aperture size and the presence of other structures. \n",
    "\n",
    "* Use the function below to create images of circular apertures of different radii.  Make a plot of the aperture to see what it looks like. \n",
    "* Take the Fourier transform of your aperture image, and plot its squared amplitude (the 2D version of its power spectrum).  Before plotting, apply `np.fft.fftshift` to your FT so that the lowest frequencies are at the center.  This is a common way to plot 2D FTs, akin to putting the origin 0,0 at the center of an $x, y$ plot. \n",
    "* You may find it useful to plot the logarithm of power spectrum in order to better see the structure. \n",
    "\n",
    "Comment on the affect of aperture size on the resulting diffraction pattern.   If you have time, try turning on the options in the function that imitate the presence of a secondary mirror in a telescope (which obstructs part of the light), and the supports for that structure. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def aperture(primary_r, secondary_r=0, supports=False, n=1024):\n",
    "    '''Return a 2D array representing a telescope aperture. The \n",
    "    only required input is primary_r, the radius of the primary mirror \n",
    "    in pixels.  This must be smaller than half of the grid size n, \n",
    "    which by default is 1024.  Setting secondary_r to a positive value \n",
    "    adds a central obstruction, akin to a telescope secondary mirror, and \n",
    "    setting support=True add mirror support struts. '''\n",
    "    \n",
    "    import numpy as np\n",
    "    \n",
    "    # Check the inputs:\n",
    "    assert primary_r < n/2.0, \"primary_r must be less than n/2\"\n",
    "    \n",
    "    # Make an empty grid:\n",
    "    aperture = np.zeros((n, n))\n",
    "\n",
    "    # Make a circle at the center:\n",
    "    mid = int(n/2.0)\n",
    "\n",
    "    xcenter, ycenter = mid, mid\n",
    "    r = primary_r\n",
    "\n",
    "    y,x = np.ogrid[-xcenter:n-xcenter, -ycenter:n-ycenter]\n",
    "    mask = x**2 + y**2 <= r**2\n",
    "\n",
    "    aperture[mask] = 1\n",
    "\n",
    "    if secondary_r > 0:\n",
    "        assert secondary_r < primary_r, 'secondary must be smaller than primary'\n",
    "        \n",
    "        #Now add a dark spot in the middle, secondary obstruction:\n",
    "        r2 = secondary_r\n",
    "        y,x = np.ogrid[-xcenter:n-xcenter, -ycenter:n-ycenter]\n",
    "        mask2 = x**2 + y**2 <= r2**2\n",
    "        aperture[mask2] = 0\n",
    "\n",
    "    if supports:\n",
    "        # Add mirror supports:\n",
    "        support_width = 5\n",
    "        aperture[:,mid-support_width:mid+support_width] = 0\n",
    "        aperture[mid-support_width:mid+support_width, :] = 0\n",
    "    \n",
    "    return aperture"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python [py36]",
   "language": "python",
   "name": "Python [py36]"
  },
  "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.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}