lostanlen_waspaa2019_fig1.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import librosa\n",
    "import numpy as np\n",
    "from librosa.display import specshow\n",
    "import matplotlib\n",
    "from matplotlib import pyplot as plt\n",
    "import glob\n",
    "import os\n",
    "import datetime\n",
    "import IPython.display\n",
    "import pandas as pd\n",
    "import collections\n",
    "work_dir = '/beegfs/vl1019/waspaa2019_data/coni_knight'\n",
    "paths = glob.glob(work_dir + \"/*.wav\")\n",
    "\n",
    "%matplotlib inline\n",
    "\n",
    "# Define PCEN settings.\n",
    "settings = {\n",
    "    \"fmin\": 2000,\n",
    "    \"fmax\": 11025,\n",
    "    \"hop_length\": 32,\n",
    "    \"n_fft\": 1024,\n",
    "    \"n_mels\": 128,\n",
    "    \"sr\": 22050,\n",
    "    \"win_length\": 256,\n",
    "    \"window\": \"flattop\",\n",
    "    \"T\": 1.0,\n",
    "    \"alpha\": 1.0,\n",
    "    \"delta\": 0.0,\n",
    "    \"r\": 1.0,\n",
    "    \"eps\": 1e-6}\n",
    "\n",
    "plt.rcParams[\"font.family\"] = \"serif\"\n",
    "\n",
    "def parse_coni_knight(path):\n",
    "    meta = {}\n",
    "    name = os.path.split(path)[1]\n",
    "    meta[\"file_name\"] = name\n",
    "    meta[\"distance\"] = int(name.split(\"_\")[0])\n",
    "    meta['vocalization_id'] = int(name.split(\"_\")[1])\n",
    "    meta['nominal_distance'] = int(name.split(\"_\")[2].split(\"-\")[2])\n",
    "    meta['datetime'] = datetime.datetime(int(name.split(\"_\")[3][:4]), int(name.split(\"_\")[3][4:6]), int(name.split(\"_\")[3][7:]), int(name.split(\"_\")[4][:2]), int(name.split(\"_\")[4][2:4]))\n",
    "    return meta\n",
    "\n",
    "data = map(parse_coni_knight, paths)\n",
    "df = pd.DataFrame(list(data))\n",
    "\n",
    "voc_ids_counter = collections.Counter(df[\"vocalization_id\"])\n",
    "voc_ids = sorted(voc_ids_counter.keys())\n",
    "voc_id = voc_ids[-3]\n",
    "\n",
    "names_dict = {}\n",
    "for voc_id in voc_ids:\n",
    "    filtered_df = df[df[\"vocalization_id\"] == voc_id]\n",
    "    sorted_df = filtered_df.sort_values(by=['nominal_distance'])\n",
    "    names = list(sorted_df[\"file_name\"])\n",
    "    names_dict[voc_id] = names\n",
    "\n",
    "\n",
    "# Define PCEN settings.\n",
    "settings = {\n",
    "    \"fmin\": 2000,\n",
    "    \"fmax\": 11025,\n",
    "    \"hop_length\": 32,\n",
    "    \"n_fft\": 1024,\n",
    "    \"n_mels\": 128,\n",
    "    \"sr\": 22050,\n",
    "    \"win_length\": 256,\n",
    "    \"window\": \"flattop\",\n",
    "    \"T\": 1.0,\n",
    "    \"alpha\": 1.0,\n",
    "    \"delta\": 0.0,\n",
    "    \"r\": 1.0,\n",
    "    \"eps\": 1e-6}\n",
    "\n",
    "\n",
    "voc_id = voc_ids[-5]\n",
    "#vod_id = 3634\n",
    "\n",
    "wav_paths = []\n",
    "for path in names_dict[voc_id]:\n",
    "    wav_path = work_dir + os.path.expanduser(\"/\") + path\n",
    "    wav_paths.append(wav_path)\n",
    "    \n",
    "\n",
    "subset = range(11)\n",
    "#subset = [0, 1, 3, 7, 9, 10]\n",
    "wav_paths = [wav_paths[i] for i in subset]\n",
    "\n",
    "Es, PCENs = [], []\n",
    "lms_ranges, pcen_ranges = [], []\n",
    "\n",
    "\n",
    "for wav_path in wav_paths:\n",
    "    # Load, resample, rescale, and pad waveform.\n",
    "    waveform, sample_rate = librosa.load(wav_path, sr=None)\n",
    "    waveform = librosa.resample(\n",
    "        waveform, sample_rate, settings[\"sr\"])\n",
    "    \n",
    "    waveform = waveform * (2**31)\n",
    "    waveform = np.concatenate((waveform, waveform, waveform[::-1]))\n",
    "\n",
    "    # Compute short-term Fourier transform (STFT)\n",
    "    stft = librosa.stft(\n",
    "        waveform,\n",
    "        n_fft=settings[\"n_fft\"],\n",
    "        win_length=settings[\"win_length\"],\n",
    "        hop_length=settings[\"hop_length\"],\n",
    "        window=settings[\"window\"])\n",
    "    \n",
    "    # Compute the squared complex modulus (\"abs2\") of the STFT.\n",
    "    abs2_stft = (stft.real*stft.real) + (stft.imag*stft.imag)\n",
    "\n",
    "    # Transform the STFT into a mel-frequency spectrogram.\n",
    "    E = librosa.feature.melspectrogram(\n",
    "        y=None,\n",
    "        S=abs2_stft,\n",
    "        sr=settings[\"sr\"],\n",
    "        n_fft=settings[\"n_fft\"],\n",
    "        n_mels=settings[\"n_mels\"],\n",
    "        htk=True,\n",
    "        fmin=settings[\"fmin\"],\n",
    "        fmax=settings[\"fmax\"])\n",
    "\n",
    "    # Apply per-channel energy normalization (PCEN).\n",
    "    PCEN = librosa.pcen(\n",
    "        S=E,\n",
    "        time_constant=settings[\"T\"],\n",
    "        sr=settings[\"sr\"],\n",
    "        hop_length=settings[\"hop_length\"],\n",
    "        gain=settings[\"alpha\"],\n",
    "        power=settings[\"r\"],\n",
    "        bias=settings[\"delta\"],\n",
    "        eps=settings[\"eps\"])\n",
    "    \n",
    "    # Unpad and append to list.\n",
    "    E = E[:,(E.shape[1]//3):]\n",
    "    PCEN = PCEN[:,(PCEN.shape[1]//3):]\n",
    "    Es.append(E)\n",
    "    PCENs.append(PCEN)\n",
    "    \n",
    "    spectral_flux = np.sum(np.maximum(0, np.diff(np.log(E[:, 10:(E.shape[1]//2 - 1)]), axis=1)), axis=0)\n",
    "    lms_range = np.max(spectral_flux) - np.min(spectral_flux)\n",
    "    lms_ranges.append(lms_range)\n",
    "    \n",
    "    pcen_range = np.max(PCEN[:, 10:(PCEN.shape[1]//2 - 1)]) - np.min(PCEN[:, 10:(PCEN.shape[1]//2 - 1)])\n",
    "    pcen_ranges.append(pcen_range)\n",
    "\n",
    "\n",
    "plt.figure(figsize=(5, 15))\n",
    "\n",
    "start_col = 5\n",
    "max_length = max([(E.shape[1]//2) for E in Es])\n",
    "logE_tensor = np.log1p(np.stack([E[:, start_col:max_length] for E in Es]))\n",
    "logE_tensor = (logE_tensor - np.min(logE_tensor)) / (np.max(logE_tensor) - np.min(logE_tensor))\n",
    "\n",
    "log1pPCEN_tensor = np.log1p(np.stack([PCEN[:, start_col:max_length] for PCEN in PCENs]))\n",
    "log1pPCEN_tensor = (log1pPCEN_tensor - np.min(log1pPCEN_tensor)) / (np.max(log1pPCEN_tensor) - np.min(log1pPCEN_tensor))\n",
    "cmap = \"magma\"\n",
    "\n",
    "# Like in BirdVoxDetect, we trim the top 8 frequency bands to avoid the effect of the Nyquist anti-aliasing filter.\n",
    "\n",
    "top_bin = 120 \n",
    "\n",
    "\n",
    "lower_sf_bin = 15\n",
    "upper_sf_bin = 90\n",
    "\n",
    "for i in range(len(wav_paths)):\n",
    "\n",
    "    # Load distance.\n",
    "    dist_str = str(int(names_dict[voc_id][subset[i]].split(\"_\")[2].split(\"-\")[2])) + \" m\" + \" \" * 11\n",
    "    \n",
    "    # Display logmelspec on left column.\n",
    "    plt.subplot(11, 2, 2*i+1)\n",
    "    specshow(logE_tensor[i, :top_bin, :], cmap=cmap)\n",
    "    argmax_row, argmax_col = np.unravel_index(\n",
    "        np.argmax(np.diff(logE_tensor[i, lower_sf_bin:upper_sf_bin, start_col:])),\n",
    "        logE_tensor[i, lower_sf_bin:upper_sf_bin, :].shape)\n",
    "    plt.scatter([start_col+argmax_col], [lower_sf_bin+argmax_row], s=80, facecolors='white')\n",
    "    plt.clim(0, 1)\n",
    "    h = plt.ylabel(dist_str, rotation=0)\n",
    "    if i == 0:\n",
    "        plt.title(\"Pointwise logarithm\")\n",
    "\n",
    "    # Display PCEN.\n",
    "    plt.subplot(11, 2, 2*i+2)\n",
    "    specshow(log1pPCEN_tensor[i, :top_bin, :]**(0.5), cmap=cmap)\n",
    "    argmax_row, argmax_col = np.unravel_index(\n",
    "        np.argmax(log1pPCEN_tensor[i, lower_sf_bin:upper_sf_bin, start_col:]),\n",
    "        log1pPCEN_tensor[i, lower_sf_bin:upper_sf_bin, start_col:].shape)\n",
    "    plt.scatter([start_col+argmax_col], [lower_sf_bin+argmax_row], s=80, facecolors='white')\n",
    "    plt.clim(0, 1)\n",
    "    if i == 0:\n",
    "        plt.title(\"PCEN\")\n",
    "        \n",
    "        \n",
    "plt.savefig(\"lostanlen_waspaa2019_coni-pcengrams.eps\", bbox_inches=\"tight\")\n",
    "plt.savefig(\"lostanlen_waspaa2019_coni-pcengrams.png\", bbox_inches=\"tight\", dpi=1000)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}