{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "18d3113c-f210-4061-9a72-b79052c71a2e",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "--- PROCESSING: 20230504_081356UTC_V12 ---\n",
      "Detected total duration: 96.54 seconds\n",
      "Applying high-pass filter...\n",
      "\n",
      " -> Running PERMISSIVE strategy (Thresholds: [40 40 30])\n",
      "    Detected 1158 clicks.\n",
      "    Saved detection parameters: ../01_data/20230504_081356UTC_V12_detections_permissive.csv\n",
      "\n",
      " -> Running STRICT strategy (Thresholds: [80 80 55])\n",
      "    Detected 834 clicks.\n",
      "    Saved detection parameters: ../01_data/20230504_081356UTC_V12_detections_strict.csv\n",
      "\n",
      "--- PROCESSING: 20230504_085315UTC_V12 ---\n",
      "Detected total duration: 96.54 seconds\n",
      "Applying high-pass filter...\n",
      "\n",
      " -> Running PERMISSIVE strategy (Thresholds: [40 40 30])\n",
      "    Detected 1445 clicks.\n",
      "    Saved detection parameters: ../01_data/20230504_085315UTC_V12_detections_permissive.csv\n",
      "\n",
      " -> Running STRICT strategy (Thresholds: [80 80 55])\n",
      "    Detected 1305 clicks.\n",
      "    Saved detection parameters: ../01_data/20230504_085315UTC_V12_detections_strict.csv\n",
      "\n",
      "--- PROCESSING: 20230504_085452UTC_V12 ---\n",
      "Detected total duration: 96.54 seconds\n",
      "Applying high-pass filter...\n",
      "\n",
      " -> Running PERMISSIVE strategy (Thresholds: [40 40 30])\n",
      "    Detected 558 clicks.\n",
      "    Saved detection parameters: ../01_data/20230504_085452UTC_V12_detections_permissive.csv\n",
      "\n",
      " -> Running STRICT strategy (Thresholds: [80 80 55])\n",
      "    Detected 466 clicks.\n",
      "    Saved detection parameters: ../01_data/20230504_085452UTC_V12_detections_strict.csv\n"
     ]
    }
   ],
   "source": [
    "\"\"\"\n",
    "\n",
    "Module 01: Click detection and TDOA extraction\n",
    "Journal of the Acoustical Society of America (JASA) - Data availability\n",
    "\n",
    "(c) Walter M.X. Zimmer (2025)\n",
    "Modified by Lara Berkenbaum and Hervé Glotin (2025, 2026)\n",
    "Original code produced for the CIAN Tool Kit, associated with the manuscript: \n",
    "Berkenbaum, L., Glotin, H., Sarano, F., Zimmer, W.M.X., Sarano, V., Adam, O., Giraudet, P. (2026). \"Audiovisual diarization of overlapping click trains in sperm whale (Physeter macrocephalus) vocal sparring using a three-hydrophone array,\" J. Acoust. Soc. Am.\n",
    "\n",
    "Description:\n",
    "This script processes the raw 5-channel acoustic recordings to detect sperm whale click transients \n",
    "using a Teager-Kaiser Energy Operator (TKEO). \n",
    "It extracts the Time Difference of Arrival (TDOA) cross-correlation delays across the 3-hydrophone \n",
    "sub-array (channels 0, 1, and 4) using a two-stage asymmetric detection strategy.\n",
    "\n",
    "Inputs: \n",
    "- Raw audio files located in `../01_data/raw_audio/*.wav`\n",
    "\n",
    "Outputs: \n",
    "- `[filename]_detections_[strategy].csv`: per-click arrival times (T_H0, T_H1, T_H2), cross-correlation delays in samples \n",
    "(B01_samples, B02_samples, B12_samples), energy (energy_dB), loop-closure error (loop_error_cm), and the sampling rate (fs_hz).\n",
    "\n",
    "Two files are produced per recording: `_detections_permissive.csv` (for the clustering pipeline, Modules 02 and 03) \n",
    "and `_detections_strict.csv` (for the spatial localization pipeline, Module 04).\n",
    "\n",
    "\"\"\"\n",
    "\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "from scipy import signal\n",
    "import soundfile as sf\n",
    "import os\n",
    "import glob\n",
    "\n",
    "# Import custom mathematical functions and shared constants from the local utils module\n",
    "from utils import wav_read, leak, extract, fftCorr, quadInt, C_SOUND, HSEL\n",
    "\n",
    "# =================================================================\n",
    "# 1. CONFIGURATION & AUTOMATION\n",
    "# =================================================================\n",
    "# Target directory containing the raw .wav files\n",
    "audio_dir = \"../01_data/raw_audio/\"\n",
    "audio_files = glob.glob(os.path.join(audio_dir, \"*.wav\"))\n",
    "\n",
    "if not audio_files:\n",
    "    print(f\"No audio files found in {audio_dir}. Please check the directory structure.\")\n",
    "\n",
    "# Two-stage asymmetric detection strategy \n",
    "# Thresholds defined as: [TH_H0, TH_H1, TH_H2]\n",
    "threshold_strategies = {\n",
    "    \"permissive\": np.array([40, 40, 30]), # Maximizes track continuity for clustering\n",
    "    \"strict\": np.array([80, 80, 55])      # Strict CFAR for spatial TDOA localization\n",
    "}\n",
    "\n",
    "# =================================================================\n",
    "# 2. BATCH PROCESSING PIPELINE\n",
    "# =================================================================\n",
    "for fname in audio_files:\n",
    "    base_name = os.path.splitext(os.path.basename(fname))[0]\n",
    "    \n",
    "    audio_info = sf.info(fname)\n",
    "    total_duration = audio_info.duration\n",
    "    \n",
    "    to = 0.0\n",
    "    duration = total_duration\n",
    "    \n",
    "    print(f\"\\n--- PROCESSING: {base_name} ---\")\n",
    "    print(f\"Detected total duration: {duration:.2f} seconds\")\n",
    "    \n",
    "    # Load and pre-process audio\n",
    "    data, fs, nd = wav_read(fname, to, duration)\n",
    "    fsk = fs / 1000.0\n",
    "    \n",
    "    # Apply 5kHz high-pass filter to mitigate oceanic background noise\n",
    "    print(\"Applying high-pass filter...\")\n",
    "    sos = signal.butter(4, 5000.0, btype='highpass', fs=fs, output='sos')\n",
    "    data = signal.sosfiltfilt(sos, data, axis=0)\n",
    "    \n",
    "    # Select hardware channels [0, 1, 4] corresponding to the 3-hydrophone array\n",
    "    xx = data[:, [0, 1, 4]]\n",
    "    \n",
    "    # Apply empirical hardware correction and polarity inversion to the C75 channel\n",
    "    xx *= np.array([1, 1, -3]) \n",
    "    \n",
    "    # Apply TKEO (Teager-Kaiser Energy Operator) for transient enhancement\n",
    "    uu = np.sqrt(np.maximum(0, xx[1:-1, :]**2 - xx[2:, :] * xx[:-2, :]))\n",
    "    med_uu = np.median(uu, 0)\n",
    "    uu /= med_uu \n",
    "    \n",
    "    # Execute detection for both asymmetric strategies\n",
    "    for strategy_name, TH in threshold_strategies.items():\n",
    "        print(f\"\\n -> Running {strategy_name.upper()} strategy (Thresholds: {TH})\")\n",
    "        \n",
    "        FILE_OUT_PARAMS = f\"../01_data/{base_name}_detections_{strategy_name}.csv\"\n",
    "        \n",
    "        yy = uu > TH\n",
    "        aa = 0.01\n",
    "        cc = fsk\n",
    "        zz = (uu > TH) * cc\n",
    "        zz = leak(zz, cc, aa) \n",
    "        zz = zz > aa\n",
    "        \n",
    "        vv = np.diff(np.int32(zz), axis=0)\n",
    "        idet1 = np.where(vv > 0)\n",
    "        \n",
    "        # Isolate initial detection frames\n",
    "        i1 = idet1[1] == 0 \n",
    "        icl = idet1[0][i1]\n",
    "        \n",
    "        num_clicks = len(icl[1:])\n",
    "        print(f\"    Detected {num_clicks} clicks.\")\n",
    "        \n",
    "        ns = int(8 * fsk) # 8 ms extraction window\n",
    "        rows_params = []\n",
    "        \n",
    "        # Extract features and compute TDOA delays for each click\n",
    "        for idx, ii in enumerate(icl[1:]):\n",
    "            \n",
    "            # Signal extraction around the transient peak\n",
    "            u2 = extract(xx, ii, -500, ns) \n",
    "            \n",
    "            # Sub-sample time delay calculation via frequency-domain cross-correlation\n",
    "            v2_full = fftCorr(u2, HSEL, ns)\n",
    "            v2 = v2_full[ns//2-100 : ns//2+100, :]\n",
    "            \n",
    "            idl1 = np.argmax(v2, 0)\n",
    "            jdl2, adl = quadInt(v2, idl1)\n",
    "            jdl2 -= 100 \n",
    "            \n",
    "            # Absolute time calculation (applying standard delay convention)\n",
    "            t_h0_abs = to + (ii / fs) \n",
    "            dt_sec = jdl2 / fs\n",
    "            t_h1_abs = t_h0_abs + dt_sec[0]\n",
    "            t_h2_abs = t_h0_abs + dt_sec[1]\n",
    "            \n",
    "            # Energy metric (dB)\n",
    "            energy_linear = np.mean(u2**2)\n",
    "            energy_dB = 10 * np.log10(energy_linear + 1e-12)\n",
    "            \n",
    "            # Geometric loop error calculation\n",
    "            loop_samples = jdl2[0] - jdl2[1] + jdl2[2]\n",
    "            loop_error_cm = np.abs((loop_samples / fs) * C_SOUND * 100)\n",
    "            \n",
    "            # Append validated parameters\n",
    "            row_data = {\n",
    "                'ID_click': idx + 1,\n",
    "                'index_frame': ii,\n",
    "                'T_H0': t_h0_abs,\n",
    "                'T_H1': t_h1_abs,\n",
    "                'T_H2': t_h2_abs,\n",
    "                'B01_samples': jdl2[0],\n",
    "                'B02_samples': jdl2[1],\n",
    "                'B12_samples': jdl2[2],\n",
    "                'energy_dB': energy_dB,\n",
    "                'loop_error_cm': loop_error_cm,\n",
    "                'fs_hz': fs\n",
    "            }\n",
    "            rows_params.append(row_data)\n",
    "        \n",
    "        # Post-processing and CSV export\n",
    "        if len(rows_params) > 0:\n",
    "            df = pd.DataFrame(rows_params)\n",
    "        \n",
    "            cols_order = [\n",
    "                'ID_click', 'index_frame', 'T_H0', 'T_H1', 'T_H2', \n",
    "                'B01_samples', 'B02_samples', 'B12_samples',\n",
    "                'energy_dB', 'loop_error_cm', 'fs_hz'\n",
    "            ]\n",
    "            df = df[[c for c in cols_order if c in df.columns]]\n",
    "            df.to_csv(FILE_OUT_PARAMS, index=False, sep=';')\n",
    "            print(f\"    Saved detection parameters: {FILE_OUT_PARAMS}\")\n",
    "            \n",
    "        else:\n",
    "            print(\"    No clicks detected for this strategy.\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0fdea904-77ea-4693-867c-2e4e512301c0",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.11.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}