SlidingWindow3-AudioApplications.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h1>Audio Applications</h1><BR>\n",
    "Now that we have established the theory behind the geometry of 1D time series sliding window embeddings, we will look at our first real applications: audio signals.  First, we will import all of the necessary libraries and define the sliding window code as before.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "##Do all of the imports and setup inline plotting\n",
    "%matplotlib notebook\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from sklearn.decomposition import PCA\n",
    "from scipy.interpolate import InterpolatedUnivariateSpline\n",
    "\n",
    "from ripser import ripser\n",
    "from persim import plot_diagrams\n",
    "import scipy.io.wavfile\n",
    "\n",
    "from IPython.display import clear_output\n",
    "\n",
    "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": [
    "<h1>Biphonation Overview</h1>\n",
    "\n",
    "Biphonation refers to the presence of two or more simultaneous frequencies in a signal which are \"incommensurate\"; that is, their frequencies are linearly independent over the rational numbers.  In other words, the frequencies are \"inharmonic.\"  We saw a synthetic example in class 1 of cos(x) + cos(pi x).  Today, we will examine how this manifests itself in biology with horse whinnies that occur during states of high emotional valence.  During the steady state of a horse whinnie, biphonation is found\n",
    "\n",
    "<table>\n",
    "<tr><td>\n",
    "<img src = \"Whinnie.png\">\n",
    "</td></tr>\n",
    "<tr><td>\n",
    "<b>Figure 1</b>: Audio of a horse whinnie.  Courtesy of <a href = \"http://www.nature.com/articles/srep09989#s1\">http://www.nature.com/articles/srep09989#s1</a></td></tr>\n",
    "</table>\n",
    "\n",
    "<!--<iframe width=\"560\" height=\"315\" src=\"https://www.youtube.com/embed/f8DdGpHkzu4\" frameborder=\"0\" allowfullscreen></iframe>!-->\n",
    "<BR><BR>\n",
    "<h1>Biphonation Example with Horse Whinnies</h1><BR>\n",
    "\n",
    "Let's now load the audio from the horse whinnie example, interactively plot the audio waveform, and listen to it.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Read in the audio file.  Fs is the sample rate, and\n",
    "#X is the audio signal\n",
    "Fs, X = scipy.io.wavfile.read(\"horsewhinnie.wav\")\n",
    "plt.figure()\n",
    "plt.plot(np.arange(len(X))/float(Fs), X)\n",
    "plt.xlabel(\"Time (Seconds)\")\n",
    "plt.title(\"Horse Whinnie Waveform\")\n",
    "plt.show()\n",
    "\n",
    "from IPython.display import Audio\n",
    "# load a remote WAV file\n",
    "Audio('horsewhinnie.wav')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The code below will extract a subsection of the signal and perform a sliding window embedding + 1D persistent homology.  Using the interactive plot of the audio waveform above, find two different time ranges to plot:<BR>\n",
    "\n",
    "<ol>\n",
    "<li>A region with a pure tone (single sinusoid), which can be found towards the beginning</li>\n",
    "<li>A region with biphonation, which can be found towards the middle.  Ensure that this region has two strongly persistent classes with early birth times.  The class will compete to find the region which shows biphonation the most clearly with these statistics, and the score will be based on the <i>second largest persistence</i>, which will be indicated in the persistence diagram plot<BR></li>\n",
    "</ol>\n",
    "To interactively search for regions, use the pan icon\n",
    "<img src = \"PanIcon.png\">\n",
    "then left click and drag to translate, and right click and drag to zoom.  Once you've found a region, modify the \"time\" variable in the code below accordingly, and run it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "#These variables are used to adjust the window size\n",
    "F0 = 493 #First fundamental frequency\n",
    "G0 = 1433 #Second fundamental frequency\n",
    "\n",
    "###TODO: Modify this variable (time in seconds)\n",
    "time = 0.91\n",
    "\n",
    "#Step 1: Extract an audio snippet starting at the chosen time\n",
    "SigLen = 512 #The number of samples to take after the start time\n",
    "iStart = int(round(time*Fs))\n",
    "x = X[iStart:iStart + SigLen]\n",
    "W = int(round(Fs/G0))\n",
    "\n",
    "#Step 2: Get the sliding window embedding\n",
    "Y = getSlidingWindow(x, W, 2, 2)\n",
    "#Mean-center and normalize\n",
    "Y = Y - np.mean(Y, 1)[:, None]\n",
    "Y = Y/np.sqrt(np.sum(Y**2, 1))[:, None]\n",
    "\n",
    "#Step 3: Do the 1D rips filtration\n",
    "PDs = ripser(Y, maxdim=1)['dgms']\n",
    "PD = PDs[1]\n",
    "\n",
    "#Step 4: Figure out the second largest persistence\n",
    "sP = 0\n",
    "sPIdx = 0\n",
    "if PD.shape[0] > 1:\n",
    "    Pers = PD[:, 1] - PD[:, 0]\n",
    "    sPIdx = np.argsort(-Pers)[1]\n",
    "    sP = Pers[sPIdx]\n",
    "    \n",
    "#Step 5: Plot the results\n",
    "plt.figure(figsize=(8, 4))\n",
    "plt.subplot(121)\n",
    "plt.title(\"Starting At %g Seconds\"%time)\n",
    "plt.plot(time + np.arange(SigLen)/Fs, x)\n",
    "plt.xlabel(\"Time\")\n",
    "plt.subplot(122)\n",
    "plot_diagrams(PDs)\n",
    "plt.plot([PD[sPIdx, 0]]*2, PD[sPIdx, :], 'r')\n",
    "plt.scatter(PD[sPIdx, 0], PD[sPIdx, 1], 20, 'r')\n",
    "plt.title(\"Second Largest Persistence: %g\"%sP)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "<h1>Music Analysis</h1>\n",
    "\n",
    "<img src = \"journey.jpg\"><BR><BR>\n",
    "\n",
    "Music is full of repetition.  For instance, there is usually a hierarchy of rhythm which determines how the music \"pulses,\" or repeates itself in beat patterns.  Often, a dominant rhythm level is deemed the \"tempo\" of the music.  Typical tempos range from about 50 beats per minute to 200 beats per minute.  Let's take a moderate tempo level of 120 beats per minute, for instance, which occurs in the song \"Don't Stop Believin'\" by Journey.  This corresponds to a period of 0.5 seconds.  As we saw in the horse example, sound is sampled at 44100 samples per second.  This corresponds to an ideal sliding window interval length of 22050.  Let's try to compute the sliding window embedding of the raw audio to see if the tempo manifests itself with TDA.  First, we will load in \"Don't Stop Believing\" below:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Fs, X = scipy.io.wavfile.read(\"journey.wav\") #Don't Stop Believing\n",
    "X = X/(2.0**15) #Loaded in as 16 bit shorts, convert to float\n",
    "plt.figure()\n",
    "plt.plot(np.arange(len(X))/float(Fs), X)\n",
    "plt.xlabel(\"Time (Seconds)\")\n",
    "plt.title(\"Don't Stop Believin\")\n",
    "plt.show()\n",
    "\n",
    "Audio('journey.wav')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<BR><BR>\n",
    "Now let's do a sliding window with a window length equal to the sample rate over 2, corresponding to the fact that a beat period is a half of a second in this song.  We will have to have a large <code>Tau</code> and <code>dT</code>, since there is such a high sampling rate, because otherwise the TDA code will grind to a halt with way too many points.  We will also need to set up special sliding window code that skips the spline interpolation step, because this will also be prohibitively slow at this sampling rate.  In other words, we will assume that <code>Tau</code> and <code>dT</code> are integers.  Here's the code that does all of this on the first three seconds of audio:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Sliding window code here assumes integer x, dim, and Tau so no interpolation\n",
    "#is needed (for computational efficiency)\n",
    "def getSlidingWindowInteger(x, dim, Tau, dT):\n",
    "    N = len(x)\n",
    "    NWindows = int(np.floor((N-dim*Tau)/dT)) #The number of windows\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)) #Create a 2D array which will store all windows\n",
    "    idx = np.arange(N)\n",
    "    for i in range(NWindows):\n",
    "        #Figure out the indices of the samples in this window\n",
    "        idxx = np.array(dT*i + Tau*np.arange(dim), dtype=np.int32)\n",
    "        X[i, :] = x[idxx]\n",
    "    return X\n",
    "\n",
    "\n",
    "#Note that dim*Tau here spans a half a second of audio, \n",
    "#since Fs is the sample rate\n",
    "dim = round(Fs/200)\n",
    "Tau = 100\n",
    "dT = Fs/100\n",
    "Y = getSlidingWindowInteger(X[0:Fs*3], dim, Tau, dT)\n",
    "print(\"Y.shape = \", Y.shape)\n",
    "#Mean-center and normalize\n",
    "Y = Y - np.mean(Y, 1)[:, None]\n",
    "Y = Y/np.sqrt(np.sum(Y**2, 1))[:, None]\n",
    "\n",
    "PDs = ripser(Y, maxdim=1)['dgms']\n",
    "pca = PCA()\n",
    "Z = pca.fit_transform(Y)\n",
    "\n",
    "plt.figure(figsize=(8, 4))\n",
    "plt.subplot(121)\n",
    "plt.title(\"2D PCA\")\n",
    "plt.scatter(Z[:, 0], Z[:, 1])\n",
    "plt.subplot(122)\n",
    "plot_diagrams(PDs)\n",
    "plt.title(\"Persistence Diagram\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Unfortunately, the sample rate is just to high and the signal is just too messy for this algorithm to work.  We will have to do some more sophisticated preprocessing before applying the algorithm"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h1>Audio Novelty Functions And Music vs Speech</h1>\n",
    "One way to deal with the fact that music is both messy and at a high sampling rate is to derive something called the \"audio novelty function,\" which is designed explicitly to pick up on rhythmic events.  To see how it's motivated, let's look at the <a href = \"https://en.wikipedia.org/wiki/Spectrogram\">audio spectrogram</a> of this song\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from MusicFeatures import *\n",
    "\n",
    "#Compute the power spectrogram and audio novelty function\n",
    "winSize = 512\n",
    "hopSize = 256\n",
    "plt.figure()\n",
    "(S, novFn) = getAudioNoveltyFn(X, Fs, winSize, hopSize)\n",
    "plt.imshow(np.log(S.T), cmap = 'afmhot', aspect = 'auto')\n",
    "plt.title('Log-frequency power spectrogram')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You might notice that there are vertical streaks in a semi-periodic pattern.  These correspond to \"broadband percussive events,\" or, on other words, likely onsets for beats when drums occur.  An audio novelty function is derived from a spectrogram by looking at the difference between successive frames to try to pick up on this.  The code below extracts the audio novelty function and displays it for the same audio snippet."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(8, 4))\n",
    "#Plot the spectrogram again\n",
    "plt.subplot(211)\n",
    "plt.imshow(np.log(S.T), cmap = 'afmhot', aspect = 'auto')\n",
    "plt.ylabel('Frequency Bin')\n",
    "plt.title('Log-frequency power spectrogram')\n",
    "\n",
    "#Plot the audio novelty function\n",
    "plt.subplot(212)\n",
    "plt.plot(np.arange(len(novFn))*hopSize/float(Fs), novFn)\n",
    "plt.xlabel(\"Time (Seconds)\")\n",
    "plt.ylabel('Audio Novelty')\n",
    "plt.xlim([0, len(novFn)*float(hopSize)/Fs])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Not only is the audio novelty function a cleaner signal, but it is also at a much lower sample rate.  Since the \"hop size\" between each spectrogram window is 256 samples, the temporal resolution is coarser by that factor.  \n",
    "<h2>Music Audio Novelty Embedding</h2>\n",
    "\n",
    "Let's now try our sliding window with a snippet of the audio novelty function of the previous example instead of the raw audio:<BR>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "(S, novFn) = getAudioNoveltyFn(X, Fs, winSize, hopSize)\n",
    "\n",
    "#Take the first 3 seconds of the novelty function\n",
    "fac = int(Fs/hopSize)\n",
    "novFn = novFn[fac*4:fac*7]\n",
    "\n",
    "#Make sure the window size is half of a second, noting that\n",
    "#the audio novelty function has been downsampled by a \"hopSize\" factor\n",
    "dim = 20\n",
    "Tau = (Fs/2)/(float(hopSize)*dim)\n",
    "dT = 1\n",
    "Y = getSlidingWindowInteger(novFn, dim, Tau, dT)\n",
    "print(\"Y.shape = \", Y.shape)\n",
    "#Mean-center and normalize\n",
    "Y = Y - np.mean(Y, 1)[:, None]\n",
    "Y = Y/np.sqrt(np.sum(Y**2, 1))[:, None]\n",
    "\n",
    "PDs = ripser(Y, maxdim=1)['dgms']\n",
    "pca = PCA()\n",
    "Z = pca.fit_transform(Y)\n",
    "\n",
    "plt.figure(figsize=(12, 6))\n",
    "plt.subplot(121)\n",
    "plt.title(\"2D PCA\")\n",
    "plt.scatter(Z[:, 0], Z[:, 1])\n",
    "plt.subplot(122)\n",
    "plot_diagrams(PDs)\n",
    "plt.title(\"Persistence Diagram\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "<h2>Speech Example</h2>\n",
    "In our final experiment in this module, we will look at a sliding window embedding of the audio novelty function on a speech excerpt which does not have a clear rhythmic structure (courtesy of <a href = \"http://marsyasweb.appspot.com/download/data_sets/\">http://marsyasweb.appspot.com/download/data_sets/</a>).  Click on the cell below to load the speech audio, and click on the cell below that to run the sliding window embedding + persistent homology"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Read in the audio file.  Fs is the sample rate, and\n",
    "#X is the audio signal\n",
    "Fs, X = scipy.io.wavfile.read(\"speech.wav\")\n",
    "X = X/(2.0**15)\n",
    "(S, novFn) = getAudioNoveltyFn(X, Fs, winSize, hopSize)\n",
    "plt.figure()\n",
    "plt.plot(np.arange(len(novFn))*hopSize/float(Fs), novFn)\n",
    "plt.xlabel(\"Time (Seconds)\")\n",
    "plt.title(\"Audio Novelty Function for Speech\")\n",
    "Audio('speech.wav')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Get the novelty function for the first three seconds, and use the \n",
    "#exact same parameters as before\n",
    "novFn = novFn[0:int((Fs/hopSize)*3)]\n",
    "dim = 20\n",
    "Tau = (Fs/2)/(float(hopSize)*dim)\n",
    "dT = 1\n",
    "Y = getSlidingWindowInteger(novFn, dim, Tau, dT)\n",
    "print(\"Y.shape = \", Y.shape)\n",
    "#Mean-center and normalize\n",
    "Y = Y - np.mean(Y, 1)[:, None]\n",
    "Y = Y/np.sqrt(np.sum(Y**2, 1))[:, None]\n",
    "\n",
    "PDs = ripser(Y, maxdim=1)['dgms']\n",
    "pca = PCA()\n",
    "Z = pca.fit_transform(Y)\n",
    "\n",
    "plt.figure(figsize=(12, 6))\n",
    "plt.subplot(121)\n",
    "plt.title(\"2D PCA\")\n",
    "plt.scatter(Z[:, 0], Z[:, 1])\n",
    "plt.subplot(122)\n",
    "plot_diagrams(PDs[1], labels=['H1'])\n",
    "plt.title(\"Persistence Diagram\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h1>Summary</h1>\n",
    "<ul>\n",
    "<li>Biphonation is present in nature, for example, in horse whinnies.  In this module we used audio, sampled at 44100hz, to record said whinnies.</li>\n",
    "<li>Due to noise and other signal artifacts, it is sometimes necessary to search around to find regions of biphonation that are clean enough to be observed via persistent homology</li>\n",
    "<li>Summary features are often better than raw data in practical scenarios.  In music examples, due to the high sampling rate of audio and the messy signal, it is often advantageous to compute summary signals at a lower temporal resolution when analyzing rhythm periodicities.</li>\n",
    "<li>After the proper preprocessing, TDA on sliding window embeddings can pick up on the presence of rhythmic periodicities in music.</li>\n",
    "</ul>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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
}