Modulation Index

cfc.pac.run_mi(low_freq_data, high_freq_data, phase_window_half_size=10, phase_step_width=20)

Calculate the modulation index between a low frequency signal and a high frequency signal.

Parameters:

• low_freq_data (np.ndarray, shape(n_samples,)) – Single array of low frequency data.

• high_freq_data (np.ndarray, shape(n_samples,)) – Single array of high frequency data.

• phase_window_half_size (int) – Width of the phase window used for calculation of frequency/phase histogram. Amplitude gets added to every phase bin within the window size. Larger windows result in more smooth/increased PAC estimates.

• phase_step_width (int) – Step width/shift of the phase window used for calculation of frequency/phase histogram.

Returns:

Amount of phase amplitude coupling measured using the modulation index.

Return type:

float

The following example shows how to apply the mi to estimate PAC.

import numpy as np

import finn.cfc.pac as pac

def generate_high_frequency_signal(n, frequency_sampling, frequency_within_bursts, random_noise_strength,
                                   offset, burst_count, burst_length):
    signal = np.random.normal(0, 1, n) * random_noise_strength

    for burst_start in np.arange(offset, n, n/burst_count):
        burst_end = burst_start + (burst_length/2)
        signal[int(burst_start):int(burst_end)] =  np.sin(2 * np.pi * frequency_within_bursts * np.arange(0, (int(burst_end) - int(burst_start))) / frequency_sampling)

    return signal

def main():
    #Configure sample data
    data_range = np.arange(0, 10000)
    frequency_sampling = 1000
    frequencies_between_bursts = [2, 5, 10, 15, 50]
    tgt_frequency_between_bursts = 10
    frequency_within_bursts = 200
    high_freq_frame_offsets = [0, 20, 40]

    #Configure noise data
    random_noise_strength = 0.3

    #Generate sample data
    burst_length = frequency_sampling / tgt_frequency_between_bursts
    burst_count = len(data_range) / frequency_sampling * tgt_frequency_between_bursts

    high_freq_signals = [generate_high_frequency_signal(len(data_range), frequency_sampling, frequency_within_bursts,
                                                        random_noise_strength, high_freq_frame_offset, burst_count, burst_length) for high_freq_frame_offset in high_freq_frame_offsets]
    low_freq_signals = [np.sin(2 * np.pi * frequency_between_bursts * data_range / frequency_sampling) for frequency_between_bursts in frequencies_between_bursts]

    scores = np.zeros((len(high_freq_signals), len(low_freq_signals)));
    for (high_freq_idx, high_freq_signal) in enumerate(high_freq_signals):
        for (low_freq_idx, low_freq_signal) in enumerate(low_freq_signals):
            scores[high_freq_idx, low_freq_idx] = pac.run_mi(low_freq_signal, high_freq_signal,
                                                             phase_window_half_size = 20, phase_step_width = 5)

    print("target frequency: ", tgt_frequency_between_bursts)
    for x in range(len(frequencies_between_bursts)):
        print("%.3f" % (frequencies_between_bursts[x]), end = "\t")
    print("")
    for y in range(len(high_freq_frame_offsets)):
        for x in range(len(frequencies_between_bursts)):
            print("%.3f" % (scores[y][x],), end = "\t")
        print("")

main()

Using the direct modulation index, PAC between the low frequency signal (2/5/10/15/50Hz) signal and the high frequency signal (10Hz amplitude modulation).

2Hz	5Hz	10Hz	15Hz	50Hz
0.002	0.005	0.004	0.001	0.001
0.003	0.004	0.006	0.001	0.001
0.002	0.004	0.005	0.001	0.001