{"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_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"import 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\nimport warnings\nwarnings.filterwarnings('ignore')\n","metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","execution":{"iopub.status.busy":"2022-09-12T06:50:01.748998Z","iopub.execute_input":"2022-09-12T06:50:01.749369Z","iopub.status.idle":"2022-09-12T06:50:03.343093Z","shell.execute_reply.started":"2022-09-12T06:50:01.749337Z","shell.execute_reply":"2022-09-12T06:50:03.341948Z"},"trusted":true},"execution_count":null,"outputs":[]},{"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\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)","metadata":{"execution":{"iopub.status.busy":"2022-09-12T06:50:03.347380Z","iopub.execute_input":"2022-09-12T06:50:03.347665Z","iopub.status.idle":"2022-09-12T06:50:04.897751Z","shell.execute_reply.started":"2022-09-12T06:50:03.347636Z","shell.execute_reply":"2022-09-12T06:50:04.896716Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from IPython.display import Audio, IFrame, display\n\ndisplay(Audio(y1,rate=sr))\ndisplay(Audio(y2,rate=sr))","metadata":{"execution":{"iopub.status.busy":"2022-09-12T06:50:04.899552Z","iopub.execute_input":"2022-09-12T06:50:04.899928Z","iopub.status.idle":"2022-09-12T06:50:04.931085Z","shell.execute_reply.started":"2022-09-12T06:50:04.899890Z","shell.execute_reply":"2022-09-12T06:50:04.930104Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"librosa.display.waveplot(y1,sr=sr, x_axis='time');","metadata":{"execution":{"iopub.status.busy":"2022-09-12T06:50:04.932382Z","iopub.execute_input":"2022-09-12T06:50:04.932894Z","iopub.status.idle":"2022-09-12T06:50:05.194712Z","shell.execute_reply.started":"2022-09-12T06:50:04.932857Z","shell.execute_reply":"2022-09-12T06:50:05.193813Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"librosa.display.waveplot(y2,sr=sr, x_axis='time');","metadata":{"execution":{"iopub.status.busy":"2022-09-12T06:50:05.198143Z","iopub.execute_input":"2022-09-12T06:50:05.198751Z","iopub.status.idle":"2022-09-12T06:50:05.386162Z","shell.execute_reply.started":"2022-09-12T06:50:05.198708Z","shell.execute_reply":"2022-09-12T06:50:05.385158Z"},"trusted":true},"execution_count":null,"outputs":[]},{"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');","metadata":{"execution":{"iopub.status.busy":"2022-09-12T06:50:05.389128Z","iopub.execute_input":"2022-09-12T06:50:05.389775Z","iopub.status.idle":"2022-09-12T06:50:05.566861Z","shell.execute_reply.started":"2022-09-12T06:50:05.389732Z","shell.execute_reply":"2022-09-12T06:50:05.566026Z"},"trusted":true},"execution_count":null,"outputs":[]},{"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');","metadata":{"execution":{"iopub.status.busy":"2022-09-12T06:50:05.568331Z","iopub.execute_input":"2022-09-12T06:50:05.568942Z","iopub.status.idle":"2022-09-12T06:50:05.733728Z","shell.execute_reply.started":"2022-09-12T06:50:05.568899Z","shell.execute_reply":"2022-09-12T06:50:05.732843Z"},"trusted":true},"execution_count":null,"outputs":[]},{"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","metadata":{"execution":{"iopub.status.busy":"2022-09-12T06:50:05.735057Z","iopub.execute_input":"2022-09-12T06:50:05.735603Z","iopub.status.idle":"2022-09-12T06:50:05.743450Z","shell.execute_reply.started":"2022-09-12T06:50:05.735566Z","shell.execute_reply":"2022-09-12T06:50:05.742378Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"yf1 = f_high(y1, sr)\nyf2 = f_high(y2, sr)","metadata":{"execution":{"iopub.status.busy":"2022-09-12T06:50:05.744894Z","iopub.execute_input":"2022-09-12T06:50:05.745156Z","iopub.status.idle":"2022-09-12T06:50:05.760186Z","shell.execute_reply.started":"2022-09-12T06:50:05.745132Z","shell.execute_reply":"2022-09-12T06:50:05.759547Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Let's see...","metadata":{}},{"cell_type":"code","source":"librosa.display.waveplot(y1,sr=sr, x_axis='time');\nlibrosa.display.waveplot(yf1,sr=sr, x_axis='time');","metadata":{"execution":{"iopub.status.busy":"2022-09-12T06:50:05.761866Z","iopub.execute_input":"2022-09-12T06:50:05.762499Z","iopub.status.idle":"2022-09-12T06:50:06.004810Z","shell.execute_reply.started":"2022-09-12T06:50:05.762458Z","shell.execute_reply":"2022-09-12T06:50:06.003865Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"librosa.display.waveplot(y2,sr=sr, x_axis='time');\nlibrosa.display.waveplot(yf2,sr=sr, x_axis='time');","metadata":{"execution":{"iopub.status.busy":"2022-09-12T06:50:06.006290Z","iopub.execute_input":"2022-09-12T06:50:06.006878Z","iopub.status.idle":"2022-09-12T06:50:06.212430Z","shell.execute_reply.started":"2022-09-12T06:50:06.006837Z","shell.execute_reply":"2022-09-12T06:50:06.211473Z"},"trusted":true},"execution_count":null,"outputs":[]},{"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');","metadata":{"execution":{"iopub.status.busy":"2022-09-12T06:50:06.214261Z","iopub.execute_input":"2022-09-12T06:50:06.214656Z","iopub.status.idle":"2022-09-12T06:50:06.381467Z","shell.execute_reply.started":"2022-09-12T06:50:06.214616Z","shell.execute_reply":"2022-09-12T06:50:06.380506Z"},"trusted":true},"execution_count":null,"outputs":[]},{"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');","metadata":{"execution":{"iopub.status.busy":"2022-09-12T06:50:06.383223Z","iopub.execute_input":"2022-09-12T06:50:06.383574Z","iopub.status.idle":"2022-09-12T06:50:06.546053Z","shell.execute_reply.started":"2022-09-12T06:50:06.383536Z","shell.execute_reply":"2022-09-12T06:50:06.545144Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"display(Audio(yf1,rate=sr))\ndisplay(Audio(yf2,rate=sr))","metadata":{"execution":{"iopub.status.busy":"2022-09-12T06:50:06.547957Z","iopub.execute_input":"2022-09-12T06:50:06.548383Z","iopub.status.idle":"2022-09-12T06:50:06.576680Z","shell.execute_reply.started":"2022-09-12T06:50:06.548341Z","shell.execute_reply":"2022-09-12T06:50:06.575636Z"},"trusted":true},"execution_count":null,"outputs":[]},{"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)","metadata":{"execution":{"iopub.status.busy":"2022-09-12T06:50:06.578400Z","iopub.execute_input":"2022-09-12T06:50:06.578721Z","iopub.status.idle":"2022-09-12T06:50:06.599661Z","shell.execute_reply.started":"2022-09-12T06:50:06.578683Z","shell.execute_reply":"2022-09-12T06:50:06.598880Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"librosa.display.specshow(Dp1, x_axis='time', y_axis='mel');","metadata":{"execution":{"iopub.status.busy":"2022-09-12T06:50:06.601236Z","iopub.execute_input":"2022-09-12T06:50:06.601588Z","iopub.status.idle":"2022-09-12T06:50:06.754230Z","shell.execute_reply.started":"2022-09-12T06:50:06.601553Z","shell.execute_reply":"2022-09-12T06:50:06.753440Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"librosa.display.specshow(Dp2, x_axis='time', y_axis='mel');","metadata":{"execution":{"iopub.status.busy":"2022-09-12T06:50:06.755575Z","iopub.execute_input":"2022-09-12T06:50:06.756015Z","iopub.status.idle":"2022-09-12T06:50:06.948084Z","shell.execute_reply.started":"2022-09-12T06:50:06.755974Z","shell.execute_reply":"2022-09-12T06:50:06.947182Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"yp1 = librosa.feature.inverse.mel_to_audio(Dp1)\nyp2 = librosa.feature.inverse.mel_to_audio(Dp2)","metadata":{"execution":{"iopub.status.busy":"2022-09-12T06:50:06.949689Z","iopub.execute_input":"2022-09-12T06:50:06.950367Z","iopub.status.idle":"2022-09-12T06:50:11.403671Z","shell.execute_reply.started":"2022-09-12T06:50:06.950314Z","shell.execute_reply":"2022-09-12T06:50:11.402751Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"librosa.display.waveplot(yp1,sr=sr, x_axis='time');\nlibrosa.display.waveplot(y1,sr=sr, x_axis='time');","metadata":{"execution":{"iopub.status.busy":"2022-09-12T06:50:11.405148Z","iopub.execute_input":"2022-09-12T06:50:11.405524Z","iopub.status.idle":"2022-09-12T06:50:11.601485Z","shell.execute_reply.started":"2022-09-12T06:50:11.405485Z","shell.execute_reply":"2022-09-12T06:50:11.600542Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"librosa.display.waveplot(yp2,sr=sr, x_axis='time');\nlibrosa.display.waveplot(y2,sr=sr, x_axis='time');","metadata":{"execution":{"iopub.status.busy":"2022-09-12T06:50:11.602918Z","iopub.execute_input":"2022-09-12T06:50:11.603578Z","iopub.status.idle":"2022-09-12T06:50:11.788951Z","shell.execute_reply.started":"2022-09-12T06:50:11.603532Z","shell.execute_reply":"2022-09-12T06:50:11.787792Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"display(Audio(yp1,rate=sr))\ndisplay(Audio(yp2,rate=sr))","metadata":{"execution":{"iopub.status.busy":"2022-09-12T06:50:11.790737Z","iopub.execute_input":"2022-09-12T06:50:11.791141Z","iopub.status.idle":"2022-09-12T06:50:11.817906Z","shell.execute_reply.started":"2022-09-12T06:50:11.791101Z","shell.execute_reply":"2022-09-12T06:50:11.816957Z"},"trusted":true},"execution_count":null,"outputs":[]},{"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    if verbose:\n        start = time.time()\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    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    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    mask_gain_dB = np.min(_amp_to_db(np.abs(sig_stft)))\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    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    sig_mask = sig_stft_db < db_thresh\n    if verbose:\n        print(\"Masking:\", td(seconds=time.time() - start))\n        start = time.time()\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    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","metadata":{"execution":{"iopub.status.busy":"2022-09-12T06:50:11.819236Z","iopub.execute_input":"2022-09-12T06:50:11.819753Z","iopub.status.idle":"2022-09-12T06:50:11.854564Z","shell.execute_reply.started":"2022-09-12T06:50:11.819717Z","shell.execute_reply":"2022-09-12T06:50:11.853286Z"},"trusted":true},"execution_count":null,"outputs":[]},{"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)","metadata":{"execution":{"iopub.status.busy":"2022-09-12T06:50:11.856136Z","iopub.execute_input":"2022-09-12T06:50:11.856579Z","iopub.status.idle":"2022-09-12T06:50:12.086779Z","shell.execute_reply.started":"2022-09-12T06:50:11.856543Z","shell.execute_reply":"2022-09-12T06:50:12.085949Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"librosa.display.waveplot(y1,sr=sr, x_axis='time');\nlibrosa.display.waveplot(yg1,sr=sr, x_axis='time');","metadata":{"execution":{"iopub.status.busy":"2022-09-12T06:50:12.088233Z","iopub.execute_input":"2022-09-12T06:50:12.088684Z","iopub.status.idle":"2022-09-12T06:50:12.526706Z","shell.execute_reply.started":"2022-09-12T06:50:12.088645Z","shell.execute_reply":"2022-09-12T06:50:12.526020Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"librosa.display.waveplot(y2,sr=sr, x_axis='time');\nlibrosa.display.waveplot(yg2,sr=sr, x_axis='time');","metadata":{"execution":{"iopub.status.busy":"2022-09-12T06:50:12.528078Z","iopub.execute_input":"2022-09-12T06:50:12.528429Z","iopub.status.idle":"2022-09-12T06:50:12.712954Z","shell.execute_reply.started":"2022-09-12T06:50:12.528399Z","shell.execute_reply":"2022-09-12T06:50:12.712020Z"},"trusted":true},"execution_count":null,"outputs":[]},{"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');","metadata":{"execution":{"iopub.status.busy":"2022-09-12T06:50:12.714731Z","iopub.execute_input":"2022-09-12T06:50:12.715236Z","iopub.status.idle":"2022-09-12T06:50:12.880990Z","shell.execute_reply.started":"2022-09-12T06:50:12.715198Z","shell.execute_reply":"2022-09-12T06:50:12.880007Z"},"trusted":true},"execution_count":null,"outputs":[]},{"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');","metadata":{"execution":{"iopub.status.busy":"2022-09-12T06:50:12.883808Z","iopub.execute_input":"2022-09-12T06:50:12.884446Z","iopub.status.idle":"2022-09-12T06:50:13.052124Z","shell.execute_reply.started":"2022-09-12T06:50:12.884398Z","shell.execute_reply":"2022-09-12T06:50:13.051305Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"display(Audio(yg1,rate=sr))\ndisplay(Audio(yg2,rate=sr))","metadata":{"execution":{"iopub.status.busy":"2022-09-12T06:50:13.053502Z","iopub.execute_input":"2022-09-12T06:50:13.054031Z","iopub.status.idle":"2022-09-12T06:50:13.081142Z","shell.execute_reply.started":"2022-09-12T06:50:13.053989Z","shell.execute_reply":"2022-09-12T06:50:13.080417Z"},"trusted":true},"execution_count":null,"outputs":[]}]}