Plot Ica Blind Source SeparationΒΆ

===================================== Blind source separation using FastICAΒΆ

An example of estimating sources from noisy data.

ref:

ICA is used to estimate sources given noisy measurements. Imagine 3 instruments playing simultaneously and 3 microphones recording the mixed signals. ICA is used to recover the sources ie. what is played by each instrument. Importantly, PCA fails at recovering our instruments since the related signals reflect non-Gaussian processes.

Imports for Blind Source Separation with FastICAΒΆ

The cocktail party problem: Independent Component Analysis (ICA) solves the blind source separation problem: given mixed observations X = S * A^T (where S contains unknown independent source signals and A is an unknown mixing matrix), ICA recovers both S and A without any knowledge of either. FastICA uses a fixed-point iteration algorithm that maximizes the non-Gaussianity of the estimated sources, exploiting the Central Limit Theorem principle that a mixture of independent signals tends toward Gaussian distribution. The three source signals here – sinusoidal, square wave (via np.sign), and sawtooth (via signal.sawtooth) – are all highly non-Gaussian, making them ideal for ICA recovery.

Why PCA fails at source separation: PCA finds orthogonal directions of maximum variance, but the true independent sources need not be orthogonal in the observation space. PCA’s components are ordered by variance, not by statistical independence, so it rotates the data to align with variance axes but cannot undo the non-orthogonal mixing. The whiten="arbitrary-variance" parameter in FastICA first spheres the data (whitening), then finds the rotation that maximizes non-Gaussianity. The mixing_ attribute recovers the estimated mixing matrix A, and the assertion np.allclose(X, S_ @ A_.T + ica.mean_) verifies that the decomposition is self-consistent.

# Authors: The scikit-learn developers
# SPDX-License-Identifier: BSD-3-Clause

# %%
# Generate sample data
# --------------------

import numpy as np
from scipy import signal

np.random.seed(0)
n_samples = 2000
time = np.linspace(0, 8, n_samples)

s1 = np.sin(2 * time)  # Signal 1 : sinusoidal signal
s2 = np.sign(np.sin(3 * time))  # Signal 2 : square signal
s3 = signal.sawtooth(2 * np.pi * time)  # Signal 3: saw tooth signal

S = np.c_[s1, s2, s3]
S += 0.2 * np.random.normal(size=S.shape)  # Add noise

S /= S.std(axis=0)  # Standardize data
# Mix data
A = np.array([[1, 1, 1], [0.5, 2, 1.0], [1.5, 1.0, 2.0]])  # Mixing matrix
X = np.dot(S, A.T)  # Generate observations

# %%
# Fit ICA and PCA models
# ----------------------

from sklearn.decomposition import PCA, FastICA

# Compute ICA
ica = FastICA(n_components=3, whiten="arbitrary-variance")
S_ = ica.fit_transform(X)  # Reconstruct signals
A_ = ica.mixing_  # Get estimated mixing matrix

# We can `prove` that the ICA model applies by reverting the unmixing.
assert np.allclose(X, np.dot(S_, A_.T) + ica.mean_)

# For comparison, compute PCA
pca = PCA(n_components=3)
H = pca.fit_transform(X)  # Reconstruct signals based on orthogonal components

# %%
# Plot results
# ------------

import matplotlib.pyplot as plt

plt.figure()

models = [X, S, S_, H]
names = [
    "Observations (mixed signal)",
    "True Sources",
    "ICA recovered signals",
    "PCA recovered signals",
]
colors = ["red", "steelblue", "orange"]

for ii, (model, name) in enumerate(zip(models, names), 1):
    plt.subplot(4, 1, ii)
    plt.title(name)
    for sig, color in zip(model.T, colors):
        plt.plot(sig, color=color)

plt.tight_layout()
plt.show()