Direct Modulation Index
=======================


.. currentmodule:: cfc.pac
.. autofunction:: run_dmi


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

.. code-block:: python
   
   import numpy as np
   import matplotlib
   matplotlib.use("Qt5agg")
   import matplotlib.pyplot as plt 
   
   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 visualize(scores, best_fits, amplitude_signals, frequencies_between_bursts, tgt_frequency_between_bursts, high_freq_frame_offsets):
      if(len(scores) < 5):
      (_, axes) = plt.subplots(2, 2)
      elif(len(scores) < 10):
      (_, axes) = plt.subplots(3, 3)
      elif(len(scores) < 17):
      (_, axes) = plt.subplots(4, 4)
      elif(len(scores) < 26):
      (_, axes) = plt.subplots(5, 5)
      else:
      raise NotImplementedError
      
      axes = np.reshape(axes, -1)
      
      for (ax_idx, ax) in enumerate(axes):
      ax.plot(np.arange(0, 1, 1/len(amplitude_signals[ax_idx])), amplitude_signals[ax_idx], label = "original data")
      ax.plot(np.arange(0, 1, 1/len(amplitude_signals[ax_idx])), best_fits[ax_idx], label = "fitted curve")
      if (frequencies_between_bursts[int(ax_idx % len(high_freq_frame_offsets))] < tgt_frequency_between_bursts):
         ax.set_title("No PAC, low frequency too low. PAC score: %2.2f" % (scores[ax_idx],))
      elif (frequencies_between_bursts[int(ax_idx % len(high_freq_frame_offsets))] > tgt_frequency_between_bursts):
         ax.set_title("No PAC, low frequency too high. PAC score: %2.2f" % (scores[ax_idx],))
      else:
         ax.set_title("Good PAC, low frequency matches. PAC score: %2.2f" % (scores[ax_idx],))

   def main():
      #Configure sample data
      data_range = np.arange(0, 10000)
      frequency_sampling = 1000
      frequencies_between_bursts = [5, 10, 15]
      tgt_frequency_between_bursts = 10
      frequency_within_bursts = 200
      high_freq_frame_offsets = [0, 20, 40]
      
      #Configure noise data
      random_noise_strength = 1
      
      #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 = list(); best_fits = list(); amplitude_signals = list()
      for high_freq_signal in high_freq_signals:
      for low_freq_signal in low_freq_signals:
         tmp = pac.run_dmi(low_freq_signal, high_freq_signal, phase_window_half_size = 20, phase_step_width = 5)
         scores.append(tmp[0]); best_fits.append(tmp[1]); amplitude_signals.append(tmp[2])
      
      #visualization
      visualize(scores, best_fits, amplitude_signals, frequencies_between_bursts, tgt_frequency_between_bursts, high_freq_frame_offsets)
      
      plt.tight_layout()
      plt.legend()
      figManager = plt.get_current_fig_manager()
      figManager.window.showMaximized()
      plt.show(block = True)

   main()

Using the direct modulation index, PAC has been detected only where there was a match between the high frequency component and the low frequency component:
    
.. image:: img/dmi.png