-
Notifications
You must be signed in to change notification settings - Fork 0
/
FPPwrapper.py
44 lines (39 loc) · 1.83 KB
/
FPPwrapper.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
# wrapper for an object that opens a patient file finds the power specturm
# then runs a simulation given C and rate and finds the returns the residues
# of the difference of the PSDs
#
# Kristian Weegink, Nov 2013
import numpy as np
from scipy import signal
from matplotlib import pylab
import FPP
import BrainSim as BS
import math
import scipy.io.wavfile as wav
import Logger as lg
class Patient:
# intialize by finding the patient PSD
def __init__(self,patient_file,logger=lg.logger(0)):
self.log = logger
self.sr,self.Vt = wav.read(patient_file)
self.log.info('patient loaded')
self.nfft=2**int(math.log(len(self.Vt),2))+1
self.Pxx,self.freqs=pylab.psd(x=self.Vt,Fs=self.sr,NFFT=self.nfft,window=pylab.window_none, noverlap=1)
self.bartlet = self.Pxx
self.log.info('patient initialized')
# find the L2 difference between the FPP simulation and patient PSDs
def fpp(self,rate,C):
self.log.info('C = ' + str(C) + '; rate = ' + str(rate) )
Vt,t = FPP.FPP(log=lg.logger(0), N=10000, dt=1./24000, distributionParameter=[C,rate], plotAll=False,efield=False)
Pxi,freqs = pylab.psd(x=Vt,Fs=self.sr,NFFT=self.nfft,window=pylab.window_none, noverlap=1)
self.log.info('L2 = ' + str(np.sqrt((np.square(np.subtract(self.Pxx[370:3300],10000*Pxi[370:3300])).sum()))))
#pylab.show()
#pylab.loglog(10000*Pxi[37:3300])
#pylab.loglog(self.Pxx[37:3300])
#pylab.show()
#pylab.plot(self.freqs)
#pylab.show()
return np.subtract(self.Pxx[37:3300],10000*Pxi[37:3300])
def recalcBartlet(self,current):
curentPSD,freqs = pylab.psd(x=current,Fs=self.sr,NFFT=self.nfft,window=pylab.window_none, noverlap=1)
self.bartlet = np.divide(np.divide(self.Pxx,np.sum(self.Pxx)),np.divide(curentPSD,np.sum(curentPSD)))