Guest User

Untitled

a guest
Nov 21st, 2017
91
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
text 1.41 KB | None | 0 0
  1. # to read the list of stcs for all the epochs
  2. with open('/home/daniel/Dropbox/F[...]', 'rb') as f:
  3. label_ts = pickle.load(f)
  4.  
  5. x = np.asarray(label_ts)
  6. nfft = 512
  7. n_freqs = nfft/2+1
  8. n_epochs = len(x) # in this case there are 120 epochs
  9. channels = 68
  10. sfreq = 1017.25
  11.  
  12. def compute_mean_psd_csd(x, n_epochs, nfft, sfreq):
  13. '''Computes mean of PSD and CSD for signals.'''
  14.  
  15. Rxy = np.zeros((n_epochs, channels, channels, n_freqs), dtype=complex)
  16. Rxx = np.zeros((n_epochs, channels, channels, n_freqs))
  17. Ryy = np.zeros((n_epochs, channels, channels, n_freqs))
  18. for i in xrange(0, n_epochs):
  19. print('computing connectivity for epoch %s'%(i+1))
  20. for j in xrange(0, channels):
  21. for k in xrange(0, channels):
  22. Rxy[i,j,k], freqs = mlab.csd(x[j], x[k], NFFT=nfft, Fs=sfreq)
  23. Rxx[i,j,k], _____ = mlab.psd(x[j], NFFT=nfft, Fs=sfreq)
  24. Ryy[i,j,k], _____ = mlab.psd(x[k], NFFT=nfft, Fs=sfreq)
  25.  
  26. Rxy_mean = np.mean(Rxy, axis=0, dtype=np.float32)
  27. Rxx_mean = np.mean(Rxx, axis=0, dtype=np.float32)
  28. Ryy_mean = np.mean(Ryy, axis=0, dtype=np.float32)
  29.  
  30. return freqs, Rxy, Rxy_mean, np.real(Rxx_mean), np.real(Ryy_mean)
  31.  
  32. if x[j] not in cache_psd:
  33. cache_psd[x[j]], ____ = mlab.psd(x[j], NFFT=nfft, Fs=sfreq)
  34. Rxx[i,j,k] = cache_psd[x[j]]
  35.  
  36. if x[k] not in cache_psd:
  37. cache_psd[x[k]], ____ = mlab.psd(x[k], NFFT=nfft, Fs=sfreq)
  38. Ryy[i,j,k] = cache_psd[x[k]]
Add Comment
Please, Sign In to add comment