In [1]:

    

import sys
sys.path.insert(0, '..')
import time
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec
%matplotlib notebook 
import numpy as np
from scipy import signal
from scipy import stats

from Configuration import Configuration
from MotorUnitPool import MotorUnitPool


    

In [2]:

    

conf = Configuration('confEMGCoherenceTA.rmto')
conf.simDuration_ms = 5000 # Here I change simulation duration without changing the Configuration file.


    

In [3]:

    

t = np.arange(0.0, conf.simDuration_ms, conf.timeStep_ms)


    

In [4]:

    

pools = dict()
pools[0] = MotorUnitPool(conf, 'SOL')
pools[1] = MotorUnitPool(conf, 'SOL')


    

    




Muscle spindle from muscle SOL built.
Motor Unit Pool SOL built
Muscle spindle from muscle SOL built.
Motor Unit Pool SOL built

In [5]:

    

tic = time.clock()
for i in xrange(0, len(t)-1):
    current1 = 4 + np.random.randn(1)
    current2 = 4 + np.random.randn(1)
    for j in xrange(len(pools[0].unit)):
        pools[0].iInjected[2*j+1] = (current1 + 0.7*current2)/1.2
    for j in xrange(len(pools[1].unit)):
        pools[1].iInjected[2*j+1] = (0.7*current1 + current2)/1.2
    pools[0].atualizeMotorUnitPool(t[i])
    pools[1].atualizeMotorUnitPool(t[i])
toc = time.clock()
print str(toc - tic) + ' seconds'


    

    




2415.655879 seconds

In [6]:

    

pools[0].listSpikes()
pools[1].listSpikes()


    

In [7]:

    

plt.figure()
plt.plot(pools[0].poolTerminalSpikes[:, 0],
         pools[0].poolTerminalSpikes[:, 1] + 1, '.')
plt.xlabel('t (ms)')
plt.ylabel('Motor Unit index')


    

    















    












    
Out[7]:





<matplotlib.text.Text at 0x7f6bde5498d0>

In [8]:

    

plt.figure()
plt.plot(pools[1].poolTerminalSpikes[:, 0],
         pools[1].poolTerminalSpikes[:, 1] + 1, '.')
plt.xlabel('t (ms)')
plt.ylabel('Motor Unit index')


    

    















    












    
Out[8]:





<matplotlib.text.Text at 0x7f6bde461d10>

In [9]:

    

plt.figure()
plt.plot(t, pools[0].Muscle.force, '-')
plt.xlabel('t (ms)')
plt.ylabel('Muscle Force (N)')


    

    















    












    
Out[9]:





<matplotlib.text.Text at 0x7f6bde3f8810>

In [10]:

    

plt.figure()
plt.plot(t, pools[1].Muscle.force, '-')
plt.xlabel('t (ms)')
plt.ylabel('Muscle Force (N)')


    

    















    












    
Out[10]:





<matplotlib.text.Text at 0x7f6bde3800d0>

In [11]:

    

pools[0].getMotorUnitPoolEMG()


    

In [12]:

    

plt.figure()
plt.plot(t, pools[0].emg, '-')
plt.xlabel('t (ms)')
plt.ylabel('EMG (mV)')


    

    















    












    
Out[12]:





<matplotlib.text.Text at 0x7f6bde182e90>

In [13]:

    

pools[1].getMotorUnitPoolEMG()


    

In [14]:

    

plt.figure()
plt.plot(t, pools[1].emg, '-')
plt.xlabel('t (ms)')
plt.ylabel('EMG (mV)')


    

    















    












    
Out[14]:





<matplotlib.text.Text at 0x7f6bde4700d0>

In [15]:

    

window = 2000
fs = 1/0.00005
f, coherence = signal.coherence(signal.detrend(pools[0].emg[int(round(1000/0.05)):-1]), 
                                signal.detrend(pools[1].emg[int(round(1000/0.05)):-1]), 
                                fs=fs, window='hann', nperseg=window, noverlap=None, nfft=window, detrend='constant', axis=0)
divisions = len(pools[0].emg[int(round(1000/0.05)):-1])/window;
erroCoh = stats.f.ppf(0.95,2, 2*(divisions-1))/(divisions - 1 + stats.f.ppf(0.95,2,2*(divisions-1)));


    

In [16]:

    

plt.figure()
plt.plot(f, coherence,'-', f, erroCoh*np.ones(len(f)),'--')
plt.xlabel('f (Hz)')
plt.ylabel('EMG (mV)')
plt.xlim((0,500))


    

    















    












    
Out[16]:





(0, 500)

In [ ]: