notebooks/1_signal_proc_pipeline.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "77341481",
   "metadata": {},
   "source": [
    "# Signal processing pipeline\n",
    "\n",
    "This notebook processes the audio signal in the following steps:\n",
    "1. Loads the 10-second audio signal at a sampling rate of 16k and converts it to mono using `librosa`.\n",
    "2. Normalizes the root-mean-square of the entire signal to 0.1.\n",
    "3. Uses `noisereduce` to segregate foreground (FG) and background (BG) audio signals. The following steps apply to each of the original (RAW), FG and BG signal.\n",
    "4. Uses `scikit-maad`, a package for quantifying environmental audio recordings, to extract the indices in both time and frequency domains (e.g., Bioacoustics Index, Acoustic Complexity Index) for each signal.\n",
    "5. Uses `soundsig` to calculate the spectrotemporal modulation power spectrum, which is essentially a 2D Fourier transformation of the audio spectrogram and is relevant to animal calling and human auditory perception.\n",
    "\n",
    "Note that extracting the embeddings of the pre-trained CNN models will be done in a separate notebook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "e6fb8e1f",
   "metadata": {},
   "outputs": [],
   "source": [
    "from __future__ import division\n",
    "def feature_eng(y, fs, vis=False):\n",
    "    \n",
    "    import librosa\n",
    "    import librosa.display\n",
    "    import numpy as np\n",
    "    import maad\n",
    "    from soundsig.sound import BioSound\n",
    "    import matplotlib.pyplot as plt\n",
    "    import pandas as pd\n",
    "    import warnings\n",
    "    warnings.simplefilter('ignore')\n",
    "    \n",
    "    # spectrum (not used)\n",
    "    ps = np.abs(np.fft.fft(y))**2\n",
    "    time_step = 1/fs\n",
    "    freqs = np.fft.fftfreq(y.size, time_step)\n",
    "    ps = ps[0:int((len(ps)/2)-1)] # take out the negative freq\n",
    "    freqs = freqs[0:int((len(freqs)/2)-1)]\n",
    "\n",
    "    # melspectrogram (not used)\n",
    "    S = librosa.feature.melspectrogram(y=y, sr=fs)\n",
    "    S_dB = librosa.power_to_db(S, ref=np.max)\n",
    "    \n",
    "    # alpha indices\n",
    "    Sxx_power,tn,fn,ext = maad.sound.spectrogram (y, fs, mode='psd')\n",
    "    df_temporal_indices = maad.features.all_temporal_alpha_indices(y, fs)\n",
    "    df_spectral_indices, _ = maad.features.all_spectral_alpha_indices(Sxx_power,tn,fn, extent=ext)\n",
    "    df_indices = pd.concat([df_temporal_indices,df_spectral_indices], axis=1)\n",
    "    \n",
    "    # soundsig for modulation power spectrum\n",
    "    \n",
    "    myBioSound = BioSound(soundWave=y, fs=fs)\n",
    "    myBioSound.mpsCalc(window=1, Norm = True)\n",
    "    \n",
    "    # reduce the dimension of MPS to one quadrant\n",
    "    len1 = int(len(myBioSound.wf-1)/2)\n",
    "    len2 = int(len(myBioSound.wt-1)/2)\n",
    "    quad1 = myBioSound.mps[len1:,len2:]\n",
    "    quad2 = np.fliplr(myBioSound.mps[len1:,:len2+1])\n",
    "    mps = (quad1+quad2)/2\n",
    "    wf = myBioSound.wf[len1:]\n",
    "    wt = myBioSound.wt[len2:]\n",
    "    \n",
    "    if vis==True:\n",
    "        plt.figure()\n",
    "        plt.plot(np.arange(0,y.size/fs,1/fs), y)\n",
    "        plt.show()\n",
    "        \n",
    "        fig, ax = plt.subplots()\n",
    "        img = librosa.display.specshow(S_dB, x_axis='time',\n",
    "                                 y_axis='mel', sr=fs,\n",
    "                                 fmax=fs/2, ax=ax)\n",
    "        fig.colorbar(img, ax=ax, format='%+2.0f dB')\n",
    "        ax.set(title='Mel-frequency spectrogram')\n",
    "        \n",
    "        DBNOISE=100\n",
    "        plt.figure()\n",
    "        plt.clf()\n",
    "        cmap = plt.get_cmap('jet')\n",
    "        ex = (myBioSound.wt.min(), myBioSound.wt.max(), myBioSound.wf.min()*1e3, myBioSound.wf.max()*1e3)\n",
    "        logMPS = 10.0*np.log10(myBioSound.mps)\n",
    "        maxMPS = logMPS.max()\n",
    "        minMPS = maxMPS-DBNOISE\n",
    "        logMPS[logMPS < minMPS] = minMPS\n",
    "        plt.imshow(logMPS, interpolation='nearest', aspect='auto', origin='lower', cmap=cmap, extent=ex)\n",
    "        plt.ylabel('Spectral Frequency (Cycles/KHz)')\n",
    "        plt.xlabel('Temporal Frequency (Hz)')\n",
    "        plt.colorbar()\n",
    "        plt.ylim((0,myBioSound.wf.max()*1e3))\n",
    "        plt.title('Modulation Power Spectrum')\n",
    "        plt.show()\n",
    "        \n",
    "    \n",
    "    return ps, freqs, S_dB, df_indices, mps, wt, wf\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "aa8a82a4",
   "metadata": {},
   "outputs": [],
   "source": [
    "def preproc(file, fs=16000, vis=False):\n",
    "    import librosa\n",
    "    import numpy as np\n",
    "    import matplotlib.pyplot as plt\n",
    "    import noisereduce as nr\n",
    "    \n",
    "    raw_y, fs = librosa.load(file, sr=fs, duration=10, mono = True)\n",
    "    y_mono_rs = raw_y - np.mean(raw_y) # remove DC\n",
    "    rms = np.sqrt(np.mean(y_mono_rs**2)) # get rms\n",
    "    y = y_mono_rs/(rms/0.1) # normalize the rms to 0.1\n",
    "    \n",
    "    fg_y = nr.reduce_noise(y=y, sr=fs)\n",
    "    bg_y = y - fg_y\n",
    "    \n",
    "    \n",
    "    if vis == True: print('++++++++++++++++++++++++ raw ++++++++++++++++++++++++')\n",
    "    ps, freqs, S_dB, df_indices, mps, wt, wf = feature_eng(y, fs, vis)\n",
    "    if vis == True: print('++++++++++++++++++++ foreground ++++++++++++++++++++')\n",
    "    ps_fg, freqs_fg, S_dB_fg, df_indices_fg, mps_fg, wt_fg, wf_fg = feature_eng(fg_y, fs, vis)\n",
    "    if vis == True: print('++++++++++++++++++++ background ++++++++++++++++++++')\n",
    "    ps_bg, freqs_fg, S_dB_bg, df_indices_bg, mps_bg, wt_bg, wf_bg = feature_eng(bg_y, fs, vis)\n",
    "    \n",
    "    output = {'fs': fs, 'y': y, 'fg_y': fg_y, 'bg_y': bg_y,\n",
    "             'df_indices': df_indices, 'mps': mps, 'wt': wt, 'wf': wf, \n",
    "             'df_indices_fg': df_indices_fg, 'mps_fg': mps_fg,\n",
    "             'df_indices_bg': df_indices_bg, 'mps_bg': mps_bg}\n",
    "    \n",
    "    return output\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "1d6c1bfb",
   "metadata": {},
   "outputs": [],
   "source": [
    "# scan data directories\n",
    "\n",
    "import glob\n",
    "import pandas as pd\n",
    "nature_file_list = []\n",
    "nature_file_list += glob.glob('../data/raw/GoogleAudioSet_unbalanced_list/Outside, rural or natural/*')\n",
    "nature_file_list += glob.glob('../data/raw/GoogleAudioSet_eval_list/Outside, rural or natural/*')\n",
    "\n",
    "city_file_list = []\n",
    "city_file_list += glob.glob('../data/raw/GoogleAudioSet_unbalanced_list/Outside, urban or manmade/*')\n",
    "city_file_list += glob.glob('../data/raw/GoogleAudioSet_eval_list/Outside, urban or manmade/*')\n",
    "\n",
    "\n",
    "remove_list_nature = glob.glob('../data/raw/GoogleAudioSet_unbalanced_list/Outside, rural or natural/AE0v7LesLZo*')\n",
    "remove_list_city = glob.glob('../data/raw/GoogleAudioSet_unbalanced_list/Outside, urban or manmade/OPyovt30GPQ*')\n",
    "remove_list_city += glob.glob('../data/raw/GoogleAudioSet_unbalanced_list/Outside, urban or manmade/CBmGYSOoeto*')\n",
    "\n",
    "\n",
    "nature_file_list = [ele for ele in nature_file_list if ele not in remove_list_nature]\n",
    "city_file_list = [ele for ele in city_file_list if ele not in remove_list_city]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "14f4b8a3",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "import time\n",
    "import pickle\n",
    "import os\n",
    "\n",
    "from joblib import Parallel, delayed\n",
    "\n",
    "def run_preproc(file_name):\n",
    "    save_file_name = '../data/interim/'+file_name[12:-4]+'.pkl'\n",
    "    if not os.path.isfile(save_file_name): # run the script only if the file does not exist\n",
    "        if not os.path.exists(save_file_name.rsplit('/', 1)[0]): # create the folder if the folder does not exist\n",
    "            os.makedirs(save_file_name.rsplit('/', 1)[0]) # extract the folder of the file path\n",
    "            \n",
    "\n",
    "        start_time = time.time()\n",
    "\n",
    "        output = preproc(file_name, vis = False)\n",
    "        f = open(save_file_name,'wb') # create a binary pickle file \n",
    "        pickle.dump(output,f)\n",
    "        f.close()\n",
    "\n",
    "        print(\"--- %s seconds ---\" % (time.time() - start_time))\n",
    "        \n",
    "\n",
    "# run the process in parallel\n",
    "Parallel(n_jobs=-1)(delayed(run_preproc)(file_name) for file_name in nature_file_list+city_file_list)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5b7844be",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# for file_name in nature_file_list+city_file_list:\n",
    "#     print(file_name)\n",
    "#     run_preproc(file_name)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.16"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}