03 - Vocal activity detection.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Part 3 - Vocal activity\n",
    "\n",
    "In this notebook, you'll implement a basic vocal activity detector.  This will build on the exercise in part 2 where you built a silence detector.\n",
    "\n",
    "## Setup\n",
    "\n",
    "Again, we begin by setting up our envrionment.  This time, we'll do it all in one block."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "import librosa\n",
    "import librosa.display\n",
    "from IPython.display import Audio"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# We'll work with an example provided in the repository.\n",
    "# There are many files in this directory, but we'll start with the first one.\n",
    "# The later examples (higher numbers) have more difficult noise conditions.\n",
    "\n",
    "y, sr = librosa.load('detection_examples/voice_detection_1.wav')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Play it back, so we know what our data sounds like\n",
    "Audio(data=y, rate=sr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Let's visualize the spectrogram to see what the data looks like\n",
    "D = librosa.stft(y)\n",
    "S = np.abs(D)\n",
    "S_db = librosa.amplitude_to_db(S, ref=np.max)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "librosa.display.specshow(S_db, x_axis='time', y_axis='log')\n",
    "plt.colorbar()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# What can we see here?\n",
    "\n",
    "- Unlike before, there's energy everywhere!\n",
    "- The background noise is \"non-stationary\", meaning that its spectrum evolves over time.  Because of this, a simple high-pass filter isn't going to work.\n",
    "- What else can we do?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### PCEN\n",
    "\n",
    "As an alternative to decibel scaling, we can use the \"per-channel energy normalization\" (PCEN) technique.  This was developed by researchers at Google for spotting keywords (\"okay google!\") in their \"Home\" product.  It turns out to work well for many problems, but it's particularly good at suppressing complex background noise.\n",
    "\n",
    "The key idea is to compute a spectrogram, and then compare each time-frequency measurement to the average energy in that frequency over the preceding few time steps.  This will amplify unpredictable changes, which often correspond to speech.\n",
    "\n",
    "We can compute this by `librosa.pcen`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "S_pcen = librosa.pcen(S)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# What does this look like?\n",
    "librosa.display.specshow(S_pcen, x_axis='time', y_axis='log')\n",
    "plt.colorbar()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Interpreting PCEN\n",
    "\n",
    "PCEN preserves the shape of the spectrogram, but the range of values is quite different.  Unlike the decibel scale, PCEN is always non-negative.\n",
    "\n",
    "After applying PCEN, we can look at the total normalized energy averaged in each time step."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# This computes the mean along the 0 (first, vertical) axis\n",
    "pcen_avg = np.mean(S_pcen, axis=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# We can then plot the curve just like we did energy in the previous examples\n",
    "times = librosa.times_like(S_pcen)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Now we can plot the data\n",
    "plt.plot(times, pcen_avg)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# What's going on here?\n",
    "\n",
    "Remember that PCEN compares each frequency measurement to its previous time steps to determine \"novelty\".  At the beginning of the signal, there's no history yet, so everything seems novel, hence the large value at `t=0`.\n",
    "\n",
    "After the filter warms up, you can look at spikes in the `pcen_avg` curve to find points of interest.\n",
    "\n",
    "We can find those points by using the \"onset detection\" method `librosa.onset.onset_detect`.  This will try to find peaks in a given curve that are well separated from each other in time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# we can get event onsets in units of \"time\"\n",
    "librosa.onset.onset_detect(onset_envelope=pcen_avg, units='time')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# or in units of \"samples\". This is useful if you want to slice `y` afterward\n",
    "onsets = librosa.onset.onset_detect(onset_envelope=pcen_avg, units='samples')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# We'll slice out from the first onset plus one second\n",
    "y_slice = y[onsets[2]:onsets[2]+sr]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Audio(data=y_slice, rate=sr)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Hmmm....\n",
    "\n",
    "Maybe some filtering can still help here!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Noise!\n",
    "\n",
    "These signals have a lot of energy in the low frequencies.\n",
    "\n",
    "Let's use scipy to filter out the high frequencies."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import scipy.signal"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# We'll make a band-pass, order-6 butterworth filter\n",
    "# Human voice lives principally in the range 300-3000 Hz\n",
    "# We can filter down to just that range, and reduce noise from other sources\n",
    "\n",
    "b, a = scipy.signal.butter(6, [300, 3000], btype='bandpass', fs=sr)\n",
    "\n",
    "# Now that we have our coefficients, we can filter the signal\n",
    "y_low = scipy.signal.filtfilt(b, a, y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Let's listen to our output\n",
    "Audio(data=y_low, rate=sr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Now let's compute the spectrogram again and plot it\n",
    "D_low = librosa.stft(y_low)\n",
    "S_low = np.abs(D_low)\n",
    "S_low_pcen = librosa.pcen(S_low)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "librosa.display.specshow(S_low_pcen, x_axis='time', y_axis='log')\n",
    "plt.colorbar()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# That looks a little better...\n",
    "\n",
    "\n",
    "But pure PCEN is not going to solve the entire problem for us.  \n",
    "Maybe we can do something better than just detecting onsets?\n",
    "\n",
    "Let's try to group frames by similar timbre.  We'll do this in two steps.\n",
    "\n",
    "First, we'll extract the \"Mel Frequency Cepstral Coefficients\" (MFCCs) from each frame of audio.  These represent the coarse shape of a spectrum without being too sensitive to pitch, and are a good first step for trying to capture timbre.\n",
    "\n",
    "Second, after we have the MFCC's, we'll put them into a clustering model to group them into two \"clusters\", which should correspond to \"voice\" and \"not voice\"."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# MFCC's are calculated by the mfcc function in librosa.feature:\n",
    "\n",
    "mfcc = librosa.feature.mfcc(y=y, sr=sr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# While we're at it, let's apply it to the filtered signal as well\n",
    "\n",
    "mfcc_low = librosa.feature.mfcc(y=y_low, sr=sr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# What do they look like?\n",
    "librosa.display.specshow(mfcc, x_axis='time')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "librosa.display.specshow(mfcc_low, x_axis='time')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Hard to see with human eyes\n",
    "\n",
    "But machines might do better!\n",
    "\n",
    "Let's now try to cluster the frames."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# For this, we'll need the \"cluster\" module from scikit-learn\n",
    "import sklearn.cluster"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# We'll build a K-means clusterer\n",
    "# With n_clusters=2, this groups each frame into one of two \"clusters\"\n",
    "# Hopefully, one cluster will get voice and one will get everything else!\n",
    "Kmeans = sklearn.cluster.KMeans(n_clusters=2)\n",
    "\n",
    "# Building the clusterer doesn't do anything yet.  We need to fit it to data.\n",
    "# scikit-learn expects each data sample to be a row of the matrix, so we have\n",
    "# to transpose our data by saying \".T\"\n",
    "\n",
    "clusters = Kmeans.fit_predict(mfcc.T)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# This cell will make two plots with a common time axis\n",
    "# This is handy for visualizing aligned features\n",
    "\n",
    "plt.figure(figsize=(12, 5))\n",
    "axes = plt.subplot(2,1,1)\n",
    "librosa.display.specshow(S_low_pcen, x_axis='time')\n",
    "plt.subplot(2,1,2, sharex=axes)\n",
    "plt.plot(times, clusters)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# There's some information here, but it could be better.\n",
    "# Let's look at the band-filtered signal instead\n",
    "\n",
    "Kmeans_band = sklearn.cluster.KMeans(n_clusters=2)\n",
    "\n",
    "clusters = Kmeans.fit_predict(mfcc_low.T)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# This cell will make two plots with a common time axis\n",
    "# This is handy for visualizing aligned features\n",
    "\n",
    "plt.figure(figsize=(12, 5))\n",
    "axes = plt.subplot(2,1,1)\n",
    "librosa.display.specshow(S_low_pcen, x_axis='time')\n",
    "plt.subplot(2,1,2, sharex=axes)\n",
    "plt.plot(times, clusters)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Much more stable.  We can do even better by adding some temporal context\n",
    "# The above was classifying each frame independently, but we can stack neighboring\n",
    "# frames together to get a bit more information\n",
    "\n",
    "mfcc_stack = librosa.feature.stack_memory(mfcc_low)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Kmeans_band = sklearn.cluster.KMeans(n_clusters=2)\n",
    "\n",
    "clusters = Kmeans.fit_predict(mfcc_stack.T)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# This cell will make two plots with a common time axis\n",
    "# This is handy for visualizing aligned features\n",
    "\n",
    "plt.figure(figsize=(12, 5))\n",
    "axes = plt.subplot(2,1,1)\n",
    "librosa.display.specshow(S_low_pcen, x_axis='time')\n",
    "plt.subplot(2,1,2, sharex=axes)\n",
    "plt.plot(times, clusters)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Okay, what now?\n",
    "\n",
    "We have a detector that works pretty well, and there are many ways to improve this process.\n",
    "\n",
    "But how can we use the cluster assignments?\n",
    "\n",
    "### First question:\n",
    "\n",
    "When does the signal transition from one cluster to the other?  These are the \"change points\".\n",
    "\n",
    "We can find these with a few `numpy` operations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Where are the frames such that\n",
    "#   clusters[t] != clusters[t+1] ?\n",
    "#\n",
    "# We can only check this up to the second-to-last frame\n",
    "# In slice notation,\n",
    "#     [:-1] means \"all but the last\"\n",
    "#     [1:] means \"all but the first\"\n",
    "#\n",
    "# np.where will give us a list of indices that satisfy the condition we're interested in\n",
    "# note: it will actually produce a list of *lists* of indices, one for each dimension\n",
    "#       since we only have one dimension here (time), we can just take the first element\n",
    "\n",
    "np.where(clusters[:-1] != clusters[1:])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "changes = np.where(clusters[:-1] != clusters[1:])[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# This array will contain the frame indices of each change point\n",
    "changes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# We can convert them to sample indices:\n",
    "change_samples = librosa.frames_to_samples(changes)\n",
    "\n",
    "# Or to time indices in seconds:\n",
    "change_time = librosa.frames_to_time(changes)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "change_time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Once we have the change-points in time, we can plot markers on top of the spectrogram\n",
    "plt.figure(figsize=(12, 5))\n",
    "librosa.display.specshow(S_low_pcen, x_axis='time')\n",
    "\n",
    "# Make vertical lines, positioned at the change times,\n",
    "# starting at 0 Hz and going all the way up to sr/2\n",
    "# and color them white\n",
    "plt.vlines(change_time, 0, sr/2, color='w')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Summary\n",
    "\n",
    "What have we done in this notebook?\n",
    "\n",
    "1. Worked with complex acoustic scenes\n",
    "2. Used PCEN to more clearly visualize the contents of the sound\n",
    "3. Used a band-pass filter to focus on human vocal range\n",
    "4. Extracted timbral features\n",
    "5. Clustered the features to separate \"voice\" from \"not voice\"\n",
    "6. Detected change-points in the signal, and inferred time markers\n",
    "\n",
    "What could we do next? Some ideas:\n",
    "\n",
    "- Source separation: can we remove noise that interferes with voice?\n",
    "- Temporally constrained detection: can we make more locally contiguous detections?\n",
    "- Modeling larger collections: this might work for one speaker, but what about more generally?"
   ]
  },
  {
   "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.6.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}