Creating AverageTFR with and without canonical frequency bands

đź’¬ Support & Discussions
1
  • MNE version: 1.9
  • operating system: Windows 10

Hey everyone. Very very new to MNE Python. I have loaded epoched data that’s been cleaned and pre-processed. I now need to run time-frequency analysis using morlet wavelets on human scalp EEG data. I have 4 conditions: hits for scenes, misses for scenes, hits for objects, and misses for objects.

I am trying to run code with two approaches: 1) I want to have pre-defined frequency bands for the AverageTFR, 2) I dont want to let the data dictate where the Average TFR object is without pre-defined frequency bands. The reason is I will run spatio-temporal permutation testing on the TFR objects and it would be ideal if I could plot significant clusters without pre-set frequency ranges. Perhaps a cluster will appear between 5-10 Hz which is across theta and alpha. So I want to capture that data is well. Below is some code I’m working with. Thanks!


    # Compute PSD over channels and epochs
    spectrum, freqs = psd_welch(epochs, n_jobs=4, 
                                fmin=3, fmax=40, n_fft=500,
                                average='median', n_per_seg=250,
                                n_overlap=125)
    mean_spectrum = spectrum.mean(axis=(0, 1))
    
    
# Compute frequency bands
# Define canonical frequency bands (in Hz)
    lower = {
        'theta': 4,
        'alpha': 8,
        'beta': 16
    }
    
    upper = {
        'theta': 7,
        'alpha': 12,
        'beta': 26
    }
    
    # Define center frequencies for each band (optional but useful for labeling)
    center = {
        'theta': 5.5,
        'alpha': 10.5,
        'beta': 21.5
    }

    
    # Run morlet analysis
    queries = {
            'scene-hit65': ("category=='scene' and study_n_responses==1 and " +
                            "test_resp in [5,6]"),
            'scene-miss65': ("category=='scene' and study_n_responses==1 and " +
                             "test_resp in [1,2,3,4]"),
            'object-hit65': ("category=='object' and study_n_responses==1 " +
                             "and test_resp in [5,6]"),
            'object-miss65': ("category=='object' and study_n_responses==1 " +
                              "and test_resp in [1,2,3,4]"),
        }
    for key, value in queries.items():
        
        # Run Morlet TFR
        power, itc = tfr_morlet(epochs[value], freqs=np.arange(3, 60, 0.25),
                                n_cycles=5, decim=2, n_jobs=4)
        
        # Power_data will create a 3D object of channels x frequencies x time points
        # There were 761 time points. decim = 2 will cut them in half so that you are only
        # left with 376 time points. 
        # Get the bands
        theta = power.copy().crop(fmin=lower['theta'], fmax=upper['theta']).data.mean(axis=1)
        theta = np.expand_dims(theta, axis=1)
        alpha = power.copy().crop(fmin=lower['alpha'], fmax=upper['alpha']).data.mean(axis=1)
        alpha = np.expand_dims(alpha, axis=1)
        beta = power.copy().crop(fmin=lower['beta'], fmax=upper['beta']).data.mean(axis=1)
        beta = np.expand_dims(beta, axis=1)
        power_data = np.concatenate([theta, alpha, beta], axis=1)
        tmp_freqs = [center[x] for x in ['theta', 'alpha', 'beta']]
        
        # Make averageTFR object
        power_red = AverageTFR(power.info, power_data, power.times, tmp_freqs, power.nave,
                               comment=key, method=power.method)
        power_red.save(out_path / f"{row['id']}_{key}_tfr.h5", overwrite=True)
2

would it work for you if you did all of the analysis and statistics on the TFR object with freqs=np.arange(3, 60, 0.25) and only at the very end average frequencies into bands?

the main “unsolved” challenge I see here is that the clusters with an associated p-value < 0.05 (don’t call them “significant clusters” in your paper: How NOT to interpret results from a cluster-based permutation test - FieldTrip toolbox) are defined as tuples of integer indices. When averaging frequencies into bands, you need to collapse these integer indices as well so they correspond to the new bands instead of the original frequencies. MNE-Python has no functionality to directly do what you want, you will have to do this yourself. I think for this task, an AI assistant would be really helpful.