{ "cells": [ { "cell_type": "code", "execution_count": 1, "id": "0956fec9", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", "--- PROCESSING SESSION: 20230504_081356 ---\n", "Step 1: Loading audio (channel 0)\n", "Audio signal pre-processed successfully.\n", "Reading timestamps from: 20230504_081356_validated_clicks.csv\n", "Filtered out 280 excluded clicks. Retained 555 clicks for CQT extraction.\n", "Loaded 555 click events for processing.\n", "Step 2: CQT patch extraction and realignment\n", "Extraction complete. Valid patches: 555\n", "Step 3: saving CQT tensors and features\n", "Saved patches (Pickle): 20230504_081356_CQT_patch_H0.pkl\n", "Saved features (CSV): 20230504_081356_CQT_feat_H0.csv\n", "\n", "--- PROCESSING SESSION: 20230504_085452 ---\n", "Step 1: Loading audio (channel 0)\n", "Audio signal pre-processed successfully.\n", "Reading timestamps from: 20230504_085452_validated_clicks.csv\n", "Filtered out 102 excluded clicks. Retained 356 clicks for CQT extraction.\n", "Loaded 356 click events for processing.\n", "Step 2: CQT patch extraction and realignment\n", "Extraction complete. Valid patches: 356\n", "Step 3: saving CQT tensors and features\n", "Saved patches (Pickle): 20230504_085452_CQT_patch_H0.pkl\n", "Saved features (CSV): 20230504_085452_CQT_feat_H0.csv\n", "\n", "--- PROCESSING SESSION: 20230504_085315 ---\n", "Step 1: Loading audio (channel 0)\n", "Audio signal pre-processed successfully.\n", "Reading timestamps from: 20230504_085315_validated_clicks.csv\n", "Filtered out 950 excluded clicks. Retained 311 clicks for CQT extraction.\n", "Loaded 311 click events for processing.\n", "Step 2: CQT patch extraction and realignment\n", "Extraction complete. Valid patches: 311\n", "Step 3: saving CQT tensors and features\n", "Saved patches (Pickle): 20230504_085315_CQT_patch_H0.pkl\n", "Saved features (CSV): 20230504_085315_CQT_feat_H0.csv\n", "\n", "--- BATCH PROCESSING COMPLETE ---\n" ] } ], "source": [ "\"\"\"\n", "\n", "Module 02: Constant-Q Transform (CQT) extraction\n", "Journal of the Acoustical Society of America (JASA) - Data availability\n", "\n", "(c) Walter MX 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 computes a Constant-Q Transform (CQT) patch for each on validated sperm whale clicks \n", "retained for clustering (i.e., detected clicks without an excluded_reason). These time-frequency patches \n", "serve as the acoustic signatures used by the unsupervised click-train clustering algorithm (Module 03).\n", "A set of scalar spectral descriptors (centroid, bandwidth, entropy, energy) is also computed and \n", "stored for reference, but is not used by the clustering itself.\n", "\n", "Inputs: \n", "- Raw audio files: `../01_data/raw_audio/*.wav`\n", "- Manually validated clicks: `../01_data/ground_truth_acoustic/*_validated_clicks.csv`\n", "\n", "Outputs: \n", "- CQT Patches (Pickle): `../01_data/ground_truth_acoustic/[PREFIX]_CQT_patch_H0.pkl`\n", "- Spectral Features (CSV): `../01_data/ground_truth_acoustic/[PREFIX]_CQT_feat_H0.csv`\n", "\n", "\"\"\"\n", "\n", "import pandas as pd\n", "import numpy as np\n", "import soundfile as sf\n", "import pickle\n", "import sys\n", "from scipy.stats import entropy\n", "from numba import jit\n", "import os\n", "import glob\n", "\n", "# ==============================================================================\n", "# 1. MATHEMATICAL TOOLS & CQT ENGINE\n", "# ==============================================================================\n", "\n", "@jit(nopython=True, cache=True)\n", "def bit_length(n):\n", " if n == 0: return 0\n", " bits = -32\n", " m = 0\n", " while n:\n", " m = n\n", " n >>= 32\n", " bits += 32\n", " while m:\n", " m >>= 1\n", " bits += 1\n", " return bits\n", "\n", "def olafilt(b, x, zi=None):\n", " L_I = b.shape[0]\n", " L_F = 2 << bit_length(L_I)\n", " L_S = L_F - L_I + 1\n", " L_sig = x.shape[0]\n", " offsets = range(0, L_sig, L_S)\n", " if np.iscomplexobj(b) or np.iscomplexobj(x):\n", " fft_func, ifft_func = np.fft.fft, np.fft.ifft\n", " res = np.zeros(L_sig+L_F, dtype=np.complex128)\n", " else:\n", " fft_func, ifft_func = np.fft.rfft, np.fft.irfft\n", " res = np.zeros(L_sig+L_F)\n", " FDir = fft_func(b, n=L_F)\n", " for n in offsets:\n", " res[n:n+L_F] += ifft_func(fft_func(x[n:n+L_S], n=L_F)*FDir)\n", " if zi is not None:\n", " res[:zi.shape[0]] += zi\n", " return res[:L_sig], res[L_sig:]\n", " else:\n", " return res[:L_sig]\n", "\n", "def CQT_custom(xx, fs, fmin, fmax, fpo):\n", " noctave = np.ceil(np.log2(fmax/fmin))\n", " nvoice = fpo\n", " nscale = int(noctave*nvoice)\n", " ytick = np.arange(nscale)/nvoice\n", " Q = 1/(2**(1/nvoice)-1)\n", " fk = fmin * 2**ytick\n", " N = np.ceil(fs/fk*Q)\n", " nx = xx.shape[0]\n", " ndat = nx\n", " S = np.zeros((nscale, ndat))\n", " W = np.zeros((nscale, int(max(N))), dtype=np.complex64)\n", " for kk in range(nscale):\n", " nw = int(N[kk])\n", " W[kk, :nw] = np.hanning(nw)/nw * np.exp(-2*np.pi*1j*Q*np.arange(nw)/nw)\n", " for kk in range(nscale):\n", " nw = int(N[kk])\n", " if nw <= 2*Q: continue\n", " win = W[kk, :nw]\n", " nwh = nw//2\n", " xx1 = np.append(xx, np.zeros(nwh))\n", " yy = olafilt(win, xx1)\n", " yy = yy[nwh:]\n", " S[kk, :] = np.abs(yy)**2\n", " ts = np.arange(ndat)/fs\n", " return S, fk, ts\n", "\n", "def extract_cqt_scalars(patch, freqs):\n", " \"\"\"Extracts spectral statistics from the CQT patch.\"\"\"\n", " if patch.size == 0: \n", " return {'cqt_centroid': 0, 'cqt_bandwidth': 0, 'cqt_entropy': 0, 'cqt_energy': -120}\n", "\n", " mean_spec = np.mean(patch, axis=1) \n", " total_energy = np.sum(mean_spec)\n", " \n", " cqt_energy = 10 * np.log10(total_energy + 1e-20)\n", " \n", " if total_energy > 0:\n", " pdf = mean_spec / total_energy\n", " centroid = np.sum(freqs * pdf)\n", " bandwidth = np.sqrt(np.sum(((freqs - centroid)**2) * pdf))\n", " spec_entropy = entropy(pdf + 1e-12)\n", " else:\n", " centroid, bandwidth, spec_entropy = 0, 0, 0\n", " \n", " return {\n", " 'cqt_centroid': centroid,\n", " 'cqt_bandwidth': bandwidth,\n", " 'cqt_entropy': spec_entropy,\n", " 'cqt_energy': cqt_energy\n", " }\n", "\n", "# ==============================================================================\n", "# 2. CONFIGURATION & AUTOMATION\n", "# ==============================================================================\n", "csv_dir = \"../01_data/ground_truth_acoustic/\"\n", "audio_dir = \"../01_data/raw_audio/\"\n", "\n", "# Ingest only the manually validated clicks for feature extraction\n", "detection_files = glob.glob(os.path.join(csv_dir, \"*_validated_clicks.csv\"))\n", "\n", "if not detection_files:\n", " print(f\"Error: No validated click files (*_validated_clicks.csv) found in {csv_dir}.\")\n", " sys.exit(0)\n", "\n", "# CQT Parameters\n", "TARGET_CHANNEL = 0\n", "F_MIN, F_MAX, F_PO = 1200, 80000, 12\n", "WIN_SIZE_CQT = 0.08\n", "PATCH_WIDTH_S = 0.009\n", "REALIGN_WIN_S = 0.004\n", "\n", "# ==============================================================================\n", "# 3. BATCH PROCESSING PIPELINE\n", "# ==============================================================================\n", "for csv_filename in detection_files:\n", " # Extract the full session prefix dynamically (e.g., \"20230504_081356\")\n", " base_name = os.path.basename(csv_filename).replace(\"_validated_clicks.csv\", \"\")\n", " \n", " # Locate corresponding audio file using the full dynamic prefix\n", " audio_candidates = glob.glob(os.path.join(audio_dir, f\"*{base_name}*.wav\"))\n", " if not audio_candidates:\n", " print(f\"Warning: No matching audio file found for session {base_name}. Skipping.\")\n", " continue\n", " \n", " filename_wav = audio_candidates[0]\n", " \n", " print(f\"\\n--- PROCESSING SESSION: {base_name} ---\")\n", " print(f\"Step 1: Loading audio (channel {TARGET_CHANNEL})\")\n", " \n", " # Load and preprocess audio (broadband differentiation)\n", " with sf.SoundFile(filename_wav) as f:\n", " fs = f.samplerate\n", " data = f.read(always_2d=True)\n", " \n", " signal = data[:, TARGET_CHANNEL].astype(np.float32)\n", " \n", " xx = np.diff(signal, n=1)\n", " xx -= np.mean(xx)\n", " max_val = np.max(np.abs(xx))\n", " if max_val > 0: xx /= max_val\n", " \n", " audio_signal = xx\n", " print(\"Audio signal pre-processed successfully.\")\n", "\n", " # Load manual validation timestamps using robust multi-separator reader\n", " print(f\"Reading timestamps from: {os.path.basename(csv_filename)}\")\n", " \n", " separators = [',', ';', r'\\s+', None]\n", " df_raw = None\n", " found = False\n", "\n", " for sep in separators:\n", " try:\n", " temp_df = pd.read_csv(csv_filename, sep=sep, engine='python')\n", " temp_df.columns = temp_df.columns.str.strip()\n", " if \"T_H0\" in temp_df.columns:\n", " df_raw = temp_df\n", " found = True\n", " break\n", " except Exception:\n", " continue\n", " \n", " if not found or df_raw is None:\n", " print(f\"Error: unable to parse '{os.path.basename(csv_filename)}' or locate target columns. Skipping.\")\n", " continue\n", " \n", " # ---------------------------------------------------------\n", " # FILTERING: excluded clicks that will not be clustered\n", " # ---------------------------------------------------------\n", " if 'excluded_reason' in df_raw.columns:\n", " initial_count = len(df_raw)\n", " df_raw = df_raw[df_raw['excluded_reason'].isna()].copy()\n", " print(f\"Filtered out {initial_count - len(df_raw)} excluded clicks. Retained {len(df_raw)} clicks for CQT extraction.\")\n", " \n", " col_target = \"T_H0\"\n", " \n", " df_times = df_raw[[col_target]].copy()\n", " df_times.columns = ['t_global']\n", " df_times['t_global'] = pd.to_numeric(df_times['t_global'], errors='coerce')\n", " \n", " if 'index' in df_raw.columns:\n", " df_times['index'] = df_raw['index']\n", " else:\n", " print(f\"Error: 'index' column missing from {os.path.basename(csv_filename)}. Crucial for matching. Skipping.\")\n", " continue\n", " \n", " print(f\"Loaded {len(df_times)} click events for processing.\")\n", "\n", " # Execute CQT extraction\n", " print(\"Step 2: CQT patch extraction and realignment\")\n", " \n", " _, fk, dummy_ts = CQT_custom(np.zeros(1000), fs, F_MIN, F_MAX, F_PO)\n", " time_res = dummy_ts[1] - dummy_ts[0]\n", " n_cols_patch = int(round(PATCH_WIDTH_S / time_res))\n", " n_samples_win = int(WIN_SIZE_CQT * fs)\n", "\n", " patches_data = [] \n", " scalars_data = [] \n", " rejected_clicks = [] \n", "\n", " for i, row in df_times.iterrows():\n", " t_clic = row['t_global']\n", " orig_idx = int(row['index'])\n", " \n", " feats = {'cqt_centroid': 0, 'cqt_bandwidth': 0, 'cqt_entropy': 0, 'cqt_energy': -120}\n", " patch_final = np.zeros((1,1))\n", " status = \"empty\"\n", "\n", " try:\n", " if pd.notna(t_clic): \n", " idx_start = int((t_clic - WIN_SIZE_CQT/2) * fs)\n", " idx_end = idx_start + n_samples_win\n", " \n", " if idx_start >= 0 and idx_end <= len(audio_signal):\n", " segment = audio_signal[idx_start:idx_end]\n", " S, _, ts_seg = CQT_custom(segment, fs, F_MIN, F_MAX, F_PO)\n", " \n", " energy_prof = np.mean(S, axis=0)\n", " center_idx = len(energy_prof) // 2\n", " win_rec_pix = int(REALIGN_WIN_S / time_res)\n", " \n", " s_r = max(0, center_idx - win_rec_pix)\n", " e_r = min(len(energy_prof), center_idx + win_rec_pix)\n", " \n", " if e_r > s_r:\n", " offset = np.argmax(energy_prof[s_r:e_r])\n", " local_peak = s_r + offset\n", " \n", " p_s = local_peak - n_cols_patch // 2\n", " p_e = p_s + n_cols_patch\n", " \n", " if p_s >= 0 and p_e <= S.shape[1]:\n", " patch_final = S[:, p_s:p_e].astype(np.float32)\n", " feats = extract_cqt_scalars(patch_final, fk)\n", " status = \"ok\"\n", " else:\n", " raise ValueError(f\"Time out of audio bounds (t={t_clic:.2f}s)\")\n", "\n", " if status == \"ok\":\n", " patches_data.append({\n", " \"index\": orig_idx, \n", " \"t_clic\": t_clic,\n", " \"patch\": patch_final,\n", " \"status\": status\n", " })\n", " \n", " row_csv = feats.copy()\n", " row_csv['index'] = orig_idx \n", " row_csv['T_H0'] = t_clic \n", " row_csv['status_cqt'] = status\n", " scalars_data.append(row_csv)\n", "\n", " except Exception as e:\n", " rejected_clicks.append({\"idx\": orig_idx, \"err\": str(e)})\n", "\n", " print(f\"Extraction complete. Valid patches: {len(patches_data)}\")\n", " \n", " if len(patches_data) == 0 and len(rejected_clicks) > 0:\n", " print(f\"Warning: no valid clicks extracted. First error: {rejected_clicks[0]['err']}\")\n", "\n", " # Save outputs dynamically using the extracted base_name\n", " print(\"Step 3: saving CQT tensors and features\")\n", " \n", " pkl_name = os.path.join(csv_dir, f\"{base_name}_CQT_patch_H{TARGET_CHANNEL}.pkl\")\n", " data_pkl = {\n", " \"patches\": patches_data,\n", " \"f_axis\": fk,\n", " \"params\": {\"fs\": fs, \"f_min\": F_MIN, \"channel\": TARGET_CHANNEL}\n", " }\n", " with open(pkl_name, \"wb\") as f:\n", " pickle.dump(data_pkl, f)\n", " print(f\"Saved patches (Pickle): {os.path.basename(pkl_name)}\")\n", " \n", " csv_name = os.path.join(csv_dir, f\"{base_name}_CQT_feat_H{TARGET_CHANNEL}.csv\")\n", " df_feats = pd.DataFrame(scalars_data)\n", " df_feats.to_csv(csv_name, index=False, sep=';')\n", " print(f\"Saved features (CSV): {os.path.basename(csv_name)}\")\n", "\n", "print(\"\\n--- BATCH PROCESSING COMPLETE ---\")" ] }, { "cell_type": "code", "execution_count": null, "id": "657c470a-ec6f-4dfa-84df-a8da492af1f4", "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 }