Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- # to read the list of stcs for all the epochs
- with open('/home/daniel/Dropbox/F[...]', 'rb') as f:
- label_ts = pickle.load(f)
- x = np.asarray(label_ts)
- nfft = 512
- n_freqs = nfft/2+1
- n_epochs = len(x) # in this case there are 120 epochs
- channels = 68
- sfreq = 1017.25
- def compute_mean_psd_csd(x, n_epochs, nfft, sfreq):
- '''Computes mean of PSD and CSD for signals.'''
- Rxy = np.zeros((n_epochs, channels, channels, n_freqs), dtype=complex)
- Rxx = np.zeros((n_epochs, channels, channels, n_freqs))
- Ryy = np.zeros((n_epochs, channels, channels, n_freqs))
- for i in xrange(0, n_epochs):
- print('computing connectivity for epoch %s'%(i+1))
- for j in xrange(0, channels):
- for k in xrange(0, channels):
- Rxy[i,j,k], freqs = mlab.csd(x[j], x[k], NFFT=nfft, Fs=sfreq)
- Rxx[i,j,k], _____ = mlab.psd(x[j], NFFT=nfft, Fs=sfreq)
- Ryy[i,j,k], _____ = mlab.psd(x[k], NFFT=nfft, Fs=sfreq)
- Rxy_mean = np.mean(Rxy, axis=0, dtype=np.float32)
- Rxx_mean = np.mean(Rxx, axis=0, dtype=np.float32)
- Ryy_mean = np.mean(Ryy, axis=0, dtype=np.float32)
- return freqs, Rxy, Rxy_mean, np.real(Rxx_mean), np.real(Ryy_mean)
- if x[j] not in cache_psd:
- cache_psd[x[j]], ____ = mlab.psd(x[j], NFFT=nfft, Fs=sfreq)
- Rxx[i,j,k] = cache_psd[x[j]]
- if x[k] not in cache_psd:
- cache_psd[x[k]], ____ = mlab.psd(x[k], NFFT=nfft, Fs=sfreq)
- Ryy[i,j,k] = cache_psd[x[k]]
Add Comment
Please, Sign In to add comment