SlidingWindow1-Basics.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Geometry of Sliding Window Embeddings\n",
    "In this module, you will interactively explore how various parameters of sliding window embeddings of 1 dimensional periodic time series data affect the geometry of their embeddings.  In each code box, press *ctrl-enter* (or *cmd-enter* on a Mac) to execute the code. Progress through the lab sequentially. As you examine the plots in each experiment, answer the questions that follow.\n",
    "\n",
    "This first box imports all of the necessary Python packages to run the code in this module. There will be a similar box at the beginning of every module in this series."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "#Do all of the imports and setup inline plotting\n",
    "import numpy as np\n",
    "from ripser import ripser\n",
    "\n",
    "%matplotlib notebook\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib import gridspec\n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "\n",
    "from sklearn.decomposition import PCA\n",
    "\n",
    "from scipy.interpolate import InterpolatedUnivariateSpline\n",
    "\n",
    "import ipywidgets as widgets\n",
    "from IPython.display import display\n",
    "\n",
    "import warnings\n",
    "warnings.filterwarnings('ignore')\n",
    "\n",
    "from IPython.display import clear_output"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Pure Sinusoid Sliding Window\n",
    "In this first experiment, you will alter the extent of the sliding window of a pure sinusoid and examine how the geometry of a 2-D embedding changes.  \n",
    "\n",
    "First, setup and plot a pure sinusoid in NumPy:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# Step 1: Setup the signal\n",
    "T = 40 # The period in number of samples\n",
    "NPeriods = 4 # How many periods to go through\n",
    "N = T*NPeriods #The total number of samples\n",
    "t = np.linspace(0, 2*np.pi*NPeriods, N+1)[:N] # Sampling indices in time\n",
    "x = np.cos(t) # The final signal\n",
    "plt.plot(x);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Sliding Window Code\n",
    "The code below performs a sliding window embedding on a 1D signal.  The parameters are as follows:\n",
    "\n",
    "| | |\n",
    "|:-:|---|\n",
    "|$x$   | The 1-D signal (numpy array)   |\n",
    "|dim|The dimension of the embedding|\n",
    "|$\\tau$   | The skip between samples in a given window  |\n",
    "|$dT$   | The distance to slide from one window to the next  |\n",
    "\n",
    "That is, along the signal given by the array $x$, the first three windows will will be $$\\begin{bmatrix} x(\\tau)\\\\ x(2\\tau) \\\\ \\ldots \\\\ x((\\mbox{dim}-1)\\cdot\\tau)\\end{bmatrix},  \\begin{bmatrix} x(dT + \\tau)\\\\ x(dT +2\\tau) \\\\ \\ldots \\\\ x(dT +(\\mbox{dim}-1)\\cdot\\tau)\\end{bmatrix},  \\begin{bmatrix} x(2dT + \\tau)\\\\ x(2dT +2\\tau) \\\\ \\ldots \\\\ x(2dT +(\\mbox{dim}-1)\\cdot\\tau)\\end{bmatrix}$$\n",
    "\n",
    "Spline interpolation is used to fill in information between signal samples, which is necessary for certain combinations of parameters, such as a non-integer $\\tau$ or $dT$.\n",
    "\n",
    "The function *getSlidingWindow* below creates an array $X$ containing the windows as its columns."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "def getSlidingWindow(x, dim, Tau, dT):\n",
    "    \"\"\"\n",
    "    Return a sliding window of a time series,\n",
    "    using arbitrary sampling.  Use linear interpolation\n",
    "    to fill in values in windows not on the original grid\n",
    "    Parameters\n",
    "    ----------\n",
    "    x: ndarray(N)\n",
    "        The original time series\n",
    "    dim: int\n",
    "        Dimension of sliding window (number of lags+1)\n",
    "    Tau: float\n",
    "        Length between lags, in units of time series\n",
    "    dT: float\n",
    "        Length between windows, in units of time series\n",
    "    Returns\n",
    "    -------\n",
    "    X: ndarray(N, dim)\n",
    "        All sliding windows stacked up\n",
    "    \"\"\"\n",
    "    N = len(x)\n",
    "    NWindows = int(np.floor((N-dim*Tau)/dT))\n",
    "    if NWindows <= 0:\n",
    "        print(\"Error: Tau too large for signal extent\")\n",
    "        return np.zeros((3, dim))\n",
    "    X = np.zeros((NWindows, dim))\n",
    "    spl = InterpolatedUnivariateSpline(np.arange(N), x)\n",
    "    for i in range(NWindows):\n",
    "        idxx = dT*i + Tau*np.arange(dim)\n",
    "        start = int(np.floor(idxx[0]))\n",
    "        end = int(np.ceil(idxx[-1]))+2\n",
    "        # Only take windows that are within range\n",
    "        if end >= len(x):\n",
    "            X = X[0:i, :]\n",
    "            break\n",
    "        X[i, :] = spl(idxx)\n",
    "    return X"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Sliding Window Result\n",
    "\n",
    "We will now perform a sliding window embedding with various choices of parameters.  [Principal component analysis](https://en.wikipedia.org/wiki/Principal_component_analysis) will be performed to project the result down to 2D for visualization.  \n",
    "\n",
    "The first two eigenvalues computed by PCA will be printed.  The closer these eigenvalues are to each other, the rounder and more close to a circle the 2D projection of the embedding is. A red vertical line will be drawn to show the product of $\\tau$ and the dimension, or \"extent\" (window length).\n",
    "\n",
    "**An important note:** we choose to project the results to 2D (or later, to 3D). Nothing in particular tells us that this is the best choice of dimension. We merely make this choice to enable visualization. In general, when doing PCA, we want to choose enough eigenvalues to account for a significant portion of [explained variance](https://stats.stackexchange.com/questions/22569/pca-and-proportion-of-variance-explained).\n",
    "\n",
    "** Exercise: ** Execute the code. Using the sliders, play around with the parameters of the sliding window embedding and examine the results.  Then answer the questions below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "def on_value_change(change):\n",
    "    execute_computation1()\n",
    "    \n",
    "dimslider = widgets.IntSlider(min=1,max=40,value=20,description='Dimension:',continuous_update=False)\n",
    "dimslider.observe(on_value_change, names='value')\n",
    "\n",
    "Tauslider = widgets.FloatSlider(min=0.1,max=5,step=0.1,value=1,description=r'\\(\\tau :\\)' ,continuous_update=False)\n",
    "Tauslider.observe(on_value_change, names='value')\n",
    "\n",
    "dTslider = widgets.FloatSlider(min=0.1,max=5,step=0.1,value=0.5,description='dT: ',continuous_update=False)\n",
    "dTslider.observe(on_value_change, names='value')\n",
    "\n",
    "display(widgets.HBox(( dimslider,Tauslider,dTslider)))\n",
    "\n",
    "plt.figure(figsize=(9.5, 3))\n",
    "def execute_computation1():    \n",
    "    plt.clf()\n",
    "    # Step 1: Setup the signal again in case x was lost\n",
    "    T = 40 # The period in number of samples\n",
    "    NPeriods = 4 # How many periods to go through\n",
    "    N = T*NPeriods # The total number of samples\n",
    "    t = np.linspace(0, 2*np.pi*NPeriods, N+1)[0:N] # Sampling indices in time\n",
    "    x = np.cos(t)  # The final signal\n",
    "    \n",
    "    # Get slider values\n",
    "    dim = dimslider.value\n",
    "    Tau = Tauslider.value\n",
    "    dT = dTslider.value\n",
    "    \n",
    "    #Step 2: Do a sliding window embedding\n",
    "    X = getSlidingWindow(x, dim, Tau, dT)\n",
    "    extent = Tau*dim\n",
    "\n",
    "    #Step 3: Perform PCA down to 2D for visualization\n",
    "    pca = PCA(n_components = 2)\n",
    "    Y = pca.fit_transform(X)\n",
    "    eigs = pca.explained_variance_\n",
    "    print(\"lambda1 = %g, lambda2 = %g\"%(eigs[0], eigs[1]))\n",
    "\n",
    "    #Step 4: Plot original signal and PCA of the embedding\n",
    "    ax = plt.subplot(121)\n",
    "    ax.plot(x)\n",
    "    ax.set_ylim((-2*max(x), 2*max(x)))\n",
    "    ax.set_title(\"Original Signal\")\n",
    "    ax.set_xlabel(\"Sample Number\")\n",
    "    yr = np.max(x)-np.min(x)\n",
    "    yr = [np.min(x)-0.1*yr, np.max(x)+0.1*yr]\n",
    "    ax.plot([extent, extent], yr, 'r')\n",
    "    ax.plot([0, 0], yr, 'r')     \n",
    "    ax.plot([0, extent], [yr[0]]*2, 'r')\n",
    "    ax.plot([0, extent], [yr[1]]*2, 'r')\n",
    "    ax2 = plt.subplot(122)\n",
    "    ax2.set_title(\"PCA of Sliding Window Embedding\")\n",
    "    ax2.scatter(Y[:, 0], Y[:, 1])\n",
    "    ax2.set_aspect('equal', 'datalim')\n",
    "    plt.tight_layout()\n",
    "    \n",
    "execute_computation1()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Questions\n",
    "1. For fixed $\\tau$:\n",
    "  1. What does varying the dimension do to the *extent* (the length of the window)?\n",
    "  1. what dimensions give eigenvalues nearest each other? (Note: dimensions! Plural!) Explain why this is the case. Explain how you might use this information to deduce the period of a signal.\n",
    "<br><br>\n",
    "1. What does varying $dT$ do to the PCA embedding? Explain this in terms of the definition of sliding windows above.\n",
    "<br><br>\n",
    "1. The command \n",
    "```python\n",
    "np.random.randn(pts)\n",
    "```\n",
    "generates an array of length *pts* filled with random values drawn from a standard normal distribution ($\\mu=0$, $\\sigma=1$). Modify the code above to add random noise to signal. \n",
    "  1. Can you still detect the period visually by inspecting the plot of the signal?\n",
    "  1. Does your method of detecting the period from the first question still work?\n",
    "  1. How does adding noise change the geometry of the PCA embedding? \n",
    "  1. Modify the amplitude of the noise (for example, by multiplying the noise-generating command by a constant) and examine the 2D projection. What feature of the 2D projection appears to imply that the signal is periodic? At what noise amplitude does this feature appear to vanish?\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Non-Periodic Signal Sliding Window\n",
    "For a contrasting example, we will now examine the sliding window embedding of a non-periodic signal which is a linear function plus Gaussian noise.  The code below sets up the signal and then does the sliding window embedding, as before.  \n",
    "\n",
    "** Exercise: ** Execute the code. Using the sliders, play around with the parameters of the sliding window embedding and examine the results.  Then answer the questions below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "noise = 0.05*np.random.randn(400)\n",
    "\n",
    "def on_value_change(change):\n",
    "    execute_computation2()\n",
    "    \n",
    "dimslider = widgets.IntSlider(min=1,max=40,value=20,description='Dimension:',continuous_update=False)\n",
    "dimslider.observe(on_value_change, names='value')\n",
    "\n",
    "Tauslider = widgets.FloatSlider(min=0.1,max=5,step=0.1,value=1,description='Tau: ',continuous_update=False)\n",
    "Tauslider.observe(on_value_change, names='value')\n",
    "\n",
    "dTslider = widgets.FloatSlider(min=0.1,max=5,step=0.1,value=0.5,description='dT: ',continuous_update=False)\n",
    "dTslider.observe(on_value_change, names='value')\n",
    "\n",
    "display(widgets.HBox(( dimslider,Tauslider,dTslider)))\n",
    "\n",
    "plt.figure(figsize=(9.5, 3))\n",
    "\n",
    "def execute_computation2():    \n",
    "    plt.clf()\n",
    "    # Step 1: Set up the signal\n",
    "    x = np.arange(400)\n",
    "    x = x/float(len(x))\n",
    "    x = x + noise # Add some noise\n",
    "    \n",
    "    # Get slider values\n",
    "    dim = dimslider.value\n",
    "    Tau = Tauslider.value\n",
    "    dT = dTslider.value\n",
    "    \n",
    "    #Step 2: Do a sliding window embedding\n",
    "    X = getSlidingWindow(x, dim, Tau, dT)\n",
    "    extent = Tau*dim\n",
    "\n",
    "    #Step 3: Perform PCA down to 2D for visualization\n",
    "    pca = PCA(n_components = 2)\n",
    "    Y = pca.fit_transform(X)\n",
    "    eigs = pca.explained_variance_\n",
    "    print(\"lambda1 = %g, lambda2 = %g\"%(eigs[0], eigs[1]))\n",
    "\n",
    "    #Step 4: Plot original signal and PCA of the embedding\n",
    "    gs = gridspec.GridSpec(1, 2)\n",
    "    ax = plt.subplot(gs[0,0])\n",
    "    ax.plot(x)\n",
    "    ax.set_ylim((-2, 2))\n",
    "    ax.set_title(\"Original Signal\")\n",
    "    ax.set_xlabel(\"Sample Number\")\n",
    "    yr = np.max(x)-np.min(x)\n",
    "    yr = [np.min(x)-0.1*yr, np.max(x)+0.1*yr]\n",
    "    ax.plot([extent, extent], yr, 'r')\n",
    "    ax.plot([0, 0], yr, 'r')     \n",
    "    ax.plot([0, extent], [yr[0]]*2, 'r')\n",
    "    ax.plot([0, extent], [yr[1]]*2, 'r')       \n",
    "    ax2 = plt.subplot(gs[0, 1])\n",
    "    ax2.set_title(\"PCA of Sliding Window Embedding\")\n",
    "    ax2.scatter(Y[:, 0], Y[:, 1])\n",
    "    ax2.set_aspect('equal', 'datalim')\n",
    "\n",
    "execute_computation2()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Questions\n",
    "\n",
    "1. Notice how changing the window extent doesn't have the same impact as it did in the periodic example above.  Why might this be?\n",
    "<br><br>\n",
    "1. Why is the second eigenvalue always tiny?\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## Multiple Sines Sliding Window\n",
    "We will now go back to periodic signals, but this time we will increase the complexity by adding two sines together.  If the ratio between the two sinusoids is a rational number, then they are called *harmonics* of each other.  For example, $\\sin(t)$ and $\\sin(3t)$ are harmonics of each other. By contrast, if the ratio of the two frequencies is irrational, then the sinusoids are called *incommensurate*.  For example, $\\sin(t)$ and $\\sin(\\pi t)$ are incommensurate.\n",
    "\n",
    "The plots below are initialized with \n",
    "\n",
    "$$f(t) = \\sin(\\omega t) + \\sin(3\\omega t),$$\n",
    "\n",
    "a sum of two harmonics.  \n",
    "\n",
    "This time, the eigenvalues of PCA will be plotted (up to the first 10), in addition to the red line showing the extent of the window.  Also, 3D PCA will be displayed instead of 2D PCA, and you can click and drag your mouse to view it from different angles.  Colors will be drawn to indicate the position of the window in time, with cool colors (greens and yellows) indicating earlier windows and hot colors (oranges and reds) indicating later windows (using the \"jet\" colormap).\n",
    "\n",
    "** Exercise: ** Execute the code. Then play with the sliders, as well as the embedding dimension (note that for the 3-D projection, you can change the view by dragging around. Do!). Then, try changing the second sinusoid to be another multiple of the first. Try both harmonic and incommensurate values. Once you have gotten a feel for the geometries and the eigenvalues, answer the questions below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "def on_value_change(change):\n",
    "    execute_computation3()\n",
    "\n",
    "embeddingdimbox = widgets.Dropdown(options=[2, 3],value=3,description='Embedding Dimension:',disabled=False)\n",
    "embeddingdimbox.observe(on_value_change,names='value')\n",
    "\n",
    "secondfreq = widgets.Dropdown(options=[2, 3, np.pi],value=3,description='Second Frequency:',disabled=False)\n",
    "secondfreq.observe(on_value_change,names='value')\n",
    "\n",
    "noiseampslider = widgets.FloatSlider(min=0,max=6,step=0.5,value=0,description='Noise Amplitude',continuous_update=False)\n",
    "noiseampslider.observe(on_value_change, names='value')\n",
    "\n",
    "dimslider = widgets.IntSlider(min=1,max=100,value=30,description='Dimension:',continuous_update=False)\n",
    "dimslider.observe(on_value_change, names='value')\n",
    "\n",
    "Tauslider = widgets.FloatSlider(min=0.1,max=5,step=0.1,value=1,description='Tau: ',continuous_update=False)\n",
    "Tauslider.observe(on_value_change, names='value')\n",
    "\n",
    "dTslider = widgets.FloatSlider(min=0.1,max=5,step=0.1,value=0.5,description='dT: ',continuous_update=False)\n",
    "dTslider.observe(on_value_change, names='value')\n",
    "\n",
    "display(widgets.HBox(( secondfreq,embeddingdimbox,noiseampslider)))\n",
    "display(widgets.HBox((dimslider,Tauslider,dTslider)))\n",
    "\n",
    "noise = np.random.randn(10000)\n",
    "\n",
    "fig = plt.figure(figsize=(9.5, 4))\n",
    "def execute_computation3():\n",
    "    plt.clf()\n",
    "    \n",
    "    # Step 1: Setup the signal\n",
    "    T1 = 20 # The period of the first sine in number of samples\n",
    "    T2 = T1*secondfreq.value # The period of the second sine in number of samples\n",
    "    NPeriods = 10 # How many periods to go through, relative to the second sinusoid\n",
    "    N = T2*NPeriods # The total number of samples\n",
    "    t = np.arange(N) # Time indices\n",
    "    x = np.cos(2*np.pi*(1.0/T1)*t) # The first sinusoid\n",
    "    x += np.cos(2*np.pi*(1.0/T2)*t) # Add the second sinusoid\n",
    "    x += noiseampslider.value*noise[:len(x)]\n",
    "    \n",
    "    # Get widget values\n",
    "    dim = dimslider.value\n",
    "    Tau = Tauslider.value\n",
    "    dT = dTslider.value\n",
    "    embeddingdim = embeddingdimbox.value\n",
    "    \n",
    "    # Step 2: Do a sliding window embedding\n",
    "    X = getSlidingWindow(x, dim, Tau, dT)\n",
    "    extent = Tau*dim\n",
    "\n",
    "    # Step 3: Perform PCA down to dimension chosen for visualization\n",
    "    pca = PCA(n_components = 10)\n",
    "    Y = pca.fit_transform(X)\n",
    "    eigs = pca.explained_variance_\n",
    "\n",
    "    # Step 4: Plot original signal and PCA of the embedding \n",
    "    gs = gridspec.GridSpec(2, 2,width_ratios=[1, 2])\n",
    "    \n",
    "    # Plot the signal\n",
    "    ax = plt.subplot(gs[0,0])\n",
    "    ax.plot(x)\n",
    "    yr = np.max(x)-np.min(x)\n",
    "    yr = [np.min(x)-0.1*yr, np.max(x)+0.1*yr]\n",
    "    ax.plot([extent, extent], yr, 'r')\n",
    "    ax.plot([0, 0], yr, 'r')     \n",
    "    ax.plot([0, extent], [yr[0]]*2, 'r')\n",
    "    ax.plot([0, extent], [yr[1]]*2, 'r')\n",
    "    ax.set_title(\"Original Signal\")\n",
    "    ax.set_xlabel(\"Sample Number\")\n",
    "\n",
    "    c = plt.get_cmap('jet')\n",
    "    C = c(np.array(np.round(np.linspace(0, 255, Y.shape[0])), dtype=np.int32))\n",
    "    C = C[:, 0:3]\n",
    "\n",
    "    # Plot the PCA embedding\n",
    "    if embeddingdim == 3:\n",
    "        ax2 = plt.subplot(gs[:,1],projection='3d')\n",
    "        ax2.scatter(Y[:, 0], Y[:, 1], Y[:, 2], c=C)\n",
    "        ax2.set_aspect('equal', 'datalim')\n",
    "    else:\n",
    "        ax2 = plt.subplot(gs[:,1])\n",
    "        ax2.scatter(Y[:, 0], Y[:, 1],c=C)\n",
    "\n",
    "    ax2.set_title(\"PCA of Sliding Window Embedding\")\n",
    "    ax2.set_aspect('equal', 'datalim')\n",
    "\n",
    "    # Plot the eigenvalues as bars\n",
    "    ax3 = plt.subplot(gs[1,0])\n",
    "    eigs = eigs[0:min(len(eigs), 10)]\n",
    "    ax3.bar(np.arange(len(eigs)), eigs)\n",
    "    ax3.set_xlabel(\"Eigenvalue Number\")\n",
    "    ax3.set_ylabel(\"Eigenvalue\")\n",
    "    ax3.set_title(\"PCA Eigenvalues\")\n",
    "\n",
    "    plt.tight_layout()\n",
    "\n",
    "    plt.show();\n",
    "\n",
    "execute_computation3()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Questions\n",
    "\n",
    "1. Comment on the relationship between the eigenvalues and the extent (width) of the window.\n",
    "<br><br>\n",
    "1. When are the eigenvalues near each other? When are they not?\n",
    "<br><br>\n",
    "1. Comment on the change in geometry when the second sinusoid is incommensurate to the first. Specifically, comment on the *intrinsic dimension* of the object in the projection. Can you name the shape in the 3-D projection in the incommensurate case?\n",
    "<br><br>\n",
    "1. Try adding noise in like you did in the single frequency case. \n",
    "  1. Can you distinguish in the projection between the incommensurate case and the noisy, but harmonic one with second frequency 3? Explain.\n",
    "  1. What can you say about the eigenvalues in the two cases? Explain your answer.\n",
    "\n",
    "It seems reasonable to ask what the ideal dimension to embed into is. While that question may be answerable, it would be better to bypass the question altogether. Similarly, it seems that beyond detecting the largest period, these tools are limited in detecting the secondary ones. Topological tools that we will see beginning in the next lab will allow us to make some progress toward that goal."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Power Spectrum\n",
    "\n",
    "We saw above that for a rather subtle change in frequency changing the second sinusoid from harmonic to noncommensurate, there is a marked change in the geometry.  By contrast, the power spectral density functions are very close between the two, as shown below.  Hence, it appears that geometric tools are more appropriate for telling the difference between these two types of signals"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "T = 20 #The period of the first sine in number of samples\n",
    "NPeriods = 10 #How many periods to go through, relative to the faster sinusoid\n",
    "N = T*NPeriods*3 #The total number of samples\n",
    "t = np.arange(N) #Time indices\n",
    "\n",
    "#Make the harmonic signal cos(t) + cos(3t)\n",
    "xH = np.cos(2*np.pi*(1.0/T)*t) + np.cos(2*np.pi*(1.0/(3*T)*t))\n",
    "                                        \n",
    "#Make the incommensurate signal cos(t) + cos(pi*t)\n",
    "xNC = np.cos(2*np.pi*(1.0/T)*t) + np.cos(2*np.pi*(1.0/(np.pi*T)*t))\n",
    "\n",
    "plt.figure()\n",
    "P1 = np.abs(np.fft.fft(xH))**2\n",
    "P2 = np.abs(np.fft.fft(xNC))**2\n",
    "plt.plot(np.arange(len(P1)), P1)\n",
    "plt.plot(np.arange(len(P2)), P2)\n",
    "plt.xlabel(\"Frequency Index\")\n",
    "plt.legend({\"Harmonic\", \"Noncommensurate\"})\n",
    "plt.xlim([0, 50])\n",
    "plt.show();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "# Summary\n",
    "\n",
    "* Signals can be transformed into geometric objects via embeddings.\n",
    "<br><br>\n",
    "* Signal properties are captured by the geometry of the sliding window embedding.  Periodicity corresponds to circularity, period length over window size corresponds to roundness, number of incommensurate frequencies corresponds to intrinsic dimension.\n",
    "<br><br>\n",
    "* The window extent is one of the most important parameters for determining roundness.\n",
    "<br><br>\n",
    "* Adding noise makes things a little trickier to see what's going on by inspection of projections."
   ]
  }
 ],
 "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.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}