{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "cell-0",
   "metadata": {},
   "source": [
    "# Physiological Features\n",
    "\n",
    "This notebook demonstrates extraction of frequency-domain, nonlinear, and entropy-based physiological features from ECG and PPG 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 signal length: {len(signal_col)} samples ({len(signal_col)/fs:.1f}s at {fs} Hz)\")\n",
    "print(f\"PPG signal length: {len(ppg_col)} samples\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-3",
   "metadata": {},
   "source": [
    "## Nonlinear Features"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-4",
   "metadata": {},
   "outputs": [],
   "source": [
    "from vitalDSP.physiological_features.nonlinear import NonlinearFeatures\n",
    "\n",
    "nl = NonlinearFeatures(signal_col, fs=fs)\n",
    "\n",
    "sample_entropy = nl.compute_sample_entropy(m=2, r=0.2)\n",
    "app_entropy    = nl.compute_approximate_entropy(m=2, r=0.2)\n",
    "dfa_alpha      = nl.compute_dfa()\n",
    "lyapunov       = nl.compute_lyapunov_exponent()\n",
    "\n",
    "print(\"Nonlinear Features:\")\n",
    "print(f\"  Sample Entropy      : {sample_entropy:.4f}\")\n",
    "print(f\"  Approximate Entropy : {app_entropy:.4f}\")\n",
    "print(f\"  DFA Alpha           : {dfa_alpha:.4f}\")\n",
    "print(f\"  Lyapunov Exponent   : {lyapunov:.4f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-5",
   "metadata": {},
   "source": [
    "## Poincar\u00c3\u00a9 Plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-6",
   "metadata": {},
   "outputs": [],
   "source": [
    "poincare = nl.compute_poincare_features()\n",
    "print(\"Poincar\u00c3\u00a9 Features:\", poincare)\n",
    "\n",
    "# Visualize Poincar\u00c3\u00a9 plot\n",
    "fig = go.Figure()\n",
    "fig.add_trace(go.Scatter(\n",
    "    x=signal_col[:-1], y=signal_col[1:],\n",
    "    mode='markers', marker=dict(size=3, opacity=0.5),\n",
    "    name='Poincar\u00c3\u00a9'\n",
    "))\n",
    "fig.update_layout(\n",
    "    title='Poincar\u00c3\u00a9 Plot',\n",
    "    xaxis_title='x(n)', yaxis_title='x(n+1)'\n",
    ")\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-7",
   "metadata": {},
   "source": [
    "## Symbolic Dynamics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-8",
   "metadata": {},
   "outputs": [],
   "source": [
    "from vitalDSP.physiological_features.symbolic_dynamics import SymbolicDynamics\n",
    "\n",
    "sd = SymbolicDynamics(signal_col, n_symbols=4, word_length=3)\n",
    "\n",
    "perm_entropy   = sd.compute_permutation_entropy()\n",
    "shannon_entropy = sd.compute_shannon_entropy()\n",
    "word_dist      = sd.compute_word_distribution()\n",
    "\n",
    "print(\"Symbolic Dynamics:\")\n",
    "print(f\"  Permutation Entropy : {perm_entropy:.4f}\")\n",
    "print(f\"  Shannon Entropy     : {shannon_entropy:.4f}\")\n",
    "\n",
    "# Word distribution bar chart\n",
    "words = list(word_dist.keys())\n",
    "counts = list(word_dist.values())\n",
    "fig = go.Figure(go.Bar(x=[str(w) for w in words], y=counts))\n",
    "fig.update_layout(title='Symbolic Word Distribution', xaxis_title='Word', yaxis_title='Count')\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-9",
   "metadata": {},
   "source": [
    "## Envelope Detection"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-10",
   "metadata": {},
   "outputs": [],
   "source": [
    "from vitalDSP.physiological_features.envelope_detection import EnvelopeDetection\n",
    "\n",
    "ed = EnvelopeDetection(signal_col)\n",
    "hilbert_env = ed.hilbert_envelope()\n",
    "rms_env     = ed.moving_average_envelope(window_size=32)\n",
    "peak_env    = ed.peak_envelope()\n",
    "\n",
    "fig = go.Figure()\n",
    "fig.add_trace(go.Scatter(y=signal_col[:512], mode='lines', name='ECG Signal', opacity=0.4))\n",
    "fig.add_trace(go.Scatter(y=hilbert_env[:512], mode='lines', name='Hilbert Envelope'))\n",
    "fig.add_trace(go.Scatter(y=rms_env[:512], mode='lines', name='RMS Envelope'))\n",
    "fig.add_trace(go.Scatter(y=peak_env[:512], mode='lines', name='Peak Envelope'))\n",
    "fig.update_layout(title='Signal Envelopes', xaxis_title='Sample', yaxis_title='Amplitude')\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-11",
   "metadata": {},
   "source": [
    "## Trend Analysis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-12",
   "metadata": {},
   "outputs": [],
   "source": [
    "from vitalDSP.physiological_features.trend_analysis import TrendAnalysis\n",
    "\n",
    "ta = TrendAnalysis(signal_col)\n",
    "ema_trend    = ta.compute_exponential_smoothing(alpha=0.1)\n",
    "ma_trend     = ta.compute_moving_average(window_size=16)\n",
    "linear_trend = ta.compute_linear_trend()\n",
    "\n",
    "fig = go.Figure()\n",
    "fig.add_trace(go.Scatter(y=signal_col[:512], mode='lines', name='ECG Signal', opacity=0.4))\n",
    "fig.add_trace(go.Scatter(y=ema_trend[:512], mode='lines', name='EMA Trend'))\n",
    "fig.add_trace(go.Scatter(y=ma_trend[:512], mode='lines', name='Moving Average'))\n",
    "fig.add_trace(go.Scatter(y=linear_trend[:512], mode='lines', name='Linear Trend'))\n",
    "fig.update_layout(title='Signal Trend Analysis', xaxis_title='Sample', yaxis_title='Amplitude')\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-13",
   "metadata": {},
   "source": [
    "## Signal Change Detection"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-14",
   "metadata": {},
   "outputs": [],
   "source": [
    "from vitalDSP.physiological_features.signal_change_detection import SignalChangeDetection\n",
    "\n",
    "scd = SignalChangeDetection(signal_col)\n",
    "zcr       = scd.zero_crossing_rate()\n",
    "var_pts   = scd.variance_based_detection(window_size=64, threshold=2.0)\n",
    "\n",
    "print(f\"Zero Crossing Rate          : {zcr:.4f}\")\n",
    "print(f\"Variance-based change points: {len(var_pts)}\")\n",
    "\n",
    "fig = go.Figure()\n",
    "fig.add_trace(go.Scatter(y=signal_col[:512], mode='lines', name='ECG Signal'))\n",
    "cp_in_range = [p for p in var_pts if p < 512]\n",
    "if cp_in_range:\n",
    "    fig.add_trace(go.Scatter(\n",
    "        x=cp_in_range,\n",
    "        y=[signal_col[p] for p in cp_in_range],\n",
    "        mode='markers', marker=dict(color='red', size=8),\n",
    "        name='Change Points'\n",
    "    ))\n",
    "fig.update_layout(title='Signal Change Detection', xaxis_title='Sample', yaxis_title='Amplitude')\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-15",
   "metadata": {},
   "source": [
    "## Signal Segmentation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-16",
   "metadata": {},
   "outputs": [],
   "source": [
    "from vitalDSP.physiological_features.signal_segmentation import SignalSegmentation\n",
    "\n",
    "seg = SignalSegmentation(signal_col)\n",
    "\n",
    "# Fixed-size segments\n",
    "fixed_segs = seg.fixed_size_segmentation(segment_size=128)\n",
    "print(f\"Fixed-size segments (128 samples): {len(fixed_segs)}\")\n",
    "\n",
    "# Variance-based adaptive segmentation\n",
    "var_segs = seg.variance_based_segmentation(window_size=64, threshold_factor=1.5)\n",
    "print(f\"Variance-based segments: {len(var_segs)}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-17",
   "metadata": {},
   "source": [
    "## Energy Analysis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-18",
   "metadata": {},
   "outputs": [],
   "source": [
    "from vitalDSP.physiological_features.energy_analysis import EnergyAnalysis\n\nea = EnergyAnalysis(signal_col, fs=fs)\ntotal_energy  = ea.compute_total_energy()\nband_energy   = ea.compute_band_energy(low_freq=0.5, high_freq=40.0)\nspectral_energy = ea.compute_spectral_energy()\n\nprint(\"Energy Analysis:\")\nprint(f\"  Total Energy          : {total_energy:.4f}\")\nprint(f\"  Band Energy (0.5-40Hz): {band_energy:.4f}\")\nprint(f\"  Spectral Energy       : {spectral_energy:.4f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-19",
   "metadata": {},
   "source": [
    "## Signal Power Analysis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-20",
   "metadata": {},
   "outputs": [],
   "source": [
    "from vitalDSP.physiological_features.signal_power_analysis import SignalPowerAnalysis\n",
    "\n",
    "spa = SignalPowerAnalysis(signal_col)\n",
    "total_power = spa.compute_total_power()\n",
    "rms_power   = spa.compute_mean_square()\n",
    "snr         = spa.compute_snr()\n",
    "\n",
    "print(\"Signal Power Analysis:\")\n",
    "print(f\"  Total Power : {total_power:.4f}\")\n",
    "print(f\"  RMS Power   : {rms_power:.4f}\")\n",
    "print(f\"  SNR (dB)    : {snr:.4f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Nonlinear Features \u2014 Fractal Dimension and Recurrence"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from vitalDSP.physiological_features.nonlinear import NonlinearFeatures\n\nnl = NonlinearFeatures(signal_col, fs=fs)\n\nfractal_dim = nl.compute_fractal_dimension()\nrecurrence  = nl.compute_recurrence_features(threshold=0.2, sample_size=200)\n\nprint(f'Fractal Dimension  : {fractal_dim:.4f}')\nprint('Recurrence Features:')\nfor k, v in recurrence.items():\n    print(f'  {k}: {v:.4f}' if isinstance(v, float) else f'  {k}: {v}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Symbolic Dynamics \u2014 Full Pipeline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from vitalDSP.physiological_features.symbolic_dynamics import SymbolicDynamics\n\nsd = SymbolicDynamics(signal_col, n_symbols=4, word_length=3)\n\nsymbols      = sd.symbolize()\nrenyi        = sd.compute_renyi_entropy(alpha=2)\nforbidden    = sd.detect_forbidden_words()\ntrans_matrix = sd.compute_transition_matrix()\nall_features = sd.compute_symbolic_features()\n\nprint(f'Symbolized sequence (first 20): {symbols[:20]}')\nprint(f'Renyi Entropy (alpha=2)        : {renyi:.4f}')\nprint(f'Forbidden words count          : {len(forbidden)}')\nprint(f'Transition matrix shape        : {trans_matrix.shape}')\nprint('All symbolic features:', all_features)\n\nfig = go.Figure(go.Heatmap(z=trans_matrix, colorscale='Blues', colorbar=dict(title='Prob')))\nfig.update_layout(title='Symbol Transition Matrix',\n                  xaxis_title='Next Symbol', yaxis_title='Current Symbol')\nfig.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Envelope Detection \u2014 Extended Methods"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from vitalDSP.physiological_features.envelope_detection import EnvelopeDetection\n\ned = EnvelopeDetection(signal_col)\nabs_env     = ed.absolute_value_envelope(smoothing_factor=16)\nwavelet_env = ed.wavelet_envelope(wavelet_name='db4', level=2)\n\nfig = go.Figure()\nfig.add_trace(go.Scatter(y=signal_col[:512], mode='lines', name='ECG Signal', opacity=0.4))\nfig.add_trace(go.Scatter(y=abs_env[:512], mode='lines', name='Absolute Value Envelope'))\nfig.add_trace(go.Scatter(y=wavelet_env[:512], mode='lines', name='Wavelet Envelope (db4)'))\nfig.update_layout(title='Extended Envelope Methods',\n                  xaxis_title='Sample', yaxis_title='Amplitude')\nfig.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Trend Analysis \u2014 Extended Methods"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from vitalDSP.physiological_features.trend_analysis import TrendAnalysis\n\nta = TrendAnalysis(signal_col)\n\npoly_trend     = ta.compute_polynomial_trend(degree=3)\nmomentum       = ta.compute_momentum(window_size=16)\ntrend_strength = ta.compute_trend_strength()\nreversals      = ta.detect_trend_reversal(window_size=32)\n\nprint(f'Trend Strength  : {trend_strength:.6f}')\nprint(f'Trend Reversals : {len(reversals)} detected')\n\nfig = go.Figure()\nfig.add_trace(go.Scatter(y=signal_col[:512], mode='lines', name='ECG Signal', opacity=0.4))\nfig.add_trace(go.Scatter(y=poly_trend[:512], mode='lines', name='Polynomial Trend (deg=3)'))\nfig.add_trace(go.Scatter(y=momentum[:512], mode='lines', name='Momentum'))\nrev_in = [r for r in reversals if r < 512]\nif rev_in:\n    fig.add_trace(go.Scatter(x=rev_in, y=[signal_col[r] for r in rev_in],\n                             mode='markers', marker=dict(color='red', size=8),\n                             name='Trend Reversals'))\nfig.update_layout(title='Trend Analysis Extended',\n                  xaxis_title='Sample', yaxis_title='Amplitude')\nfig.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Signal Change Detection \u2014 Extended Methods"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from vitalDSP.physiological_features.signal_change_detection import SignalChangeDetection\n\nscd = SignalChangeDetection(signal_col)\n\nabs_diff   = scd.absolute_difference()\nenergy_pts = scd.energy_based_detection(window_size=64)\nadapt_pts  = scd.adaptive_threshold_detection(threshold_factor=1.5, window_size=64)\n\nprint(f'Absolute difference mean  : {abs_diff.mean():.4f}')\nprint(f'Energy-based change points: {len(energy_pts)}')\nprint(f'Adaptive threshold points : {len(adapt_pts)}')\n\nfig = go.Figure()\nfig.add_trace(go.Scatter(y=abs_diff[:512], mode='lines', name='Absolute Difference'))\nfig.update_layout(title='Absolute Difference Signal',\n                  xaxis_title='Sample', yaxis_title='|delta|')\nfig.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Signal Segmentation \u2014 Extended Methods"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from vitalDSP.physiological_features.signal_segmentation import SignalSegmentation\n\nseg = SignalSegmentation(signal_col)\n\nthresh_segs = seg.threshold_based_segmentation(threshold=0.5)\npeak_segs   = seg.peak_based_segmentation(min_distance=32, height=0.5)\n\nprint(f'Threshold-based segments : {len(thresh_segs)}')\nprint(f'Peak-based segments      : {len(peak_segs)}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Energy Analysis \u2014 Extended Methods"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from vitalDSP.physiological_features.energy_analysis import EnergyAnalysis\n\nea = EnergyAnalysis(signal_col, fs=fs)\n\nseg_energy = ea.compute_segment_energy(start_idx=0, end_idx=256)\nqrs_energy = ea.compute_qrs_energy(r_peaks=[100, 220, 340, 460])\n\nprint(f'Segment energy (first 256 samples) : {seg_energy:.4f}')\nprint(f'QRS energy (4 synthetic peaks)     : {qrs_energy:.4f}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Signal Power Analysis \u2014 Extended Methods"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from vitalDSP.physiological_features.signal_power_analysis import SignalPowerAnalysis\n\nspa = SignalPowerAnalysis(signal_col)\n\nrmse       = spa.compute_rmse()\npeak_power = spa.compute_peak_power()\nfreqs, psd = spa.compute_psd(fs=fs, nperseg=256)\nband_power = spa.compute_band_power(band=(0.5, 40.0), fs=fs, nperseg=256)\n\nprint(f'RMSE                  : {rmse:.4f}')\nprint(f'Peak Power            : {peak_power:.4f}')\nprint(f'Band Power (0.5-40Hz) : {band_power:.4f}')\n\nfig = go.Figure()\nfig.add_trace(go.Scatter(x=freqs[:60], y=psd[:60], mode='lines', name='PSD'))\nfig.update_layout(title='Power Spectral Density',\n                  xaxis_title='Frequency (Hz)', yaxis_title='Power')\nfig.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Frequency Domain HRV Features"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from vitalDSP.physiological_features.frequency_domain import FrequencyDomainFeatures\nfrom vitalDSP.physiological_features.beat_to_beat import BeatToBeatAnalysis\n\nbb = BeatToBeatAnalysis(signal_col, fs=fs, signal_type='ECG')\nrr_intervals = bb.compute_rr_intervals()\nprint(f'RR intervals: {len(rr_intervals)} beats')\n\nfdf = FrequencyDomainFeatures(rr_intervals)\npsd_dict = fdf.compute_psd()\n\nprint('Frequency Domain HRV Features:')\nprint(f'  ULF power   : {fdf.compute_ulf():.4f} ms2')\nprint(f'  VLF power   : {fdf.compute_vlf():.4f} ms2')\nprint(f'  LF power    : {fdf.compute_lf():.4f} ms2')\nprint(f'  HF power    : {fdf.compute_hf():.4f} ms2')\nprint(f'  LF/HF ratio : {fdf.compute_lf_hf_ratio():.4f}')\nprint(f'  Total power : {fdf.compute_total_power():.4f} ms2')\nprint(f'  LFnu        : {fdf.compute_lfnu():.4f}')\nprint(f'  HFnu        : {fdf.compute_hfnu():.4f}')\n\nfig = go.Figure()\nfig.add_trace(go.Scatter(x=psd_dict['frequencies'], y=psd_dict['psd'],\n                         mode='lines', name='HRV PSD'))\nfig.update_layout(title='HRV Power Spectral Density',\n                  xaxis_title='Frequency (Hz)', yaxis_title='Power (ms2)')\nfig.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "name": "python",
   "version": "3.9.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}