fix a typo

bf6c6372 · Olivier Leblanc · 6e6f4dcb · bf6c6372
--- a/getting_started_student.ipynb
+++ b/getting_started_student.ipynb
@@ -302,7 +302,7 @@
        "\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",
+        "stft4 = np.abs(librosa.stft(y, 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",

 %% Cell type:markdown id: tags:

 # Hands-On 1: Audio Feature Extraction

 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>
 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.
 > Karol J. Piczak, 2015, "ESC: Dataset for Environmental Sound Classification"
 > https://doi.org/10.7910/DVN/YDEPUT, Harvard Dataverse, V2 ([Available here](https://dataverse.harvard.edu/dataset.xhtml?persistentId=doi:10.7910/DVN/YDEPUT))

 Useful functions to select, read and play the dataset sounds are provided in the ``utils`` folder. The following required packages should be installed:
 - soundfile
 - scipy
 - sounddevice
 - librosa
 - seaborn

 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>

 To ensure you are catching the content of this notebook, we leave you with an infinitesimal amount of **code to write**.

 You will find the zones to be briefly filled  with a ``### TO COMPLETE`` in the cells below.

 %% Cell type:code id: tags:

 ``` 
 '''
 !pip install -U soundfile
 !pip install -U sounddevice
 !pip install -U pygame
 !pip install -U librosa
 !pip install -U seaborn
 '''
 import numpy as np
 import matplotlib.pyplot as plt
 import soundfile as sf
 from scipy import signal
 import sounddevice as sd
 import librosa # For audio signal computations as MFCC
 from scipy.fftpack import dct

 "Self created functions"
 from utils_ import getclass, getname, gen_allpath, plot_audio, plot_specgram
 ```

 %% Cell type:markdown id: tags:

 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.

 %% Cell type:code id: tags:

 ``` 
 classnames, allpath = gen_allpath() # Note: this function contains an implicit input "folder", where you can change the path to ESC-50 dataset.
 print('The classes are : \n')
 for ind, elem in enumerate(classnames) :
    print(ind, elem)
 ```

 %% Cell type:markdown id: tags:

 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:

 %% Cell type:code id: tags:

 ``` 
 sound = allpath[3,29]
 x, fs = sf.read(sound)
 print(f"Playing a \"{getclass(sound)}\"")
 sd.play(x, fs)
 ```

 %% Cell type:markdown id: tags:

 We now ask you to complete the cells below.

 %% Cell type:markdown id: tags:

 ## 1) Resampling and filtering

 Most probably your circuit board will sample the analog audio signal at a frequency $f_s = 11025$ Hz. <br>
 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:
 - Rewrite a new dataset with the downsampled audio signals.
 - Downsample each audio which is read.

 We provide you with the second one.

 ***
 #### <u> The following derivations are not necessary for the rest of this notebook, but are still provided for the curious students... </u>

 Let us consider one original audio signal from the dataset and denote it $x[n]$, for $n=0,\dots,N-1$.

 The downsampled signal $y$ can be written as

 $$
    y[m] = w[mM],\quad \text{with}\ w[k] = (h \ast x)[k] = \sum_{n=-\infty}^{\infty} h[n]x[k-n],
 $$

 where $h$ is a discrete low-pass filter and $M$ is the downsampling factor, here $M=4$. <br>

 We can expand both $y$ and $w$ according to their Fourier series (DTFT) $Y$ and $W$, respectively, so that:

 $$
    y[m] = \frac{1}{2\pi} \int_0^{2\pi} Y(e^{j\Omega}) e^{jm\Omega} d\Omega \tag{1}
 $$
 $$
    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}
    $$

 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

 $$
   \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.
 $$


 With a final change of variable $\Omega \leftarrow M\Omega$, we get
 $$
    \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.
 $$

 And by identifying (1) and (2), this yields:
 $$
    Y(e^{j\Omega}) = \frac{1}{M} \sum_{k=0}^{M-1} W[e^{j(\Omega -2\pi k)/M}]
 $$

 ***

 In practice, here is the observed phenomenon when downsampling the spectrum above with factor $M=2$.

 <center> <img src="images/downsampling.png" alt=""  width="600" height="300"/> </center>

 %% Cell type:code id: tags:

 ``` 
 fs_down = 11025              # Target sampling frequency
 sound = allpath[14,29]       # Sound choice (default should be Shirping Birds)

 print('Playing and showing data for : ', getname(sound) )
 x, fs = sf.read(sound)

 M = fs//fs_down # Downsampling factor
 print('Downsampling factor: ', M)

 ### TO COMPLETE
 ### Downsample "audio"
 x_naive_down = ...

 plot_audio(x,x_naive_down,fs,fs_down) # Function call
 ```

 %% Cell type:markdown id: tags:

 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?

 In order to avoid it, the original signal should be low-pass filtered prior to downsampling (as presented in the mathematics above).

 %% Cell type:code id: tags:

 ``` 
 "Low-pass filtering before downsampling"
 N = 100 # number of taps
 taps = signal.firwin(numtaps=N, cutoff=fs_down/2, window='hamming', fs=fs)
 x_filt = np.convolve(x,taps,mode='full')

 ### TO COMPLETE
 ### Downsample ``audio_filt``
 x_filt_down = ...

 plot_audio(x,x_filt_down,fs,fs_down)
 ```

 %% Cell type:markdown id: tags:

 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:

 %% Cell type:code id: tags:

 ``` 
 help(signal.resample)
 ```

 %% Cell type:markdown id: tags:

 In the following, we use this function.

 %% Cell type:code id: tags:

 ``` 
 y = signal.resample(x, int(len(x)/M))
 L = len(y)

 plot_audio(x,y,fs,fs_down)
 ```

 %% Cell type:markdown id: tags:

 Can you hear the differences between the downsampled versions of the audio signal and the original one?

 %% Cell type:code id: tags:

 ``` 
 sd.play(x, fs, blocksize=1024)
 #sd.play(x_naive_down, fs_down, blocksize=1024)
 #sd.play(x_filt_down, fs_down, blocksize=1024)
 #sd.play(y, fs_down, blocksize=1024)
 ```

 %% Cell type:markdown id: tags:

 You can also try sounds from different classes by running again the code above and changing the choice of the variable ``sound``.

 Now we are working with sound signals with same sampling frequency as for the project, we can go on.

 %% Cell type:markdown id: tags:

 ## 2) Windowing and spectrogram computation

 A very intuitive way to represent an audio signal is with a time-frequency analysis.
 The spectrogram of a signal consists in applying an FFT on successive subpieces of it, and thus obtaining a spectral content evolving with time.

 Find an illustration of the idea here below.

 <center> <img src="images/melspecgram.jpg" alt=""  width="1000" height="500"/> </center>

 %% Cell type:code id: tags:

 ``` 
 Nft = 512 # Number of samples by FFT

 # Homemade computation of stft
 "Crop the signal such that its length is a multiple of Nft"
 y = y[:L-L%Nft]
 L = len(y)
 "Reshape the signal with a piece for each row"
 audiomat = np.reshape(y, (L//Nft,Nft))
 audioham = audiomat*np.hamming(Nft) # Windowing. Hamming, Hanning, Blackman,..
 z = np.reshape(audioham,-1) # y windowed by pieces
 "FFT row by row"
 stft = np.fft.fft(audioham, axis=1)
 stft = np.abs(stft[:,:Nft//2].T) # Taking only positive frequencies and computing the magnitude

 "Library Librosa computing stft"
 stft2 = librosa.stft(x, n_fft=Nft, hop_length=Nft, window='hamm', center='False') # without downsampling the signal
-stft4 = np.abs(librosa.stft(z, n_fft=Nft, hop_length=Nft, window='hamm', center=False))
+stft4 = np.abs(librosa.stft(y, n_fft=Nft, hop_length=Nft, window='hamm', center=False))

 print("Note: You can eventually add a \"+1\" in the \"np.log\" to get positive dB. This will look differently.")

 "Plots"
 fig = plt.figure(figsize=(9, 3))
 ax1 = fig.add_axes([0.0, 0.0, 0.42, 0.9])
 ax2 = fig.add_axes([0.54, 0.0, 0.42, 0.9])

 ax1.plot(np.arange(L)/fs_down, y, 'b', label='Original')
 ax1.plot(np.arange(L)/fs_down, z, 'r', label='Hamming windowed by pieces')
 ax1.set_xlabel('Time [s]')
 ax1.legend()

 plot_specgram(np.log(np.abs(stft2)), ax2, title='Specgram obtained with librosa.stft (full signal)', tf=len(x)/fs)
 plt.show()

 "Comparing the spectrograms"
 fig2 = plt.figure(figsize=(9, 3))
 ax3 = fig2.add_axes([0.0, 0.0, 0.42, 0.9])
 ax4 = fig2.add_axes([0.54, 0.0, 0.42, 0.9])
 plot_specgram(np.log(np.abs(stft)), ax3, title='Homemade specgram', tf=len(y)/fs_down)
 plot_specgram(np.log(np.abs(stft4)), ax4, title='Specgram obtained with librosa.stft', tf=len(y)/fs_down)
 plt.show()
 ```

 %% Cell type:markdown id: tags:

 What differences can you notice between the upper spectrogram and the two others at the bottom?

 <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.

 %% Cell type:markdown id: tags:

 ## 3) From Hz to Melspectrogram

 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>

 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>

 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.

 %% Cell type:code id: tags:

 ``` 
 Nmel = 20

 "Obtain the Hz2Mel transformation matrix"
 mels = librosa.filters.mel(sr=fs_down, n_fft=Nft, n_mels=Nmel)
 mels = mels[:,:-1]

 ### TO COMPLETE
 ### Normalize the mels matrix such that its maximum value is one.
 mels = mels

 "Plot"
 plt.figure(figsize=(5,4))
 plt.imshow(mels, aspect='auto')
 plt.gca().invert_yaxis()
 plt.colorbar()
 plt.title('Hz2Mel transformation matrix')
 plt.xlabel('$N_{FT}$')
 plt.ylabel('$N_{Mel}$')
 plt.show()
 ```

 %% Cell type:code id: tags:

 ``` 
 "Melspectrogram computation"
 ### TO COMPLETE
 ###  Perform the matrix multiplication between the Hz2Mel matrix and stft.
 melspec = ...

 "Plot"
 fig = plt.figure(figsize=(9, 3))
 ax1 = fig.add_axes([0.0, 0.0, 0.42, 0.9])
 ax2 = fig.add_axes([0.54, 0.0, 0.42, 0.9])
 plot_specgram(np.log(np.abs(stft)), ax=ax1, title='Specgram', tf=len(y)/fs_down)
 plot_specgram(np.log(np.abs(melspec)), ax=ax2, is_mel=True, title='Melspecgram', tf=len(y)/fs_down)
 plt.show()
 ```

 %% Cell type:markdown id: tags:

 Do these two spectrogram look similar? :) <br>
 What is the gain in the number of coefficients?

 %% Cell type:markdown id: tags:

 ## 4) Creating black boxes

 Now you have seen how to make the computations. <br>
 A universal procedure consists in writing functions that will serve as working blocks and hide the computation details.
 We can then gradually increase the abstraction. <br>

 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.

 %% Cell type:code id: tags:

 ``` 
 def resample(x, M=4):
    """   [description]

    Inputs
      x:  [type, size, description]
      M:  [type, size, description]

    Outputs
      y:   [type, size, description]

    """

    ### TO COMPLETE
    return 0

 def specgram(y, Nft=512):
    """   [description]

    Inputs
      y:  [type, size, description]
      Nft:  [type, size, description]

    Outputs
      stft:   [type, size, description]

    """

    ### TO COMPLETE

    return stft


 def melspecgram(x, Nmel=20, Nft=512, fs=44100, fs_down=11025):
    """   [description]

    Inputs
      x:  [type, size, description]
      Nmel:  [type, size, description]
      Nft:  [type, size, description]
      fs:  [type, size, description]
      fs_down:  [type, size, description]

    Outputs
      melspec:   [type, size, description]

    """

    ### TO COMPLETE, using the functions resample() and specgram() defined above

    return melspec
 ```

 %% Cell type:markdown id: tags:

 ## 5) Show us your skills
 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?

 %% Cell type:code id: tags:

 ``` 
 ### TO COMPLETE
 ### Choose 3 sounds from different classes to observe how their mel spectrograms differ
 sound1 = ...
 sound2 = ...
 sound3 = ...


 "Compute the melspecgrams"
 x1, _ = sf.read(sound1)
 x2, _ = sf.read(sound2)
 x3, _ = sf.read(sound3)
 melspec1 = melspecgram(x1)
 melspec2 = melspecgram(x2)
 melspec3 = melspecgram(x3)

 print("Note: Notice that here we added the \"+1\" for the visualization!")

 "Plot"
 fig = plt.figure(figsize=(12, 3))
 ax1 = fig.add_axes([0.0, 0.0, 0.28, 0.9])
 ax2 = fig.add_axes([0.33, 0.0, 0.28, 0.9])
 ax3 = fig.add_axes([0.66, 0.0, 0.28, 0.9])
 plot_specgram(np.log(melspec1+1), ax=ax1, is_mel=True, title=getclass(sound1), tf=len(y)/fs_down)
 plot_specgram(np.log(melspec2+1), ax=ax2, is_mel=True, title=getclass(sound2), tf=len(y)/fs_down)
 plot_specgram(np.log(melspec3+1), ax=ax3, is_mel=True, title=getclass(sound3), tf=len(y)/fs_down)
 plt.show()
 ```

 %% Cell type:code id: tags:

 ``` 
 ### TO COMPLETE
 ### Briefly comment what is intuitive for you in the content of these 3 spectrograms respectively with the corresponding classes.
 ### How can you differentiate them?
 ```