Modulation Index

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

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

Parameters:

• low_freq_data – Single array of low frequency data.

• high_freq_data – Single array of high frequency data.

• phase_window_half_size – 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 – 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.

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