getting_started_student.ipynb

{
  "cells": [
    {
      "cell_type": "markdown",
      "source": [
        "# Hands-On 1: Audio Feature Extraction\r\n",
        "\r\n",
        "For this first hands-on session, we are going to investigate the extraction of audio features. This is the first step for designing a classification model (in future hands-on sessions). <br>\r\n",
        "As the recording of audio signals using a microphone will also be covered in a future hands-on session, we will start here by using sounds available in the ``ESC-50`` dataset.\r\n",
        "> Karol J. Piczak, 2015, \"ESC: Dataset for Environmental Sound Classification\"  \r\n",
        "> https://doi.org/10.7910/DVN/YDEPUT, Harvard Dataverse, V2 ([Available here](https://dataverse.harvard.edu/dataset.xhtml?persistentId=doi:10.7910/DVN/YDEPUT))\r\n",
        "\r\n",
        "Useful functions to select, read and play the dataset sounds are provided in the ``utils`` folder. The following required packages should be installed:\r\n",
        "- soundfile\r\n",
        "- scipy\r\n",
        "- sounddevice\r\n",
        "- librosa\r\n",
        "- seaborn\r\n",
        "\r\n",
        "You can install them with ``!pip install <packagename>`` (decomment the lines in the following cell). In case of an error with the utils folder (folder not found), you may need to launch Jupyter with the directory where the code to execute is located. To do so, open the Anaconda Prompt (if you are using Anaconda) and type ``jupyter notebook --notebook-dir=$YOUR PATH$``. <br>\r\n",
        "\r\n",
        "To ensure you are catching the content of this notebook, we leave you with an infinitesimal amount of **code to write**. \r\n",
        "\r\n",
        "You will find the zones to be briefly filled  with a ``### TO COMPLETE`` in the cells below."
      ],
      "metadata": {}
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "source": [
        "'''\r\n",
        "!pip install -U soundfile\r\n",
        "!pip install -U sounddevice\r\n",
        "!pip install -U pygame\r\n",
        "!pip install -U librosa\r\n",
        "!pip install -U seaborn\r\n",
        "'''\r\n",
        "import numpy as np\r\n",
        "import matplotlib.pyplot as plt\r\n",
        "import soundfile as sf\r\n",
        "from scipy import signal\r\n",
        "import sounddevice as sd\r\n",
        "import librosa # For audio signal computations as MFCC\r\n",
        "from scipy.fftpack import dct\r\n",
        "\r\n",
        "\"Self created functions\"\r\n",
        "from utils_ import getclass, getname, gen_allpath, plot_audio, plot_specgram"
      ],
      "outputs": [],
      "metadata": {}
    },
    {
      "cell_type": "markdown",
      "source": [
        "We then create a matrix with path names of height H=50 classes and W=40 sounds per class. This will give you simple access to any sound from the dataset."
      ],
      "metadata": {}
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "source": [
        "classnames, allpath = gen_allpath() # Note: this function contains an implicit input \"folder\", where you can change the path to ESC-50 dataset.\r\n",
        "print('The classes are : \\n')\r\n",
        "for ind, elem in enumerate(classnames) :\r\n",
        "    print(ind, elem)"
      ],
      "outputs": [],
      "metadata": {}
    },
    {
      "cell_type": "markdown",
      "source": [
        "You can now select a sound from a given class using ``allpath[class_index,sound_index]``. For example, the first sound of the ``Cow`` class is accessed with ``allpath[3,0]`` (note that this index may change from one operating system to the other) and the following cell plays the sound:"
      ],
      "metadata": {}
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "source": [
        "sound = allpath[3,29]\r\n",
        "x, fs = sf.read(sound)\r\n",
        "print(f\"Playing a \\\"{getclass(sound)}\\\"\")\r\n",
        "sd.play(x, fs)"
      ],
      "outputs": [],
      "metadata": {}
    },
    {
      "cell_type": "markdown",
      "source": [
        "We now ask you to complete the cells below."
      ],
      "metadata": {}
    },
    {
      "cell_type": "markdown",
      "source": [
        "## 1) Resampling and filtering\r\n",
        "\r\n",
        "Most probably your circuit board will sample the analog audio signal at a frequency $f_s = 11025$ Hz. <br>\r\n",
        "However, the audio provided in the ``ESC-50`` dataset are sampled with $f_s = 44100$ Hz, you should thus downsample each audio signal to keep coherency with your real setup. There are 2 solutions:\r\n",
        "- Rewrite a new dataset with the downsampled audio signals.\r\n",
        "- Downsample each audio which is read.\r\n",
        "\r\n",
        "We provide you with the second one.\r\n",
        "\r\n",
        "***\r\n",
        "#### <u> The following derivations are not necessary for the rest of this notebook, but are still provided for the curious students... </u>\r\n",
        "\r\n",
        "Let us consider one original audio signal from the dataset and denote it $x[n]$, for $n=0,\\dots,N-1$.\r\n",
        "\r\n",
        "The downsampled signal $y$ can be written as \r\n",
        "\r\n",
        "$$\r\n",
        "    y[m] = w[mM],\\quad \\text{with}\\ w[k] = (h \\ast x)[k] = \\sum_{n=-\\infty}^{\\infty} h[n]x[k-n],\r\n",
        "$$\r\n",
        "\r\n",
        "where $h$ is a discrete low-pass filter and $M$ is the downsampling factor, here $M=4$. <br>\r\n",
        "\r\n",
        "We can expand both $y$ and $w$ according to their Fourier series (DTFT) $Y$ and $W$, respectively, so that:\r\n",
        "\r\n",
        "$$\r\n",
        "    y[m] = \\frac{1}{2\\pi} \\int_0^{2\\pi} Y(e^{j\\Omega}) e^{jm\\Omega} d\\Omega \\tag{1}\r\n",
        "$$\r\n",
        "$$\r\n",
        "    w[mM] = \\frac{1}{2\\pi} \\int_{0}^{2\\pi} W(e^{j\\Omega}) e^{jmM\\Omega} d\\Omega = \\frac{1}{2\\pi} \\sum_{k=0}^{M-1} \\int_{2\\pi k/M}^{2\\pi(k+1)/M} W(e^{j\\Omega}) e^{jmM\\Omega} d\\Omega \\tag{2}\r\n",
        "    $$\r\n",
        "\r\n",
        "Regarding $w$, applying the change of variable $\\Omega \\leftarrow \\Omega-2\\pi k/M$ to each integral of the sum, and changing $k \\leftarrow M - k$, we can further write\r\n",
        "\r\n",
        "$$\r\n",
        "   \\textstyle w[mM] = \\frac{1}{2\\pi} \\int_0^{2\\pi/M} \\sum_{k=0}^{M-1} W[e^{j(\\Omega - 2\\pi k/M)}] e^{jmM\\Omega} d\\Omega.\r\n",
        "$$\r\n",
        "\r\n",
        "\r\n",
        "With a final change of variable $\\Omega \\leftarrow M\\Omega$, we get\r\n",
        "$$\r\n",
        "    \\textstyle w[mM] = \\frac{1}{2\\pi} \\int_0^{2\\pi} \\frac{1}{M} \\sum_{k=0}^{M-1} W[e^{j(\\Omega - 2\\pi k)/M}] e^{jm\\Omega} d\\Omega.\r\n",
        "$$\r\n",
        "\r\n",
        "And by identifying (1) and (2), this yields:\r\n",
        "$$\r\n",
        "    Y(e^{j\\Omega}) = \\frac{1}{M} \\sum_{k=0}^{M-1} W[e^{j(\\Omega -2\\pi k)/M}]\r\n",
        "$$\r\n",
        "\r\n",
        "***\r\n",
        "\r\n",
        "In practice, here is the observed phenomenon when downsampling the spectrum above with factor $M=2$.\r\n",
        "\r\n",
        "<center> <img src=\"images/downsampling.png\" alt=\"\"  width=\"600\" height=\"300\"/> </center>"
      ],
      "metadata": {}
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "source": [
        "fs_down = 11025              # Target sampling frequency\r\n",
        "sound = allpath[14,29]       # Sound choice (default should be Shirping Birds)  \r\n",
        "\r\n",
        "print('Playing and showing data for : ', getname(sound) )\r\n",
        "x, fs = sf.read(sound)\r\n",
        "\r\n",
        "M = fs//fs_down # Downsampling factor\r\n",
        "print('Downsampling factor: ', M)\r\n",
        "\r\n",
        "### TO COMPLETE\r\n",
        "### Downsample \"audio\"\r\n",
        "x_naive_down = ...\r\n",
        "\r\n",
        "plot_audio(x,x_naive_down,fs,fs_down) # Function call"
      ],
      "outputs": [],
      "metadata": {
        "scrolled": true
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "Is your downsampling working properly in the time domain (verify using the zoom on the temporal signal)? What can you observe on the spectrum of the downsampled signal? How is this phenomenon named and what is its origin?\r\n",
        "\r\n",
        "In order to avoid it, the original signal should be low-pass filtered prior to downsampling (as presented in the mathematics above)."
      ],
      "metadata": {}
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "source": [
        "\"Low-pass filtering before downsampling\"\r\n",
        "N = 100 # number of taps\r\n",
        "taps = signal.firwin(numtaps=N, cutoff=fs_down/2, window='hamming', fs=fs)\r\n",
        "x_filt = np.convolve(x,taps,mode='full')\r\n",
        "\r\n",
        "### TO COMPLETE\r\n",
        "### Downsample ``audio_filt``\r\n",
        "x_filt_down = ...\r\n",
        "\r\n",
        "plot_audio(x,x_filt_down,fs,fs_down)"
      ],
      "outputs": [],
      "metadata": {}
    },
    {
      "cell_type": "markdown",
      "source": [
        "The obtained spectrum of the downsampled signal should not suffer from aliasing anymore. In fact, there is a built-in function in ``scipy.signal`` that performs the downsampling, including a low-pass filter: ``scipy.signal.resample``. Its docstring is:"
      ],
      "metadata": {}
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "source": [
        "help(signal.resample)"
      ],
      "outputs": [],
      "metadata": {}
    },
    {
      "cell_type": "markdown",
      "source": [
        "In the following, we use this function."
      ],
      "metadata": {}
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "source": [
        "y = signal.resample(x, int(len(x)/M))\r\n",
        "L = len(y)\r\n",
        "\r\n",
        "plot_audio(x,y,fs,fs_down)"
      ],
      "outputs": [],
      "metadata": {}
    },
    {
      "cell_type": "markdown",
      "source": [
        "Can you hear the differences between the downsampled versions of the audio signal and the original one?"
      ],
      "metadata": {}
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "source": [
        "sd.play(x, fs, blocksize=1024)\r\n",
        "#sd.play(x_naive_down, fs_down, blocksize=1024)\r\n",
        "#sd.play(x_filt_down, fs_down, blocksize=1024)\r\n",
        "#sd.play(y, fs_down, blocksize=1024)"
      ],
      "outputs": [],
      "metadata": {}
    },
    {
      "cell_type": "markdown",
      "source": [
        "You can also try sounds from different classes by running again the code above and changing the choice of the variable ``sound``.\r\n",
        "\r\n",
        "Now we are working with sound signals with same sampling frequency as for the project, we can go on."
      ],
      "metadata": {}
    },
    {
      "cell_type": "markdown",
      "source": [
        "## 2) Windowing and spectrogram computation\r\n",
        "\r\n",
        "A very intuitive way to represent an audio signal is with a time-frequency analysis. \r\n",
        "The spectrogram of a signal consists in applying an FFT on successive subpieces of it, and thus obtaining a spectral content evolving with time.\r\n",
        "\r\n",
        "Find an illustration of the idea here below.\r\n",
        "\r\n",
        "<center> <img src=\"images/melspecgram.jpg\" alt=\"\"  width=\"1000\" height=\"500\"/> </center>"
      ],
      "metadata": {}
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "source": [
        "Nft = 512 # Number of samples by FFT\r\n",
        "\r\n",
        "# Homemade computation of stft\r\n",
        "\"Crop the signal such that its length is a multiple of Nft\"\r\n",
        "y = y[:L-L%Nft]\r\n",
        "L = len(y)\r\n",
        "\"Reshape the signal with a piece for each row\"\r\n",
        "audiomat = np.reshape(y, (L//Nft,Nft))\r\n",
        "audioham = audiomat*np.hamming(Nft) # Windowing. Hamming, Hanning, Blackman,..\r\n",
        "z = np.reshape(audioham,-1) # y windowed by pieces\r\n",
        "\"FFT row by row\"\r\n",
        "stft = np.fft.fft(audioham, axis=1)\r\n",
        "stft = np.abs(stft[:,:Nft//2].T) # Taking only positive frequencies and computing the magnitude\r\n",
        "\r\n",
        "\"Library Librosa computing stft\"\r\n",
        "stft2 = librosa.stft(x, n_fft=Nft, hop_length=Nft, window='hamm', center='False') # without downsampling the signal\r\n",
        "stft4 = np.abs(librosa.stft(z, n_fft=Nft, hop_length=Nft, window='hamm', center=False))\r\n",
        "\r\n",
        "print(\"Note: You can eventually add a \\\"+1\\\" in the \\\"np.log\\\" to get positive dB. This will look differently.\")\r\n",
        "\r\n",
        "\"Plots\"\r\n",
        "fig = plt.figure(figsize=(9, 3))\r\n",
        "ax1 = fig.add_axes([0.0, 0.0, 0.42, 0.9])\r\n",
        "ax2 = fig.add_axes([0.54, 0.0, 0.42, 0.9])\r\n",
        "\r\n",
        "ax1.plot(np.arange(L)/fs_down, y, 'b', label='Original')\r\n",
        "ax1.plot(np.arange(L)/fs_down, z, 'r', label='Hamming windowed by pieces')\r\n",
        "ax1.set_xlabel('Time [s]')\r\n",
        "ax1.legend()\r\n",
        "\r\n",
        "plot_specgram(np.log(np.abs(stft2)), ax2, title='Specgram obtained with librosa.stft (full signal)', tf=len(x)/fs)\r\n",
        "plt.show()\r\n",
        "\r\n",
        "\"Comparing the spectrograms\" \r\n",
        "fig2 = plt.figure(figsize=(9, 3))\r\n",
        "ax3 = fig2.add_axes([0.0, 0.0, 0.42, 0.9])\r\n",
        "ax4 = fig2.add_axes([0.54, 0.0, 0.42, 0.9])\r\n",
        "plot_specgram(np.log(np.abs(stft)), ax3, title='Homemade specgram', tf=len(y)/fs_down)\r\n",
        "plot_specgram(np.log(np.abs(stft4)), ax4, title='Specgram obtained with librosa.stft', tf=len(y)/fs_down)\r\n",
        "plt.show()"
      ],
      "outputs": [],
      "metadata": {}
    },
    {
      "cell_type": "markdown",
      "source": [
        "What differences can you notice between the upper spectrogram and the two others at the bottom?\r\n",
        "\r\n",
        "<u> Remark:</u> Although the spectrograms obtained with the homemade version and the librosa library look similar, trying to compute and show their difference will not work. Indeed, there are subtle differences in the computation process of ``librosa.stft``, a slightly different window, frames which are computed are different as well. The priority for you is that the energy distribution in time and frequency remains similar, qualitatively speaking. This issue will be handled in the classification model.  "
      ],
      "metadata": {}
    },
    {
      "cell_type": "markdown",
      "source": [
        "## 3) From Hz to Melspectrogram\r\n",
        "\r\n",
        "Now we have done the major part of the job. But recall that this information will have to be transmitted wirelessly from your circuit board (transmitter) to a base station (receiver). It is thus good practive to try synthetizing a bit the content of this spectrogram. <br>\r\n",
        "\r\n",
        "A popular approach is to transform the frequency axis from Hz to Mel unit. The intuition behind this transformation is that the human ear will more easily distinguish between $100$ and $200$ Hz than between $3000$ and $3100$ Hz. So higher frequencies are more likely to be put together in very fewer coefficients. <br>\r\n",
        "\r\n",
        "This last step will thus consist in replacing each column of the spectrogram ``stft`` with size $N_{FT}$ by a shorter column with size $N_{mel} \\ll N_{FT}$. To do so, we will use an Hz to Mel (``Hz2Mel``) transformation matrix provided by ``librosa``, and apply a matrix multiplication for each column."
      ],
      "metadata": {}
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "source": [
        "Nmel = 20\r\n",
        "\r\n",
        "\"Obtain the Hz2Mel transformation matrix\" \r\n",
        "mels = librosa.filters.mel(sr=fs_down, n_fft=Nft, n_mels=Nmel)\r\n",
        "mels = mels[:,:-1]\r\n",
        "\r\n",
        "### TO COMPLETE\r\n",
        "### Normalize the mels matrix such that its maximum value is one.\r\n",
        "mels = mels\r\n",
        "\r\n",
        "\"Plot\"\r\n",
        "plt.figure(figsize=(5,4))\r\n",
        "plt.imshow(mels, aspect='auto')\r\n",
        "plt.gca().invert_yaxis()\r\n",
        "plt.colorbar()\r\n",
        "plt.title('Hz2Mel transformation matrix')\r\n",
        "plt.xlabel('$N_{FT}$')\r\n",
        "plt.ylabel('$N_{Mel}$')\r\n",
        "plt.show()"
      ],
      "outputs": [],
      "metadata": {}
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "source": [
        "\"Melspectrogram computation\"\r\n",
        "### TO COMPLETE\r\n",
        "###  Perform the matrix multiplication between the Hz2Mel matrix and stft.\r\n",
        "melspec = ...\r\n",
        "\r\n",
        "\"Plot\" \r\n",
        "fig = plt.figure(figsize=(9, 3))\r\n",
        "ax1 = fig.add_axes([0.0, 0.0, 0.42, 0.9])\r\n",
        "ax2 = fig.add_axes([0.54, 0.0, 0.42, 0.9])\r\n",
        "plot_specgram(np.log(np.abs(stft)), ax=ax1, title='Specgram', tf=len(y)/fs_down)\r\n",
        "plot_specgram(np.log(np.abs(melspec)), ax=ax2, is_mel=True, title='Melspecgram', tf=len(y)/fs_down)\r\n",
        "plt.show()"
      ],
      "outputs": [],
      "metadata": {}
    },
    {
      "cell_type": "markdown",
      "source": [
        "Do these two spectrogram look similar? :) <br> \r\n",
        "What is the gain in the number of coefficients?"
      ],
      "metadata": {}
    },
    {
      "cell_type": "markdown",
      "source": [
        "## 4) Creating black boxes\r\n",
        "\r\n",
        "Now you have seen how to make the computations. <br>\r\n",
        "A universal procedure consists in writing functions that will serve as working blocks and hide the computation details.\r\n",
        "We can then gradually increase the abstraction. <br>\r\n",
        "\r\n",
        "As any programmer should do, you are strongly encourage to ``fill your functions with a clear and concise docstring``. This will help you later this year when you will want to make improvements to some parts of your code."
      ],
      "metadata": {}
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "source": [
        "def resample(x, M=4):\r\n",
        "    \"\"\"   [description]\r\n",
        "    \r\n",
        "    Inputs\r\n",
        "      x:  [type, size, description]\r\n",
        "      M:  [type, size, description]\r\n",
        "\r\n",
        "    Outputs\r\n",
        "      y:   [type, size, description] \r\n",
        "\r\n",
        "    \"\"\"\r\n",
        "    \r\n",
        "    ### TO COMPLETE\r\n",
        "    return 0\r\n",
        "\r\n",
        "def specgram(y, Nft=512):\r\n",
        "    \"\"\"   [description]\r\n",
        "    \r\n",
        "    Inputs\r\n",
        "      y:  [type, size, description]\r\n",
        "      Nft:  [type, size, description]\r\n",
        "\r\n",
        "    Outputs\r\n",
        "      stft:   [type, size, description] \r\n",
        "\r\n",
        "    \"\"\"\r\n",
        "    \r\n",
        "    ### TO COMPLETE\r\n",
        "\r\n",
        "    return stft\r\n",
        "\r\n",
        "\r\n",
        "def melspecgram(x, Nmel=20, Nft=512, fs=44100, fs_down=11025):\r\n",
        "    \"\"\"   [description]\r\n",
        "    \r\n",
        "    Inputs\r\n",
        "      x:  [type, size, description]\r\n",
        "      Nmel:  [type, size, description]\r\n",
        "      Nft:  [type, size, description]\r\n",
        "      fs:  [type, size, description]\r\n",
        "      fs_down:  [type, size, description]\r\n",
        "\r\n",
        "    Outputs\r\n",
        "      melspec:   [type, size, description] \r\n",
        "\r\n",
        "    \"\"\"\r\n",
        "    \r\n",
        "    ### TO COMPLETE, using the functions resample() and specgram() defined above\r\n",
        "\r\n",
        "    return melspec"
      ],
      "outputs": [],
      "metadata": {}
    },
    {
      "cell_type": "markdown",
      "source": [
        "## 5) Show us your skills\r\n",
        "You are now encouraged to apply the functions you created above to sounds from at least 3 different classes. Observe their spectrograms and comment. Is it easy to differentiate sounds from the classes you chosed?"
      ],
      "metadata": {}
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "source": [
        "### TO COMPLETE\r\n",
        "### Choose 3 sounds from different classes to observe how their mel spectrograms differ\r\n",
        "sound1 = ...\r\n",
        "sound2 = ...\r\n",
        "sound3 = ...\r\n",
        "\r\n",
        "\r\n",
        "\"Compute the melspecgrams\"\r\n",
        "x1, _ = sf.read(sound1)\r\n",
        "x2, _ = sf.read(sound2)\r\n",
        "x3, _ = sf.read(sound3)\r\n",
        "melspec1 = melspecgram(x1)\r\n",
        "melspec2 = melspecgram(x2)\r\n",
        "melspec3 = melspecgram(x3)\r\n",
        "\r\n",
        "print(\"Note: Notice that here we added the \\\"+1\\\" for the visualization!\")\r\n",
        "\r\n",
        "\"Plot\" \r\n",
        "fig = plt.figure(figsize=(12, 3))\r\n",
        "ax1 = fig.add_axes([0.0, 0.0, 0.28, 0.9])\r\n",
        "ax2 = fig.add_axes([0.33, 0.0, 0.28, 0.9])\r\n",
        "ax3 = fig.add_axes([0.66, 0.0, 0.28, 0.9])\r\n",
        "plot_specgram(np.log(melspec1+1), ax=ax1, is_mel=True, title=getclass(sound1), tf=len(y)/fs_down)\r\n",
        "plot_specgram(np.log(melspec2+1), ax=ax2, is_mel=True, title=getclass(sound2), tf=len(y)/fs_down)\r\n",
        "plot_specgram(np.log(melspec3+1), ax=ax3, is_mel=True, title=getclass(sound3), tf=len(y)/fs_down)\r\n",
        "plt.show()"
      ],
      "outputs": [],
      "metadata": {}
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "source": [
        "### TO COMPLETE\r\n",
        "### Briefly comment what is intuitive for you in the content of these 3 spectrograms respectively with the corresponding classes.\r\n",
        "### How can you differentiate them?"
      ],
      "outputs": [],
      "metadata": {}
    }
  ],
  "metadata": {
    "interpreter": {
      "hash": "2c58dd59b1473525b9c2ce9ae8c4e1b85447983fecced33eb1cd38f3ca9c72e6"
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3.8.8 64-bit ('base': conda)"
    },
    "language_info": {
      "name": "python",
      "version": "3.8.8",
      "mimetype": "text/x-python",
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "pygments_lexer": "ipython3",
      "nbconvert_exporter": "python",
      "file_extension": ".py"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 2
}