{"cells":[{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"# !pip install librosa\n# !pip install ipywidgets","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"This notebook covers the exploratory data analysis as well as different features that can be extracted from the libraries like [librosa](https://librosa.org/doc/latest/index.html). Ther has already been several notebooks based on EDA for this challenge. But since this is my first competition based on audio recognition, so I took the full opportunity to start from scratch and has been a great learning experience so far. Hope you'll like it."},{"metadata":{"trusted":true},"cell_type":"code","source":"# Necessary imports\n%matplotlib notebook\n\nimport os\nimport librosa\nimport numpy as np\nimport pandas as pd\nimport matplotlib as mpl\nimport matplotlib.pyplot as plt\nimport ipywidgets as widgets\nimport seaborn as sns\nfrom sklearn import preprocessing\n\n# plotting style\nsns.set_style(\"whitegrid\")\nsns.set_palette(\"muted\")\n\n# format plot options\nplt.rcParams[\"xtick.labelsize\"] = 14\nplt.rcParams[\"ytick.labelsize\"] = 14\nplt.rcParams[\"patch.force_edgecolor\"] = True","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## A little background\nA wave, in general, can be thought of as some kind of disturbance propagating in space and time. For light waves, these disturbances are nothing but changing electric and magnetic field with position as well as time. Similarly, for sound waves it is the change in the air pressure with time and position. \n\nMathematically, all of these waves (*at some fixed observation point*) can be represented as a function of time as $$sin(2 \\pi f t)$$ where $f$ is the frequency of the wave.\n\n## What does adding two sin waves do?\nOne property of this wave function is that it obeys principle of superposition meaning summing two or more sin waves with different amplitude will give us a sin wave with some frequency. This is precisely the case when people add more than one instrument and create music. These different instruments often have different frequencies and hence the resulting sin wave has frequency different from the individual frequencies. \n\nEach one of the audio files in this challenge contain sound waves which are a sum of lot of sin waves with different frequencies."},{"metadata":{"trusted":true,"_kg_hide-input":true,"_kg_hide-output":true},"cell_type":"code","source":"# help from: https://kapernikov.com/ipywidgets-with-matplotlib/\n# output = widgets.Output()\n\nfig, ax = plt.subplots(3, 1, figsize=(6, 10), constrained_layout=True)\n\n# generate X-values\nx = np.linspace(0, 2 * np.pi, 100)\n\ndef sin_fn1(x, w1, a1):\n    return a1 * np.sin(w1 * x)\n\ndef sin_fn2(x, w2, a2):\n    return a2 * np.sin(w2 * x)\n\ndef update_fn1(w1=1.0, a1=1.0):\n    [l.remove() for l in ax[0].lines]\n    ax[0].set_ylim([-4, 4])\n    ax[0].plot(x, (sin_fn1(x, w1, a1)), color='royalblue')\n    ax[0].set_title(r'$a_{1} \\sin(w_{1} x)$')\n\ndef update_fn2(w2=1.0, a2=1.0):\n    [l.remove() for l in ax[1].lines]\n    ax[1].set_ylim([-4, 4])\n    ax[1].plot(x, sin_fn2(x, w2, a2), color='crimson')\n    ax[1].set_title(r'$a_{2} \\sin(w_{2} x)$')\n\n@widgets.interact(w1=(0, 10, 1), w2=(0, 10, 1), a1=(0, 5, 1), a2=(0, 5, 1))\ndef update(w1=1.0, w2=1.0, a1=1.0, a2=1.0):\n    [l.remove() for l in ax[2].lines]\n    update_fn1(w1, a1)\n    update_fn2(w2, a2)\n    ax[2].set_ylim([-6, 6])\n    ax[2].plot(x, (sin_fn1(x, w1, a1) + sin_fn2(x, w2, a2)), color='forestgreen')\n    ax[2].set_title(r'$a_{1} \\sin(w_{1} x) + a_{2} \\sin (w_{2} x)$')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"_Unfortunately, the rendered notebook does not show the ipywidgets interactive plot. So had to hide the code as well as the plot._"},{"metadata":{"trusted":true},"cell_type":"code","source":"# define various paths\nroot_dir = '../input/rfcx-species-audio-detection'\ntrain_audio = os.path.join(root_dir, 'train')\ntest_audio = os.path.join(root_dir, 'test')\n\ntrain_tp = os.path.join(root_dir, 'train_tp.csv')\ntrain_fp = os.path.join(root_dir, 'train_fp.csv')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Some EDA"},{"metadata":{"trusted":true},"cell_type":"code","source":"# train true positive dataset\ntp_df = pd.read_csv(train_tp)\ntp_df.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"tp_df.info()","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"%matplotlib inline\n# explore the target column, i.e. species_id\nfig, ax = plt.subplots(nrows=2, ncols=1, figsize=(12, 10))\nsns.countplot(ax=ax[0], x='species_id', data=tp_df, alpha=0.6, color='navy')\nax[0].set_title('Distribution of Species for true positive data', fontsize=15)\n\nsns.countplot(ax=ax[1], x='species_id', hue='songtype_id', data=tp_df, alpha=0.7)\nax[1].set_title('Distribution of Species for true positive data w.r.t songtype', fontsize=15)\nfig.tight_layout()\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"As evident from above plots, species are almost uniformly distributed in the true positive dataset. Also, other than species with `species_id` 16, 17, and 23, all the species have just one `songtype` which basically is the type of sounds produced by a given species."},{"metadata":{"trusted":true},"cell_type":"code","source":"# train false positive dataset\nfp_df = pd.read_csv(train_fp)\nfp_df.head()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"fp_df.info()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# explore the target column, i.e. species_id\nfig, ax = plt.subplots(nrows=2, ncols=1, figsize=(12, 10))\nsns.countplot(ax=ax[0], x='species_id', data=fp_df, alpha=0.6, color='navy')\nax[0].set_title('Distribution of Species for true positive data', fontsize=15)\n\nsns.countplot(ax=ax[1], x='species_id', hue='songtype_id', data=fp_df, alpha=0.7)\nax[1].set_title('Distribution of Species for true positive data w.r.t songtype', fontsize=15)\nfig.tight_layout()\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Although the labels in false positive dataset are verified by the experts to _not_ contain the [flagged species/songtype](https://www.kaggle.com/c/rfcx-species-audio-detection/discussion/197866) but still the overall trend stays more or less the same as true positive."},{"metadata":{"trusted":true},"cell_type":"code","source":"# playing audio\nfrom IPython.display import Audio","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# select some random samples from train set\ndef get_random_samples(data_path, num_samples=1):\n    data_list = os.listdir(data_path)\n    indices = np.random.choice(len(data_list), num_samples)\n    sample_audio = [os.path.join(data_path, data_list[idx]) for idx in indices]\n    \n    return sample_audio","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# get random train samples\naudio_samples_train = get_random_samples(train_audio)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"Audio(audio_samples_train[0])","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"from librosa import display\ntrain_amp, train_sr = librosa.load(audio_samples_train[0])\nrecord_id = os.path.basename(audio_samples_train[0]).split('.')[0]\n\nprint(f'Recording ID: {record_id}')\nprint(f'Total number of samples in the recording: {len(train_amp)}')\nprint(f'Sampling rate of the recording: {train_sr}')\n\nprint('Checking the data corresponding to the recording ID')\nrecord_df = tp_df[tp_df[\"recording_id\"] == record_id]\n\nif record_df.empty:\n    record_df = fp_df[fp_df[\"recording_id\"] == record_id]\nprint(record_df)\n\nprint('Now plotting the recording waveform...')\nplt.figure(figsize=(14, 6))\nlibrosa.display.waveplot(y=train_amp, sr=train_sr, color='navy', alpha=0.5)\nplt.xlabel('Time (seconds)-->')\nplt.ylabel('<-- Amplitude -->')\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"`train_amp` is just an array of amplitudes whose length gives total number of samples whereas `train_sr` is the sampling rate.\n\n__Sampling rate__ or __Sampling frequency__ is the rate at which we are capturing amplitudes. In other words, it is just the number of data points recorded per second. These amplitudes can be electric or magnetic fileds in the case of light waves, current or voltage in the case of digital signals and pressure or displacement values for the case of sound waves.\n\nThe above plot gives the time-domain representation of the signal. Higher amplitude just corresponds to loudness in the signal and the points of zero amplitude represent silence. Other than giving information about the variation of the loudness of the audio signal with time, time-domain does not convey other important information about the signal.\n\nWe can rather try to understand the signal in frequency domain which reveals much more information about the signal. The reason for that is the sound that we are hearing results from the superposition of many different audio signals with different frequencies. \n\nThis is our motivation for working in the frequency domain since we can decompose the wave into its constituent frequencies. The mathematical tools used for this is called __Fourier Transform__. Since our signal is composed of several discrete samples, we will use Fast Fourier Transform instead of usual Fourier Transform.\n\nLibrosa has a special method which performs Fourier transform, it's called `stft` (Short Time Fourier Transform). Using `stft` we can generate a plot which shows the signal in frequency-time space, this plot is commonly referred to as __Spectrogram__."},{"metadata":{},"cell_type":"markdown","source":"## Spectrograms\n\nUsing Librosa, we can easily create spectrograms using Short Time Fourier Transform. Following two plots show spectrograms for train audio on both Hertz (linear) as well as log scale."},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"train_spec = librosa.stft(train_amp)\nprint(train_spec.shape)\n\n# convert amplitude into decibel scale\ntrain_db = librosa.amplitude_to_db(abs(train_spec))\nprint(train_db.shape)\n\nfig, ax = plt.subplots(nrows=2, ncols=1, figsize=(14, 12), sharex=True)\nimg = librosa.display.specshow(train_db, sr=train_sr, x_axis='time', \n                         y_axis='hz', ax=ax[0], cmap='twilight')\nax[0].set(title='Linear-frequency power spectrogram')\nax[0].label_outer()\n\nlibrosa.display.specshow(train_db, sr=train_sr, x_axis='time', \n                         y_axis='log', ax=ax[1], cmap='twilight')\nax[1].set(title='Log-frequency power spectrogram')\nax[1].label_outer()\nfig.colorbar(img, ax=ax, format='%+2.f dB')\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Clearly, frequencies in`log` scale show much more detail than the linear scale. However, in the linear scale, the perceptual distance between 300 Hz and 500 Hz might not seem equal to the distance between 11400 Hz and 11600Hz to human ears even though the difference is the same. So it is quite common so describe the spectrogram in __mel__ scale which is similar to the `log` scale and it is defined as \n$$\nf_{m} = 2595\\: \\mathrm{log}_{10}\\left( 1 + \\frac{f}{700}\\right)\n$$\nMel scale is just a result of the above non-linear transformation. It renders the frequencies, that are at equal distances from each other, also felt by humans as if they are at equal distances. This is because human ear does not perceive frequencies on a linear scale.\n\n## Mel-Spectrograms"},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"plt.figure(figsize=(14, 6))\nlibrosa.display.specshow(train_db, sr=train_sr, x_axis='time', \n                         y_axis='mel', cmap='twilight')\nplt.colorbar(format='%+2.0f dB')\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Spectral Centroids\nTo quote the wikipedia, spectral centroid gives us the \"impression of brightness of sound\". Just like in physics where center of mass gives the location of the point where the whole mass of the body can be thought to be concentrated, spectral centroid can be visualized as the as the center of mass of the frequencies in a sound spectrum (typically a spectrogram). Mathematically, it is the weighted sum of the frequencies present in the signal."},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"spectral_centroids = librosa.feature.spectral_centroid(train_amp, sr=train_sr)[0]\nprint(spectral_centroids.shape)\n\n# extract the time and frame indices\nframes = range(len(spectral_centroids))\nt = librosa.frames_to_time(frames)\n\n\nplt.figure(figsize=(14, 6))\nlibrosa.display.specshow(train_db, sr=train_sr, x_axis='time', \n                         y_axis='log', cmap='twilight')\n\nplt.plot(t, spectral_centroids, color='yellow', \n         label='Spectral Centroid')\nplt.title(\"Spectral Centroids for Train sample\")\nplt.colorbar(format='%+2.0f dB')\nplt.legend(loc='upper right', fontsize=16, facecolor='gray')\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Spectral Roll-off\nWe will calculate some more features about the dataset. As described in Librosa documentation, it corresponds to that frequency for which a given spectrogram frequency bin has at least a given percent of spectrum energy (called roll_percent, 0.85 by default) in this stft frame (2584 in above spectrogram) is contained in this bin and the bins below. In simpler words, the roll-off frequency is defined as the frequency under which some percentage (cutoff) of the total energy of the spectrum is contained. Setting the roll percent to approximately 1 (or 0) can be used to find the maximum (or minimum) frequency."},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"plt.figure(figsize=(14, 6))\nspectral_rolloff_full = librosa.feature.spectral_rolloff(train_amp, sr=train_sr,\n                                                         roll_percent=0.95)[0]\nspectral_rolloff_empty = librosa.feature.spectral_rolloff(train_amp, sr=train_sr, \n                                                          roll_percent=0.01)[0]\nlibrosa.display.specshow(train_db, sr=train_sr, x_axis='time', \n                         y_axis='log', cmap='twilight')\nplt.plot(t, spectral_rolloff_full, color='white', \n         label='Roll-off frequency (0.95)')\nplt.plot(t, spectral_rolloff_empty, color='yellow', \n         label='Roll-off frequency (0.01)')\nplt.title(\"Spectral Roll-off for Train sample\")\nplt.colorbar(format='%+2.0f dB')\nplt.legend(loc='lower right', fontsize=16, facecolor='gray')\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## MFCC (Mel-Frequency Cepstral Coefficients)\n\nMFCCs of a signal are a small set of features (usually 15-35) which are used to describe the overall shape of the spectral envelope for each time frame. Consider the signal for each time frame as a histogram showing the distribution of frquencies. MFCCs basically are the bins of that histogram. It is explained in more detail [here](http://www.speech.cs.cmu.edu/15-492/slides/03_mfcc.pdf)."},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"mfccs = librosa.feature.mfcc(y=train_amp, sr=train_sr)\nprint(mfccs.shape)\n\n# display the scaled mfccs\nplt.figure(figsize=(14, 6))\nlibrosa.display.specshow(mfccs, sr=train_sr, x_axis='time', \n                         cmap='twilight')\nplt.title('MFCCs', fontsize=15)\nplt.colorbar()\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Chroma STFT\nChroma features is a set of 12 elements feature vectors. This 12 element feature corresponds to the energy contained in each pitch class  (7 major notes + 5 sharps). This is mainly used to detect similarity between music and ASR."},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"def plot_feature(signal, y_axis='chroma', x_axis='time', \n                 title='Feature space', format=None):\n    plt.figure(figsize=(14,6))\n    librosa.display.specshow(signal, y_axis=y_axis, \n                             x_axis=x_axis, cmap='twilight')\n    plt.colorbar(format=format)\n    plt.title(title)\n    return plt.show()","execution_count":null,"outputs":[]},{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"# chroma stft\n# get the energy spectrum\nsignal = np.abs(librosa.stft(train_amp))\nprint(signal.shape)\n\nchroma_stft = librosa.feature.chroma_stft(S=signal, sr=train_sr)\nprint(chroma_stft.shape)\n\n# now plot the feature\nplot_feature(chroma_stft, title='Chroma STFT')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Chroma-CQT\nLike Mel scale, the constant Q-transform uses a logarithm scale for the frequencies."},{"metadata":{"trusted":true},"cell_type":"code","source":"# chroma cqt\nchroma_cqt = librosa.feature.chroma_cqt(y=train_amp, sr=train_sr)\nprint(chroma_cqt.shape)\n\n# now plot the feature\nplot_feature(chroma_cqt, title='Chroma CQT')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Chroma-CENS\n__Chroma Energy Normalized Statistics__ is used to smoothen local deviations in tempo, pitch, etc by taking statistics over large time windows. It is also used for audio matching and similarity."},{"metadata":{"trusted":true},"cell_type":"code","source":"# chroma cens\nchroma_cens = librosa.feature.chroma_cens(y=train_amp, sr=train_sr)\nprint(chroma_cens.shape)\n\n# now plot the feature\nplot_feature(chroma_cens, title='Chroma CENS')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Poly-Features\nThes features are just the polynomial coefficients obtained by fitting polynomials to the spectral envelope in each of the time-frame of the spectrogram."},{"metadata":{"trusted":true},"cell_type":"code","source":"# poly features\nsignal_poly = librosa.feature.poly_features(S=signal, sr=train_sr, \n                                            order=2)\nprint(signal_poly.shape)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"## Tonnetz Features\nThis method of feature extraction projects the chroma features onto a 6-dimensional basis corresponding to perfect-fifths and major and minor third as two-dimensional coordinates. Like spectral centroid, it computes the tonal centroids represented as the 6-dimensional basis."},{"metadata":{"trusted":true},"cell_type":"code","source":"# tonnetz features\nsignal_tonnetz = librosa.feature.tonnetz(y=train_amp, sr=train_sr, \n                                         chroma=chroma_cqt)\nprint(signal_tonnetz.shape)\n\nplot_feature(signal_tonnetz, title='Tonal Centroids', y_axis='tonnetz')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"markdown","source":"All of these features can be aggregated over the given frame and can be used in linear and tree-based models to provide a baseline for more sophisticated methods. "},{"metadata":{},"cell_type":"markdown","source":"### References:\n1. Librosa [documentation](https://librosa.org/doc/main/feature.html)\n2. https://musicinformationretrieval.com/index.html\n3. https://haythamfayek.com/2016/04/21/speech-processing-for-machine-learning.html\n4. https://www.kdnuggets.com/2020/02/audio-data-analysis-deep-learning-python-part-1.html"},{"metadata":{"trusted":true},"cell_type":"code","source":"","execution_count":null,"outputs":[]}],"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat":4,"nbformat_minor":4}