{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "cell-0",
   "metadata": {},
   "source": [
    "# Signal Decomposition and Spectral Features\n",
    "\n",
    "This notebook demonstrates PCA/ICA signal decomposition, Wavelet-FFT fusion, and MFCC feature extraction for physiological signals using vitalDSP."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-1",
   "metadata": {},
   "source": [
    "## Setup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-2",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import plotly.io as pio\n",
    "pio.renderers.default = \"sphinx_gallery\"\n",
    "from plotly import graph_objects as go\n",
    "from vitalDSP.notebooks import load_sample_ecg_small, load_sample_ppg, plot_trace\n",
    "\n",
    "fs = 128\n",
    "signal_col, date_col = load_sample_ecg_small()\n",
    "signal_col = np.array(signal_col)\n",
    "\n",
    "ppg_col, ppg_date = load_sample_ppg()\n",
    "ppg_col = np.array(ppg_col)\n",
    "\n",
    "print(f\"ECG: {len(signal_col)} samples at {fs} Hz\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-3",
   "metadata": {},
   "source": [
    "## PCA Signal Decomposition\n",
    "\n",
    "`PCASignalDecomposition` computes eigenvectors of the signal covariance matrix. Project the multi-channel signals onto these eigenvectors to get the principal component time series."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-4",
   "metadata": {},
   "outputs": [],
   "source": [
    "from vitalDSP.transforms.pca_ica_signal_decomposition import PCASignalDecomposition\n",
    "\n",
    "# Short segment: PCA computes a (n_samples x n_samples) covariance matrix internally\n",
    "n = 256\n",
    "signals = np.stack([\n",
    "    signal_col[:n],\n",
    "    ppg_col[:n],\n",
    "    signal_col[:n] * 0.8 + ppg_col[:n] * 0.2,\n",
    "], axis=0)  # shape: (3 channels, n samples)\n",
    "\n",
    "pca = PCASignalDecomposition(signals, n_components=3)\n",
    "# compute_pca returns the loading matrix (eigenvectors) shape (n_components, n_channels)\n",
    "loading_matrix = pca.compute_pca()\n",
    "print(f\"Channels input : {signals.shape}\")\n",
    "print(f\"Loading matrix : {loading_matrix.shape}\")\n",
    "\n",
    "# Project original channels onto principal components -> (n_samples, n_components)\n",
    "pc_timeseries = signals.T @ loading_matrix.T  # (256, 3)\n",
    "print(f\"PC time series : {pc_timeseries.shape}\")\n",
    "\n",
    "fig = go.Figure()\n",
    "for i in range(pc_timeseries.shape[1]):\n",
    "    fig.add_trace(go.Scatter(y=pc_timeseries[:, i], mode='lines', name=f'PC {i+1}'))\n",
    "fig.update_layout(\n",
    "    title='PCA Principal Components',\n",
    "    xaxis_title='Sample', yaxis_title='Amplitude'\n",
    ")\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-5",
   "metadata": {},
   "source": [
    "## ICA Signal Decomposition\n",
    "\n",
    "ICA (Independent Component Analysis) recovers statistically independent sources from a linear mixture."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-6",
   "metadata": {},
   "outputs": [],
   "source": [
    "from vitalDSP.transforms.pca_ica_signal_decomposition import ICASignalDecomposition\n",
    "\n",
    "# ICA uses FastICA — whitening + deflation; result is the unmixing weight matrix\n",
    "ica = ICASignalDecomposition(signals, max_iter=500, tolerance=1e-3)\n",
    "try:\n",
    "    ica_weights = ica.compute_ica()  # shape: (n_components, n_components)\n",
    "    # Recover independent components: project whitened signals through unmixing matrix\n",
    "    whitened = (signals - signals.mean(axis=1, keepdims=True)) / (signals.std(axis=1, keepdims=True) + 1e-8)\n",
    "    ics = ica_weights @ whitened  # (n_components, n_samples)\n",
    "    print(f\"ICA weights : {ica_weights.shape}\")\n",
    "    print(f\"IC time series: {ics.shape}\")\n",
    "    fig = go.Figure()\n",
    "    for i in range(ics.shape[0]):\n",
    "        fig.add_trace(go.Scatter(y=ics[i], mode='lines', name=f'IC {i+1}'))\n",
    "    fig.update_layout(\n",
    "        title='ICA Independent Components',\n",
    "        xaxis_title='Sample', yaxis_title='Amplitude'\n",
    "    )\n",
    "    fig.show()\n",
    "except Exception as e:\n",
    "    print(f\"ICA did not converge for this data segment: {e}\")\n",
    "    print(\"Try increasing max_iter or using a longer segment.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-7",
   "metadata": {},
   "source": [
    "## Wavelet-FFT Fusion"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-8",
   "metadata": {},
   "outputs": [],
   "source": [
    "from vitalDSP.transforms.wavelet_fft_fusion import WaveletFFTfusion\n",
    "\n",
    "# Fuse wavelet time-frequency localization with FFT spectral analysis\n",
    "fusion = WaveletFFTfusion(signal_col, wavelet_type='db4', order=4)\n",
    "fusion_result = fusion.compute_fusion()\n",
    "\n",
    "print(f\"Fusion result shape: {fusion_result.shape}  ({fusion_result.shape[0]} levels x {fusion_result.shape[1]} samples)\")\n",
    "\n",
    "# Each row: wavelet detail/approx coefficients scaled by corresponding FFT magnitude\n",
    "fig = go.Figure()\n",
    "for i in range(fusion_result.shape[0]):\n",
    "    fig.add_trace(go.Scatter(\n",
    "        y=np.abs(fusion_result[i, :256]), mode='lines',\n",
    "        name=f'Level {i+1} |fusion|'\n",
    "    ))\n",
    "fig.update_layout(\n",
    "    title='Wavelet-FFT Fusion (magnitude per level)',\n",
    "    xaxis_title='Sample', yaxis_title='|Fusion Coefficient|'\n",
    ")\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-9",
   "metadata": {},
   "source": [
    "## Wavelet-FFT Fusion — Level Energy by Wavelet Type"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-10",
   "metadata": {},
   "outputs": [],
   "source": [
    "wavelets = ['haar', 'db2', 'db4']\n",
    "fig = go.Figure()\n",
    "\n",
    "for wt in wavelets:\n",
    "    f = WaveletFFTfusion(signal_col, wavelet_type=wt, order=3)\n",
    "    result = f.compute_fusion()\n",
    "    energy_per_level = [np.mean(np.abs(result[i])**2) for i in range(result.shape[0])]\n",
    "    fig.add_trace(go.Bar(\n",
    "        name=wt,\n",
    "        x=[f'L{i+1}' for i in range(len(energy_per_level))],\n",
    "        y=energy_per_level\n",
    "    ))\n",
    "\n",
    "fig.update_layout(\n",
    "    title='Wavelet-FFT Fusion Energy per Level by Wavelet Type',\n",
    "    xaxis_title='Decomposition Level', yaxis_title='Mean Energy',\n",
    "    barmode='group'\n",
    ")\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-11",
   "metadata": {},
   "source": [
    "## MFCC Feature Extraction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-12",
   "metadata": {},
   "outputs": [],
   "source": [
    "from vitalDSP.transforms.mfcc import MFCC\n",
    "\n",
    "# MFCC captures spectral envelope — here applied to ECG for demonstration\n",
    "mfcc = MFCC(signal_col, sample_rate=fs, num_filters=26, num_coefficients=13)\n",
    "mfcc_result = mfcc.compute_mfcc()\n",
    "\n",
    "print(f\"MFCC result shape: {mfcc_result.shape}  (frames x coefficients)\")\n",
    "\n",
    "fig = go.Figure(go.Heatmap(\n",
    "    z=mfcc_result.T,\n",
    "    colorscale='Viridis',\n",
    "    colorbar=dict(title='Value')\n",
    "))\n",
    "fig.update_layout(\n",
    "    title='MFCC Cepstrogram',\n",
    "    xaxis_title='Frame', yaxis_title='Coefficient Index'\n",
    ")\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-13",
   "metadata": {},
   "source": [
    "## MFCC Per-Coefficient Statistics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-14",
   "metadata": {},
   "outputs": [],
   "source": [
    "mfcc_means = np.mean(mfcc_result, axis=0)\n",
    "mfcc_stds  = np.std(mfcc_result, axis=0)\n",
    "coeff_idx  = list(range(1, len(mfcc_means) + 1))\n",
    "\n",
    "fig = go.Figure()\n",
    "fig.add_trace(go.Bar(\n",
    "    x=coeff_idx, y=mfcc_means,\n",
    "    error_y=dict(type='data', array=mfcc_stds, visible=True),\n",
    "    name='Mean ± Std'\n",
    "))\n",
    "fig.update_layout(\n",
    "    title='MFCC Coefficient Statistics',\n",
    "    xaxis_title='Coefficient Index', yaxis_title='Value'\n",
    ")\n",
    "fig.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "name": "python",
   "version": "3.9.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}