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