Applying SSP to CSD matrix

csd_linearity.py

"""
Test whether applying SSP to a CSD matrix is the same as computing the CSD
matrix on data where the SSP has already been applied.

author: Marijn van Vliet <[email protected]>
"""
import numpy as np
import mne

# Be quiet
mne.set_log_level(False)

# Get some epochs with and without SSP projections applied
path = mne.datasets.sample.data_path() + '/MEG/sample'
raw = mne.io.read_raw_fif(path + '/sample_audvis_raw.fif')
events = mne.find_events(raw)
picks = mne.pick_types(raw.info, meg='mag')
epochs_no_proj = mne.Epochs(raw, events, event_id=1, tmin=0, tmax=1,
                            picks=picks, preload=True, proj=False)
epochs_proj = epochs_no_proj.copy().apply_proj()


def test_equivalency(csd_no_proj, csd_proj):
    """Apply SSPs directly to CSD and see if the result is the same as
    having computed the CSD on data with the SSP applied."""
    
    # CSDs should start out different
    assert not np.allclose(csd_no_proj, csd_proj, atol=0)

    # Make and apply the projector
    proj, _, _ = mne.io.proj.make_projector(
        epochs_no_proj.info['projs'], epochs_no_proj.ch_names)
    csd_applied_proj = proj @ csd_no_proj @ proj.T

    # Is the result the same?
    print('csd_proj == csd_applied_proj?',
          np.allclose(csd_proj, csd_applied_proj, atol=0))


############################
# Computing CSD with morlets
############################
print('Testing csd_morlet')

# Compute the CSD on data without SSP projections
csd_no_proj = mne.time_frequency.csd_morlet(epochs_no_proj, [40]).get_data()

# Compute the CSD on data with SSP projections
csd_proj = mne.time_frequency.csd_morlet(epochs_proj, [40]).get_data()

test_equivalency(csd_no_proj, csd_proj)


############################
# Computing CSD with fourier
############################
print('Testing csd_fourier')

# Compute the CSD on data without SSP projections
csd_no_proj = mne.time_frequency.csd_fourier(
    epochs_no_proj, fmin=40, fmax=41).get_data()

# Compute the CSD on data with SSP projections
csd_proj = mne.time_frequency.csd_fourier(
    epochs_proj, fmin=40, fmax=41).get_data()

test_equivalency(csd_no_proj, csd_proj)


###############################
# Computing CSD with multitaper
###############################
print('Testing csd_multitaper')

# Compute the CSD on data without SSP projections
csd_no_proj = mne.time_frequency.csd_multitaper(
    epochs_no_proj, fmin=40, fmax=41).get_data()

# Compute the CSD on data with SSP projections
csd_proj = mne.time_frequency.csd_multitaper(
    epochs_proj, fmin=40, fmax=41).get_data()

test_equivalency(csd_no_proj, csd_proj)

Output:

Testing csd_morlet
csd_proj == csd_applied_proj? True
Testing csd_fourier
csd_proj == csd_applied_proj? True
Testing csd_multitaper
csd_proj == csd_applied_proj? True