examples/basic-audio-features.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "DSC160 Data Science and the Arts - Twomey - Spring 2020 - [dsc160.roberttwomey.com](http://dsc160.roberttwomey.com)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Basic Audio Features\n",
    "\n",
    "This notebook demonstrates how to read, write, and calculate simple time-domain statistics on digital audio. It depends on the [numpy](https://numpy.org/), [matplotlib](https://matplotlib.org/), [seaborn](https://seaborn.pydata.org/), and [LibROSA](https://librosa.github.io/librosa/) libraries. The contents are adapted for DSC160 from the Stanford Music Information Retrieval workshop tutorials on [musicinformationretrieval.com](musicinformationretrieval.com).\n",
    "\n",
    "__Sections:__ \n",
    "- Digital Audio Basics\n",
    "  - [Setup](#Setup)\n",
    "  - [Reading Audio Files](#Reading-Audio-Files)\n",
    "  - [Displaying Waveforms](#Displaying-Waveforms)\n",
    "- Temporal Features\n",
    "  - [Zero-Crossings](#Zero-Crossings)\n",
    "  - [Energy](#Energy)\n",
    "  - [Autocorrelation and Pitch Estimation](#Autocorrelation-and-Pitch-Estimation)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setup"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Tell matplotlib to render plots within the jupter frontend."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you have not previously installed LibROSA, uncomment the following line and run to install."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#!pip install librosa --user"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Basic imports"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# visualization\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# sound processing\n",
    "import librosa\n",
    "import librosa.display\n",
    "\n",
    "# to play audio inline in ipython/jupyter notebooks\n",
    "from IPython.display import Audio\n",
    "\n",
    "# numpy\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Reading Audio Files"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Read in a digital audio file, loading it as a waveform (time series of samples) `y` with sample rate `sr`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "filename = librosa.util.example_audio_file()\n",
    "y, sr = librosa.load(filename)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(type(y), type(sr))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How many samples do we have, and what is the sample rate?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(y.shape, sr)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Given the sample rate and array of samples, how long is our audio file?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "length = y.shape[0]/sr\n",
    "print(\"file is {:.4} seconds long\".format(length))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "librosa also has a method for calculating this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "librosa.get_duration(y, sr)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What data type are our samples stored in?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(y.dtype)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What is the range of values in our samples?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(min(y), max(y))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "NOTE: Typically for a floating point sample, amplitudes (values) will be between normalized between -1.0 and 1.0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Play an audio file in notebook"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It's useful to be able to hear our audio directly in notebook. Here we use the IPython Audio class to create a playable widget inline:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Audio(data=y, rate=sr)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note: `Audio` can also accept a numpy array. Let's synthesize a pure tone at 440 Hz:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sr = 22050 # sample rate\n",
    "T = 2.0 # duration in seconds\n",
    "t = np.linspace(0, T, int(T*sr), endpoint=False) # generate the time stamps for samples\n",
    "x = 0.5*np.sin(2*np.pi*440*t) # generate the amplitudes (pure sine wave at 440 Hz)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Audio(data=x, rate=sr)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And we can write it to disk (as a [WAV file](https://en.wikipedia.org/wiki/WAV)) using `librosa.output.write_wav`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "librosa.output.write_wav('tone_440.wav', x, sr)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Displaying Waveforms"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the audio file as a waveform. The horizontal axis (X) corresponds to time (minutes:seconds), and vertical axis (Y) corresponds to the instantaneous sound amplitude (pressure)/magnitude of the digital audio signal. \n",
    "\n",
    "NOTE: we are using `matplotlib.pyplt` commands and `librosa.display` interchangeably. `librosa.display` builds on matplotlib."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(12,6))\n",
    "librosa.display.waveplot(y, sr=sr)\n",
    "plt.title('Waveform')\n",
    "plt.tight_layout()\n",
    "plt.ylabel(\"amplitude\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we zoom in on a subrange of the time series (say items 2000-2300 in the array), we can see individual samples and their values:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n0 = 2000\n",
    "n1 = 2300\n",
    "plt.figure(figsize=(12, 6))\n",
    "plt.scatter(range(n0, n1), y[n0:n1])\n",
    "plt.plot(range(n0, n1), y[n0:n1])\n",
    "plt.tight_layout()\n",
    "plt.xlabel(\"sample #\")\n",
    "plt.ylabel(\"amplitude\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can convert those sample indices (`n0` and `n1` above) to times (seconds) based on the sample rate for this file:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "t0, t1 = librosa.core.samples_to_time((n0, n1), sr)\n",
    "print(\"samples from {:0.4f}s - {:0.4}s in the input audio file\".format(t0, t1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Displaying Spectogram\n",
    "\n",
    "A spectrogram is a visual representation of the spectrum of frequency components of a signal as it changes over time. [wikipedia](https://en.wikipedia.org/wiki/Spectrogram). The plot below uses the [Short-Time Fourier Transform](https://en.wikipedia.org/wiki/Short-time_Fourier_transform) to calculate the magnitudes of frequency components. We'll consider this in following lecture."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Y = librosa.stft(y)\n",
    "Ydb = librosa.amplitude_to_db(abs(Y))\n",
    "plt.figure(figsize=(14, 5))\n",
    "librosa.display.specshow(Ydb, sr=sr, x_axis='time', y_axis='hz')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "EXTENSION: repeat the above plots (waveform and spectogram) with the sine tone you generated at the beginning of the notebook."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Temporal Features of Audio Files\n",
    "\n",
    "The following sections demonstrate the calculation of basic temporal domain features on our input files."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Zero-Crossings\n",
    "\n",
    "The [zero-crossing rate](https://en.wikipedia.org/wiki/Zero-crossing_rate) indicates the number of times that an audio signal crosses the horizontal axis (switches from positive to negative, or negative to positive). This feature is used in speech recognition and music information retrieval, and is a key feature of percussive sounds (drums). For a simple monophonic, tonal signals, the zcr can serve as a rudimentary pitch estimation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First we will load an audio file:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x, sr = librosa.load('audio/simple_loop.wav')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can listen to it:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Audio(x, rate=sr)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What does the waveform look like?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(14, 5))\n",
    "librosa.display.waveplot(x, sr=sr)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we zoom in on a small subset of samples, we see a few distinct zero crossings:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n0 = 6500\n",
    "n1 = 7500\n",
    "plt.figure(figsize=(14, 5))\n",
    "plt.plot(x[n0:n1])\n",
    "plt.plot([0.0]*(n1-n0))\n",
    "plt.ylabel(\"amplitude\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "By visual inspection, there appear to be 5 zero crossings. What does librosa compute? "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "zero_crossings = librosa.zero_crossings(x[n0:n1], pad=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Interestingly, the result has as many values as our input sample"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "zero_crossings.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "That computed a binary mask where `True` indicates the presence of a zero crossing. To find the total number of zero crossings, use `sum`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(sum(zero_crossings))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To find the zero-crossing rate over time, we can use librosa's `zero_crossing_rate` method:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "zcrs = librosa.feature.zero_crossing_rate(x)\n",
    "print(zcrs.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(14, 5))\n",
    "plt.plot(zcrs[0])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The reason for the high rate near the beginning is because the silence oscillates quietly around zero:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(14, 5))\n",
    "plt.plot(x[:1000])\n",
    "plt.ylim(-0.0001, 0.0001)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A simple hack around this is to add a small constant before computing the zero crossing rate:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "zcrs = librosa.feature.zero_crossing_rate(x + 0.0001)\n",
    "plt.figure(figsize=(14, 5))\n",
    "plt.plot(zcrs[0])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "EXTENSION: try this calculation of zero-crossing rate with an audio file of your choosing. Does the zero-crossing rate capture valuable information for other kinds of sounds? What about an excerpt from a complex song?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Energy\n",
    "\n",
    "The [energy](https://en.wikipedia.org/wiki/Energy_(signal_processing)) of a signal corresponds to the total magntiude of the signal. For audio signals, that roughly corresponds to how loud the signal is. The energy in a signal is defined as\n",
    "\n",
    "$$ \\sum_n \\left| x(n) \\right|^2 $$\n",
    "\n",
    "The **root-mean-square (RMS) energy** in a signal is defined as\n",
    "\n",
    "$$ \\sqrt{ \\frac{1}{N} \\sum_n \\left| x(n) \\right|^2 } $$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Working with the `simple_loop.wav` from above, lets estimate the amount of energy in the signal at various points in time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x, sr = librosa.load('audio/simple_loop.wav')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(14, 5))\n",
    "librosa.display.waveplot(x, sr=sr)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First we will compute the short-time energy using list comprehension. \n",
    "\n",
    "\"Hop length\" below tells us how far over to \"hop\" for each successive frame of analysis (step size). \"Frame length\" tells us how wide the window is over which we wish to calculate the energy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hop_length = 256\n",
    "frame_length = 512"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "energy = np.array([\n",
    "    sum(abs(x[i:i+frame_length])**2)\n",
    "    for i in range(0, len(x), hop_length)\n",
    "])\n",
    "\n",
    "energy.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Secondly, we compute the root-mean-square (RMS) energy using [`librosa.feature.rms`](https://librosa.github.io/librosa/generated/librosa.feature.rms.html):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rmse = librosa.feature.rms(x, frame_length=frame_length, hop_length=hop_length, center=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rmse.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rmse = rmse[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot both the energy and RMSE along with the waveform:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "frames = range(len(energy))\n",
    "t = librosa.frames_to_time(frames, sr=sr, hop_length=hop_length)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(14, 5))\n",
    "librosa.display.waveplot(x, sr=sr, alpha=0.4)\n",
    "plt.plot(t, energy/energy.max(), 'r--') # normalized for visualization\n",
    "plt.plot(t[:len(rmse)], rmse/rmse.max(), color='g') # normalized for visualization\n",
    "plt.legend(('Energy', 'RMSE'))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "NOTE: The Energy plot and RMSE are slightly out of sync, temporally. This is due to the use of `center=True` in `librosa.feature.rms`, which centers the RMSE values calculated over the feature window in the output."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Autocorrelation and Pitch Estimation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The [autocorrelation](http://en.wikipedia.org/wiki/Autocorrelation) of a signal describes the similarity of a signal against a time-shifted version of itself. For a signal $x$, the autocorrelation $r$ is:\n",
    "\n",
    "$$ r(k) = \\sum_n x(n) x(n-k) $$\n",
    "\n",
    "In this equation, $k$ is often called the **lag** parameter. $r(k)$ is maximized at $k = 0$ and is symmetric about $k$.\n",
    "\n",
    "The autocorrelation is useful for finding repeated patterns in a signal. For example, at short lags, the autocorrelation can tell us something about the signal's fundamental frequency. For longer lags, the autocorrelation may tell us something about the tempo of a musical signal. \n",
    "\n",
    "This section will walk you through estimating pitch using autocorrelation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "First let's load a file:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "x, sr = librosa.load('audio/c_strum.wav')\n",
    "Audio(x, rate=sr)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "...and view its waveform"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "plt.figure(figsize=(14, 5))\n",
    "librosa.display.waveplot(x, sr)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are two ways we can compute the autocorrelation in Python. The first method is [`numpy.correlate`](http://docs.scipy.org/doc/numpy/reference/generated/numpy.correlate.html). Because the autocorrelation produces a symmetric signal, we only care about the \"right half\" (so we grab from `len(x)-1` to the end)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "r = np.correlate(x, x, mode='full')[len(x)-1:]\n",
    "print(x.shape, r.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the autocorrelation as calculated with numpy:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "plt.figure(figsize=(14, 5))\n",
    "plt.plot(r[:10000])\n",
    "plt.xlabel('Lag (samples)')\n",
    "plt.xlim(0, 10000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "notes"
    }
   },
   "source": [
    "The second method uses [`librosa.autocorrelate`](https://librosa.github.io/librosa/generated/librosa.core.autocorrelate.html):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "r = librosa.autocorrelate(x, max_size=10000)\n",
    "print(r.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "plt.figure(figsize=(14, 5))\n",
    "plt.plot(r)\n",
    "plt.xlabel('Lag (samples)')\n",
    "plt.xlim(0, 10000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These look pretty similar. The second method, `librosa.autocorrelate` conveniently only keeps one half of the autocorrelation function, since the autocorrelation is symmetric. Also, the `max_size` parameter prevents unnecessary calculations."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Pitch Estimation with Autocorrelation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The autocorrelation is used to find repeated patterns within a signal. For musical signals, a repeated pattern can correspond to a pitch period. We can therefore use the autocorrelation function to estimate the pitch in a musical signal.\n",
    "\n",
    "Here is a recording of an oboe playing a single note (C6)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x, sr = librosa.load('audio/oboe_c6.wav')\n",
    "Audio(x, rate=sr)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compute and plot the autocorrelation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "r = librosa.autocorrelate(x, max_size=5000)\n",
    "plt.figure(figsize=(14, 5))\n",
    "plt.plot(r[:200])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The autocorrelation always has a maximum at zero, i.e. at zero lag. We want to identify the maximum outside of the peak centered at zero. Therefore, we might choose only to search within a range of reasonable pitches:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# pick low and high bounds based on a reasonable range of notes\n",
    "# 120 is past the highest key on a piano, 12 is below the lowest\n",
    "midi_hi = 120.0 \n",
    "midi_lo = 12.0\n",
    "\n",
    "# convert midi pitches to frequency (hz)\n",
    "f_hi = librosa.midi_to_hz(midi_hi)\n",
    "f_lo = librosa.midi_to_hz(midi_lo)\n",
    "\n",
    "# calculate the corresponding period (seconds) of each bounding frequency\n",
    "t_lo = sr/f_hi\n",
    "t_hi = sr/f_lo"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(f_lo, f_hi)\n",
    "print(t_lo, t_hi)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Set invalid pitch candidates to zero (zero all values in the autocorrelation result below `t_lo` or above `t_hi`):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "r[:int(t_lo)] = 0\n",
    "r[int(t_hi):] = 0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "plot what remains"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(14, 5))\n",
    "plt.plot(r[:1400])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Find the location of the maximum:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "t_max = r.argmax()\n",
    "print(t_max)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, estimate the pitch in Hertz:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "float(sr)/t_max"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Indeed, that is very close to the true frequency of C6 (the note that was played by the oboe):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "librosa.midi_to_hz(84)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "EXTENSION: Repeat the above pitch estimation through autocorrelation using some other tonal (pure tone) audio recording of your choice, for instance a recording of someone singing a certain note. Can you correctly estimate the pitch?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## References\n",
    "- International Society for Music Information Retrieval (ISMIR) [https://ismir.net/](https://ismir.net/)\n",
    "- Laboratory for the Recognition and Organization of Speech and Audio at Columbia University: [LabROSA](https://labrosa.ee.columbia.edu/)\n",
    "  - LibROSA [https://librosa.github.io/librosa/](https://librosa.github.io/librosa/)\n",
    "- Brian McFee - SciPy 2015 Talk on Audio Processing and Music Information Retrieval (MIR) with LibROSA: https://www.youtube.com/watch?v=MhOdbtPhbLU\n",
    "  - [website](https://bmcfee.github.io/) [paper](https://bmcfee.github.io/papers/scipy2015_librosa.pdf)\n",
    "- Code examples above are adapted from [musicinformationretrieval.com](https://musicinformationretrieval.com/), including:\n",
    "  - [Zero-Crossing Rate](https://musicinformationretrieval.com/zcr.html)\n",
    "  - [Energy](https://musicinformationretrieval.com/energy.html)\n",
    "  - [Autocorrelation](https://musicinformationretrieval.com/autocorrelation.html)\n"
   ]
  },
  {
   "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.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}