MEG信号时间序列正交化咨询：正弦信号测试异常问题

Hey there! Let's unpack what's going on with your orthogonalization test and explore better methods tailored for MEG signal space leakage.

Why Your Current Test Shows Such Big Differences #

The key issue here is your test signals: they’re perfectly linearly dependent — sig1 is just a scaled version of sig (200/325 ≈ 0.615). The spectral.orthogonalize function (which likely uses the Gram-Schmidt process) takes the first signal as a basis, then projects the second signal onto the orthogonal complement of the first. For completely dependent signals, this results in the second orthogonalized signal being nearly zero or completely out of phase with the original, hence the stark visual difference. This isn’t a flaw in the method—it’s exactly how orthogonalization works for linearly dependent vectors!

Better Orthogonalization Methods for MEG Signal Space Leakage #

MEG space leakage happens when signals from one brain region bleed into adjacent sensors, so we need methods that preserve the original signal’s temporal/spatial morphology while eliminating cross-talk. Here are the best approaches:

1. Signal Space Separation (SSS) — MEG’s Gold Standard #

SSS is purpose-built for MEG data. It decomposes signals into brain-generated (endogenous) and external/leakage (exogenous) components, retaining the core brain signal’s shape while suppressing unwanted cross-talk. You can implement it easily with MNE-Python:

import mne

# Load your raw MEG data (example)
raw = mne.io.read_raw_fif("your_meg_data.fif", preload=True)

# Apply Maxwell filtering (SSS implementation)
raw = mne.preprocessing.maxwell_filter(raw, origin="auto")

2. Normalized Gram-Schmidt Orthogonalization #

If you want to keep the signal’s temporal morphology intact, add normalization to the Gram-Schmidt process. This ensures orthogonalized signals retain their waveform shape (only amplitude/phase adjusts slightly):

import numpy as np

def normalized_gram_schmidt(X):
    X_ortho = np.zeros_like(X)
    for i in range(X.shape[0]):
        # Subtract projections from previous orthogonal components
        proj = np.sum(X_ortho[:i] * X[i], axis=1) @ X_ortho[:i]
        X_ortho[i] = X[i] - proj
        # Normalize to unit norm
        norm = np.linalg.norm(X_ortho[i])
        if norm > 1e-10:
            X_ortho[i] /= norm
    return X_ortho

# Test with linearly independent signals (fixing your original test case)
t = np.linspace(-0.02, 0.05, 1000)
sig = 325 * np.sin(2*np.pi*50*t)
sig1 = 200 * np.sin(2*np.pi*55*t + np.pi/4)  # Different frequency + phase shift
signal = np.vstack([sig, sig1])

ortho_norm = normalized_gram_schmidt(signal)

# Plot results
import matplotlib.pyplot as plt
plt.figure()
plt.plot(sig, "b", label="Original 1")
plt.plot(sig1, "g", label="Original 2")
plt.title("Original Linearly Independent Signals")
plt.legend()

plt.figure()
plt.plot(ortho_norm[0], "r", label="Orthogonalized 1")
plt.plot(ortho_norm[1], "k", label="Orthogonalized 2")
plt.title("Normalized Gram-Schmidt Results")
plt.legend()
plt.show()

3. PCA-Based Orthogonalization #

PCA projects signals onto orthogonal principal components that retain the maximum variance of the original data. For MEG, this reduces redundancy (including space leakage) while preserving the core signal morphology. It’s ideal for high-dimensional MEG datasets:

from sklearn.decomposition import PCA
import matplotlib.pyplot as plt

# Use the same linearly independent test signals as above
pca = PCA(n_components=2)
ortho_pca = pca.fit_transform(signal.T).T  # Transpose to fit PCA's expected shape

# Plot results
plt.figure()
plt.plot(ortho_pca[0], "r", label="PCA Component 1")
plt.plot(ortho_pca[1], "k", label="PCA Component 2")
plt.title("PCA Orthogonalized Signals")
plt.legend()
plt.show()

内容的提问来源于stack exchange，提问作者Vivek