{"cells":[{"metadata":{},"cell_type":"markdown","source":"# Introduction","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"This is an example notebook for testing noise filtering methods on Birds Dataset.\nThe focus is to start to explore some techniques of preprocessing that can be used to improve bird detection models.\nThe problem is that the Cornell Birdcall Identification Dataset have many sounds of low quality, with high level background sounds and noises. Thus, I believe that applying some filtering method for noise reduction can improve the classification task.\n\nCommon detection pipelines consist of:\n* Preprocessing: To read a audio file, optionally apply signal filtering techniques, and perform feature extraction (e.g. mfccs);\n* Trainning a classification model based on features (TO DO);\n* Evaluation: To test the trainned models over a split of the dataset (TO DO).\n\nThe preprocess methods you can see here are:\n* Traditional log mel-spectogram;\n* High-Pass Filtering: Reduces low frequencies, once bird sound are commonly present on high frequencies;\n* Per-channel energy normalization (PCEN): Technique for automatic gain control, followed by nonlinear compression;\n* Spectral Gating: Common strategy for denoising music by gating the signal only on high level sounds.\n\nPlease, use this notebook as a didactical one. \nIf you enjoy it, please, leave your upvote and comments.\n\nThanks!","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"# Imports","execution_count":null},{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":true},"cell_type":"code","source":"# feature extractoring and preprocessing data\nimport librosa\nimport librosa.display\nimport pandas as pd\nimport numpy as np\nimport scipy.signal\nimport matplotlib.pyplot as plt\n%matplotlib inline\nfrom PIL import Image\nfrom pathlib import Path\nfrom pylab import rcParams\nrcParams['figure.figsize'] = 14, 6\n\nimport csv\n# Preprocessing\nfrom sklearn.model_selection import train_test_split\nfrom sklearn.preprocessing import LabelEncoder, StandardScaler\n#Reports\nfrom sklearn.metrics import classification_report, confusion_matrix\n\nimport warnings\nwarnings.filterwarnings('ignore')\n","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Preprocessing methods","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"## Reading some audio samples","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"Let's see some preprocessing techniques over the example_test_audio dataset.","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"sr = 16000\ne_file1 = '../input/birdsong-recognition/example_test_audio/BLKFR-10-CPL_20190611_093000.pt540.mp3'\ne_file2 = '../input/birdsong-recognition/example_test_audio/ORANGE-7-CAP_20190606_093000.pt623.mp3'\n\n# 10 seconds of each file\ny1,sr = librosa.load(e_file1, mono=True, sr=sr, offset=0, duration=10)\ny2,sr = librosa.load(e_file2, mono=True, sr=sr, offset=0, duration=10)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Listen to them!","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"from IPython.display import Audio, IFrame, display\n\ndisplay(Audio(y1,rate=sr))\ndisplay(Audio(y2,rate=sr))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"As you can hear, the first audio presents a high level background noise, and birds seems far from the mic. In the second audio, bird sounds are much more distinguished from the other noises. We can say that the second audio presents a better SNR (signal-noise ratio).","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"librosa.display.waveplot(y1,sr=sr, x_axis='time');","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"librosa.display.waveplot(y2,sr=sr, x_axis='time');","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"If you look to both signal waves, you can see that the first sound presents a lower difference between the background level and some sound events, while this difference is much higher in the second sound. If you do it in sync with the audio, we can see that some bird calls does not appear from the average sound level.","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"## Logmel-spectogram","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"A very common preprocessing technique in audio detection applications is to transform audios to its log mel-spectogram representation.\nSome concepts here: https://en.wikipedia.org/wiki/Mel-frequency_cepstrum","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"S1 = librosa.feature.melspectrogram(y=y1, sr=sr, n_mels=64)\nD1 = librosa.power_to_db(S1, ref=np.max)\nlibrosa.display.specshow(D1, x_axis='time', y_axis='mel');","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"S2 = librosa.feature.melspectrogram(y=y2, sr=sr, n_mels=64)\nD2 = librosa.power_to_db(S2, ref=np.max)\nlibrosa.display.specshow(D2, x_axis='time', y_axis='mel');","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"The differences become very clear in the log mel spectogram. In the fist case, you can see a lot of artefacts on low frequencies (not birds), and the birds are in levels below the background noises. Besides, background noises are higher in frequencies below 2 kHz.","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"## Filtering low-frequencies","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"As we noticed, low frequencies does not contribute to bird sounds, a first idea is to remove these low frequencies. A high pass filter helps in this task. Reference: https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.lfilter.html","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"from scipy import signal\nimport random\n\n\ndef f_high(y,sr):\n    b,a = signal.butter(10, 2000/(sr/2), btype='highpass')\n    yf = signal.lfilter(b,a,y)\n    return yf\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"yf1 = f_high(y1, sr)\nyf2 = f_high(y2, sr)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Let's see...","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"librosa.display.waveplot(y1,sr=sr, x_axis='time');\nlibrosa.display.waveplot(yf1,sr=sr, x_axis='time');","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"librosa.display.waveplot(y2,sr=sr, x_axis='time');\nlibrosa.display.waveplot(yf2,sr=sr, x_axis='time');","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"Sf1 = librosa.feature.melspectrogram(y=yf1, sr=sr, n_mels=64)\nDf1 = librosa.power_to_db(Sf1, ref=np.max)\nlibrosa.display.specshow(Df1, x_axis='time', y_axis='mel');","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"Sf2 = librosa.feature.melspectrogram(y=yf2, sr=sr, n_mels=64)\nDf2 = librosa.power_to_db(Sf2, ref=np.max)\nlibrosa.display.specshow(Df2, x_axis='time', y_axis='mel');","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"display(Audio(yf1,rate=sr))\ndisplay(Audio(yf2,rate=sr))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"In both cases, the filter helped to isolate the interesting frequencies. The second audio is in a very good quality for distincting the birds. But the first audio also have high level noises in bird frequencies. Simply removing these frequencies can make us loose important information.","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"## PCEN","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"PCEN has become a very useful strategy for acoustic event detection, and it has shown to perform better in such tasks as a frontend. Its idea is to perform non-linear compression on time-frequency channels.\n\nI am using the example shown here: https://librosa.org/doc/latest/generated/librosa.pcen.html?highlight=pcen#librosa.pcen","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"Dp1 = librosa.pcen(S1 * (2**31), sr=sr, gain=1.1, hop_length=512, bias=2, power=0.5, time_constant=0.8, eps=1e-06, max_size=2)\nDp2 = librosa.pcen(S2 * (2**31), sr=sr, gain=1.1, hop_length=512, bias=2, power=0.5, time_constant=0.8, eps=1e-06, max_size=2)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"librosa.display.specshow(Dp1, x_axis='time', y_axis='mel');","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"librosa.display.specshow(Dp2, x_axis='time', y_axis='mel');","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"yp1 = librosa.feature.inverse.mel_to_audio(Dp1)\nyp2 = librosa.feature.inverse.mel_to_audio(Dp2)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"librosa.display.waveplot(yp1,sr=sr, x_axis='time');\nlibrosa.display.waveplot(y1,sr=sr, x_axis='time');","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"librosa.display.waveplot(yp2,sr=sr, x_axis='time');\nlibrosa.display.waveplot(y2,sr=sr, x_axis='time');","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"display(Audio(yp1,rate=sr))\ndisplay(Audio(yp2,rate=sr))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"Uowwwww! It looks very promising!\nThis method does not eliminate all background, but the bird signal shapes are much more visible in the spectograms.\nThis method works on time-frequency representarions of the sounds, so I perform signal reconstitution from spectogram to audio to gives us a good idea on how it worked. This recostitution is not perfect, some artefacts are inserted.","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"## Spectral Gating","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"This is also a technique for noise reduction based on gates that monitor audio level. It is commonly used in music industry, and present in tools like Audacity (https://wiki.audacityteam.org/wiki/How_Audacity_Noise_Reduction_Works).\nI reproduced here the code made available by Tim Sainburg in his github (https://github.com/timsainb/noisereduce).","execution_count":null},{"metadata":{"trusted":true},"cell_type":"code","source":"import time\nfrom datetime import timedelta as td\n\n\ndef _stft(y, n_fft, hop_length, win_length):\n    return librosa.stft(y=y, n_fft=n_fft, hop_length=hop_length, win_length=win_length)\n\n\ndef _istft(y, hop_length, win_length):\n    return librosa.istft(y, hop_length, win_length)\n\n\ndef _amp_to_db(x):\n    return librosa.core.amplitude_to_db(x, ref=1.0, amin=1e-20, top_db=80.0)\n\n\ndef _db_to_amp(x,):\n    return librosa.core.db_to_amplitude(x, ref=1.0)\n\n\ndef plot_spectrogram(signal, title):\n    fig, ax = plt.subplots(figsize=(20, 4))\n    cax = ax.matshow(\n        signal,\n        origin=\"lower\",\n        aspect=\"auto\",\n        cmap=plt.cm.seismic,\n        vmin=-1 * np.max(np.abs(signal)),\n        vmax=np.max(np.abs(signal)),\n    )\n    fig.colorbar(cax)\n    ax.set_title(title)\n    plt.tight_layout()\n    plt.show()\n\n\ndef plot_statistics_and_filter(\n    mean_freq_noise, std_freq_noise, noise_thresh, smoothing_filter\n):\n    fig, ax = plt.subplots(ncols=2, figsize=(20, 4))\n    plt_mean, = ax[0].plot(mean_freq_noise, label=\"Mean power of noise\")\n    plt_std, = ax[0].plot(std_freq_noise, label=\"Std. power of noise\")\n    plt_std, = ax[0].plot(noise_thresh, label=\"Noise threshold (by frequency)\")\n    ax[0].set_title(\"Threshold for mask\")\n    ax[0].legend()\n    cax = ax[1].matshow(smoothing_filter, origin=\"lower\")\n    fig.colorbar(cax)\n    ax[1].set_title(\"Filter for smoothing Mask\")\n    plt.show()\n\ndef removeNoise(\n    audio_clip,\n    noise_clip,\n    n_grad_freq=2,\n    n_grad_time=4,\n    n_fft=2048,\n    win_length=2048,\n    hop_length=512,\n    n_std_thresh=1.5,\n    prop_decrease=1.0,\n    verbose=False,\n    visual=False,\n):\n    \"\"\"Remove noise from audio based upon a clip containing only noise\n\n    Args:\n        audio_clip (array): The first parameter.\n        noise_clip (array): The second parameter.\n        n_grad_freq (int): how many frequency channels to smooth over with the mask.\n        n_grad_time (int): how many time channels to smooth over with the mask.\n        n_fft (int): number audio of frames between STFT columns.\n        win_length (int): Each frame of audio is windowed by `window()`. The window will be of length `win_length` and then padded with zeros to match `n_fft`..\n        hop_length (int):number audio of frames between STFT columns.\n        n_std_thresh (int): how many standard deviations louder than the mean dB of the noise (at each frequency level) to be considered signal\n        prop_decrease (float): To what extent should you decrease noise (1 = all, 0 = none)\n        visual (bool): Whether to plot the steps of the algorithm\n\n    Returns:\n        array: The recovered signal with noise subtracted\n\n    \"\"\"\n    if verbose:\n        start = time.time()\n    # STFT over noise\n    noise_stft = _stft(noise_clip, n_fft, hop_length, win_length)\n    noise_stft_db = _amp_to_db(np.abs(noise_stft))  # convert to dB\n    # Calculate statistics over noise\n    mean_freq_noise = np.mean(noise_stft_db, axis=1)\n    std_freq_noise = np.std(noise_stft_db, axis=1)\n    noise_thresh = mean_freq_noise + std_freq_noise * n_std_thresh\n    if verbose:\n        print(\"STFT on noise:\", td(seconds=time.time() - start))\n        start = time.time()\n    # STFT over signal\n    if verbose:\n        start = time.time()\n    sig_stft = _stft(audio_clip, n_fft, hop_length, win_length)\n    sig_stft_db = _amp_to_db(np.abs(sig_stft))\n    if verbose:\n        print(\"STFT on signal:\", td(seconds=time.time() - start))\n        start = time.time()\n    # Calculate value to mask dB to\n    mask_gain_dB = np.min(_amp_to_db(np.abs(sig_stft)))\n    #print(noise_thresh, mask_gain_dB)\n    # Create a smoothing filter for the mask in time and frequency\n    smoothing_filter = np.outer(\n        np.concatenate(\n            [\n                np.linspace(0, 1, n_grad_freq + 1, endpoint=False),\n                np.linspace(1, 0, n_grad_freq + 2),\n            ]\n        )[1:-1],\n        np.concatenate(\n            [\n                np.linspace(0, 1, n_grad_time + 1, endpoint=False),\n                np.linspace(1, 0, n_grad_time + 2),\n            ]\n        )[1:-1],\n    )\n    smoothing_filter = smoothing_filter / np.sum(smoothing_filter)\n    # calculate the threshold for each frequency/time bin\n    db_thresh = np.repeat(\n        np.reshape(noise_thresh, [1, len(mean_freq_noise)]),\n        np.shape(sig_stft_db)[1],\n        axis=0,\n    ).T\n    # mask if the signal is above the threshold\n    sig_mask = sig_stft_db < db_thresh\n    if verbose:\n        print(\"Masking:\", td(seconds=time.time() - start))\n        start = time.time()\n    # convolve the mask with a smoothing filter\n    sig_mask = scipy.signal.fftconvolve(sig_mask, smoothing_filter, mode=\"same\")\n    sig_mask = sig_mask * prop_decrease\n    if verbose:\n        print(\"Mask convolution:\", td(seconds=time.time() - start))\n        start = time.time()\n    # mask the signal\n    sig_stft_db_masked = (\n        sig_stft_db * (1 - sig_mask)\n        + np.ones(np.shape(mask_gain_dB)) * mask_gain_dB * sig_mask\n    )  # mask real\n    sig_imag_masked = np.imag(sig_stft) * (1 - sig_mask)\n    sig_stft_amp = (_db_to_amp(sig_stft_db_masked) * np.sign(sig_stft)) + (\n        1j * sig_imag_masked\n    )\n    if verbose:\n        print(\"Mask application:\", td(seconds=time.time() - start))\n        start = time.time()\n    # recover the signal\n    recovered_signal = _istft(sig_stft_amp, hop_length, win_length)\n    recovered_spec = _amp_to_db(\n        np.abs(_stft(recovered_signal, n_fft, hop_length, win_length))\n    )\n    if verbose:\n        print(\"Signal recovery:\", td(seconds=time.time() - start))\n    if visual:\n        plot_spectrogram(noise_stft_db, title=\"Noise\")\n    if visual:\n        plot_statistics_and_filter(\n            mean_freq_noise, std_freq_noise, noise_thresh, smoothing_filter\n        )\n    if visual:\n        plot_spectrogram(sig_stft_db, title=\"Signal\")\n    if visual:\n        plot_spectrogram(sig_mask, title=\"Mask applied\")\n    if visual:\n        plot_spectrogram(sig_stft_db_masked, title=\"Masked signal\")\n    if visual:\n        plot_spectrogram(recovered_spec, title=\"Recovered spectrogram\")\n    return recovered_signal","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"noise1 = y1[5*sr:6*sr]\nyg1 = removeNoise(audio_clip=y1, noise_clip=noise1,\n    n_grad_freq=2,\n    n_grad_time=4,\n    n_fft=2048,\n    win_length=2048,\n    hop_length=512,\n    n_std_thresh=1.5,\n    prop_decrease=1.0,\n    verbose=False,\n    visual=False)\nnoise2 = y2[0:1*sr]\nyg2 = removeNoise(audio_clip=y2, noise_clip=noise2,\n    n_grad_freq=2,\n    n_grad_time=4,\n    n_fft=2048,\n    win_length=2048,\n    hop_length=512,\n    n_std_thresh=2.5,\n    prop_decrease=1.0,\n    verbose=False,\n    visual=False)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"librosa.display.waveplot(y1,sr=sr, x_axis='time');\nlibrosa.display.waveplot(yg1,sr=sr, x_axis='time');","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"librosa.display.waveplot(y2,sr=sr, x_axis='time');\nlibrosa.display.waveplot(yg2,sr=sr, x_axis='time');","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"Sg1 = librosa.feature.melspectrogram(y=yg1, sr=sr, n_mels=64)\nDg1 = librosa.power_to_db(Sg1, ref=np.max)\nlibrosa.display.specshow(Dg1, x_axis='time', y_axis='mel');","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"Sg2 = librosa.feature.melspectrogram(y=yg2, sr=sr, n_mels=64)\nDg2 = librosa.power_to_db(Sg2, ref=np.max)\nlibrosa.display.specshow(Dg2, x_axis='time', y_axis='mel');","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"display(Audio(yg1,rate=sr))\ndisplay(Audio(yg2,rate=sr))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"It didn't work so well on the noisy sound, but it was very good on the high SNR sample. Maybe some parameters can be adjusted to improve the results.","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"# Combining the methods","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"Yes, I know what you are thinking. What if I combine these methods? From this code it is very easy to do such tests.\n\nAlso, there are lots of parameter to try!\n\nIt is up to you.","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"# Impact on training and classification","execution_count":null},{"metadata":{},"cell_type":"markdown","source":"I am still testing these preprocessing methods on my models.\n\nBy now, I have a very initial model here: https://www.kaggle.com/mauriciofigueiredo/intro-to-filtering-process-model-submitting\n\nAs you can see, it is a very basic model based on KNN over mfcc features. I combined some of the above strategies and applied on a subset of bird classes. The accuracy increase was around 10% over a split of  the selected data.","execution_count":null}],"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}