Surface NMR processing and inversion GUI
選択できるのは25トピックまでです。 トピックは、先頭が英数字で、英数字とダッシュ('-')を使用した35文字以内のものにしてください。

decay.py 20KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550
  1. import numpy, pylab,array #,rpy2
  2. from rpy2.robjects.packages import importr
  3. import rpy2.robjects as robjects
  4. import rpy2.robjects.numpy2ri
  5. #import notch
  6. from numpy.fft import fft, fftfreq
  7. # We know/can calculate frequency peak, use this to guess where picks will be.
  8. # maybe have a sliding window that reports peak values.
  9. def peakPicker(data, omega, dt):
  10. # compute window based on omega and dt
  11. # make sure you are not aliased, grab every other peak
  12. window = (2*numpy.pi) / (omega*dt)
  13. data = numpy.array(data)
  14. peaks = []
  15. troughs = []
  16. times = []
  17. times2 = []
  18. indices = []
  19. ws = 0
  20. we = window
  21. ii = 0
  22. for i in range((int)(len(data)/window)):
  23. # initially was just returning this I think avg is better
  24. #times.append( (ws + numpy.abs(data[ws:we]).argmax()) * dt )
  25. peaks.append(numpy.max(data[ws:we]))
  26. times.append( (ws + data[ws:we].argmax()) * dt )
  27. indices.append( ii + data[ws:we].argmax() )
  28. troughs.append(numpy.min(data[ws:we]))
  29. times2.append( (ws + (data[ws:we]).argmin()) * dt )
  30. indices.append( ii + data[ws:we].argmin() )
  31. ws += window
  32. we += window
  33. ii += (int)(we-ws)
  34. #return numpy.array(peaks), numpy.array(times)
  35. # Averaging peaks does a good job of removing bias in noise
  36. return (numpy.array(peaks)-numpy.array(troughs))/2., \
  37. (numpy.array(times)+numpy.array(times2))/2., \
  38. indices
  39. #################################################
  40. # Regress for T2 using rpy2 interface
  41. def regressCurve(peaks,times,sigma2=None ,intercept=True):
  42. # TODO, if regression fails, it might be because there is no exponential
  43. # term, maybe do a second regression then on a linear model.
  44. b1 = 0 # Bias
  45. b2 = 0 # Linear
  46. rT2 = 0.3 # T2 regressed
  47. r = robjects.r
  48. # Variable shared between R and Python
  49. robjects.globalenv['b1'] = b1
  50. robjects.globalenv['b2'] = b2
  51. robjects.globalenv['rT2'] = rT2
  52. #robjects.globalenv['sigma2'] = sigma2
  53. value = robjects.FloatVector(peaks)
  54. times = robjects.FloatVector(numpy.array(times))
  55. # my_weights = robjects.RVector(value/sigma2)
  56. # robjects.globalenv['my_weigts'] = my_weights
  57. if sigma2 != None:
  58. # print ("weighting")
  59. #tw = numpy.array(peaks)/sigma2
  60. my_weights = robjects.FloatVector( sigma2 )
  61. #else:
  62. # my_weights = robjects.FloatVector(numpy.ones(len(peaks)))
  63. robjects.globalenv['my_weights'] = my_weights
  64. #print (my_weights)
  65. #print (len(peaks))
  66. if (intercept):
  67. my_list = robjects.r('list(b1=50, b2=1e2, rT2=0.03)')
  68. my_lower = robjects.r('list(b1=0, b2=0, rT2=.005)')
  69. my_upper = robjects.r('list(b1=20000, b2=2000, rT2=.700)')
  70. else:
  71. my_list = robjects.r('list(b2=1e2, rT2=0.3)')
  72. my_lower = robjects.r('list(b2=0, rT2=.005)')
  73. my_upper = robjects.r('list(b2=2000, rT2=.700)')
  74. my_cont = robjects.r('nls.control(maxiter=1000, warnOnly=TRUE, printEval=FALSE)')
  75. if (intercept):
  76. #fmla = robjects.RFormula('value ~ b1 + exp(-times/rT2)')
  77. fmla = robjects.Formula('value ~ b1 + b2*exp(-times/rT2)')
  78. #fmla = robjects.RFormula('value ~ b1 + b2*times + exp(-times/rT2)')
  79. else:
  80. fmla = robjects.Formula('value ~ b2*exp(-times/rT2)')
  81. env = fmla.getenvironment()
  82. env['value'] = value
  83. env['times'] = times
  84. # ugly, but I get errors with everything else I've tried
  85. #my_weights = robjects.r('rep(1,length(value))')
  86. #for ii in range(len(my_weights)):
  87. # my_weights[ii] *= peaks[ii]/sigma2
  88. Error = False
  89. #fit = robjects.r.nls(fmla,start=my_list,control=my_cont,weights=my_weights)
  90. if (sigma2 != None):
  91. #print("SIGMA 2")
  92. #fit = robjects.r.tryCatch(robjects.r.suppressWarnings(robjects.r.nls(fmla,start=my_list,control=my_cont,algorithm="port", \
  93. # weights=my_weights)), 'silent=TRUE')
  94. try:
  95. fit = robjects.r.tryCatch( robjects.r.nls(fmla, start=my_list, control=my_cont, weights=my_weights, algorithm="port" , \
  96. lower=my_lower,upper=my_upper))
  97. except:
  98. print("regression issue pass")
  99. Error = True
  100. # weights=my_weights))
  101. else:
  102. try:
  103. fit = robjects.r.tryCatch(robjects.r.nls(fmla,start=my_list,control=my_cont,algorithm="port",lower=my_lower,upper=my_upper))
  104. except:
  105. print("regression issue pass")
  106. Error = True
  107. # If failure fall back on zero regression values
  108. if not Error:
  109. #Error = fit[3][0]
  110. report = r.summary(fit)
  111. b1 = 0
  112. b2 = 0
  113. rT2 = 1
  114. if (intercept):
  115. if not Error:
  116. b1 = r['$'](report,'par')[0]
  117. b2 = r['$'](report,'par')[1]
  118. rT2 = r['$'](report,'par')[2]
  119. #print report
  120. #print r['$'](report,'convergence')
  121. #print r['convergence'] #(report,'convergence')
  122. #print r['$'](report,'par')[13]
  123. #print r['$'](report,'par')[14]
  124. else:
  125. print("ERROR DETECTED, regressed values set to default")
  126. b1 = 1e1
  127. b2 = 1e-2
  128. rT2 = 1e-2
  129. #print r['$'](report,'par')[0]
  130. #print r['$'](report,'par')[1]
  131. #print r['$'](report,'par')[2]
  132. return [b1,b2,rT2]
  133. else:
  134. if not Error:
  135. rT2 = r['$'](report,'par')[1]
  136. b2 = r['$'](report,'par')[0]
  137. else:
  138. print("ERROR DETECTED, regressed values set to default")
  139. return [b2, rT2]
  140. def quadratureDetect(X, Y, tt):
  141. r = robjects.r
  142. robjects.r('''
  143. Xc <- function(E0, df, tt, phi, T2) {
  144. E0 * -sin(2*pi*df*tt + phi) * exp(-tt/T2)
  145. }
  146. Yc <- function(E0, df, tt, phi, T2) {
  147. E0 * cos(2*pi*df*tt + phi) * exp(-tt/T2)
  148. }
  149. ''')
  150. # Make 0 vector
  151. Zero = robjects.FloatVector(numpy.zeros(len(X)))
  152. # Fitted Parameters
  153. E0 = 0.
  154. df = 0.
  155. phi = 0.
  156. T2 = 0.
  157. robjects.globalenv['E0'] = E0
  158. robjects.globalenv['df'] = df
  159. robjects.globalenv['phi'] = phi
  160. robjects.globalenv['T2'] = T2
  161. XY = robjects.FloatVector(numpy.concatenate((X,Y)))
  162. # Arrays
  163. tt = robjects.FloatVector(numpy.array(tt))
  164. X = robjects.FloatVector(numpy.array(X))
  165. Y = robjects.FloatVector(numpy.array(Y))
  166. Zero = robjects.FloatVector(numpy.array(Zero))
  167. #fmla = robjects.Formula('Zero ~ QI( E0, df, tt, phi, T2, X, Y )')
  168. #fmla = robjects.Formula('X ~ Xc( E0, df, tt, phi, T2 )')
  169. #fmla = robjects.Formula('Y ~ Yc( E0, df, tt, phi, T2 )')
  170. fmla = robjects.Formula('XY ~ c(Xc( E0, df, tt, phi, T2 ), Yc( E0, df, tt, phi, T2 ))')
  171. env = fmla.getenvironment()
  172. env['Zero'] = Zero
  173. env['X'] = X
  174. env['Y'] = Y
  175. env['XY'] = XY
  176. env['tt'] = tt
  177. # Bounds and control
  178. start = robjects.r('list(E0=100, df= 0 , phi= 0.00, T2=.100)')
  179. lower = robjects.r('list(E0=1, df=-13.0, phi= -3.14, T2=.005)')
  180. upper = robjects.r('list(E0=1000, df= 13.0, phi= 3.14, T2=.800)')
  181. cont = robjects.r('nls.control(maxiter=10000, warnOnly=TRUE, printEval=FALSE)')
  182. fit = robjects.r.tryCatch(robjects.r.nls(fmla, start=start, control=cont, lower=lower, upper=upper, algorithm='port')) #, \
  183. report = r.summary(fit)
  184. #print (report)
  185. E0 = r['$'](report,'par')[0]
  186. df = r['$'](report,'par')[1]
  187. phi = r['$'](report,'par')[2]
  188. T2 = r['$'](report,'par')[3]
  189. #print ( E0,df,phi,T2 )
  190. return E0,df,phi,T2
  191. #################################################
  192. # Regress for T2 using rpy2 interface
  193. def regressSpec(w, wL, X): #,sigma2=1,intercept=True):
  194. # compute s
  195. s = -1j*w
  196. # TODO, if regression fails, it might be because there is no exponential
  197. # term, maybe do a second regression then on a linear model.
  198. a = 0 # Linear
  199. rT2 = 0.1 # T2 regressed
  200. r = robjects.r
  201. # Variable shared between R and Python
  202. robjects.globalenv['a'] = a
  203. #robjects.globalenv['w'] = w
  204. robjects.globalenv['rT2'] = rT2
  205. robjects.globalenv['wL'] = wL
  206. robjects.globalenv['nb'] = 0
  207. #s = robjects.ComplexVector(numpy.array(s))
  208. w = robjects.FloatVector(numpy.array(w))
  209. XX = robjects.FloatVector(X)
  210. #Xr = robjects.FloatVector(numpy.real(X))
  211. #Xi = robjects.FloatVector(numpy.imag(X))
  212. #Xa = robjects.FloatVector(numpy.abs(X))
  213. #Xri = robjects.FloatVector(numpy.concatenate((Xr,Xi)))
  214. #my_lower = robjects.r('list(a=.001, rT2=.001, nb=.0001)')
  215. my_lower = robjects.r('list(a=.001, rT2=.001)')
  216. #my_upper = robjects.r('list(a=1.5, rT2=.300, nb =100.)')
  217. my_upper = robjects.r('list(a=1.5, rT2=.300)')
  218. #my_list = robjects.r('list(a=.2, rT2=0.03, nb=.1)')
  219. my_list = robjects.r('list(a=.2, rT2=0.03)')
  220. my_cont = robjects.r('nls.control(maxiter=5000, warnOnly=TRUE, printEval=FALSE)')
  221. #fmla = robjects.Formula('Xri ~ c(a*Re((wL) / (wL^2+(s+1/rT2)^2 )), a*Im((wL)/(wL^2 + (s+1/rT2)^2 )))') # envelope
  222. ##fmla = robjects.Formula('Xri ~ c(a*Re((wL) / (wL^2+(s+1/rT2)^2 )), a*Im((wL)/(wL^2 + (s+1/rT2)^2 )))') # envelope
  223. #fmla = robjects.Formula('XX ~ a*(wL) / (wL^2 + (s+1/rT2)^2 )') # complex
  224. #fmla = robjects.Formula('Xa ~ abs(a*(wL) / (wL^2 + (s+1/rT2)^2 )) + nb') # complex
  225. #fmla = robjects.Formula('XX ~ Re(a*( s + 1./rT2 ) / (wL^2 + (s+1/rT2)^2 ))') # complex
  226. fmla = robjects.Formula('XX ~ a*(.5/rT2) / ((1./rT2)^2 + (w-wL)^2 )')
  227. #fmla = robjects.Formula('Xa ~ (s + 1./T2) / ( wL**2 + (1/T2 + 1j*w)**2 ) ')
  228. env = fmla.getenvironment()
  229. #env['s'] = s
  230. env['w'] = w
  231. #env['Xr'] = Xr
  232. #env['Xa'] = Xa
  233. #env['Xi'] = Xi
  234. #env['Xri'] = Xri
  235. env['XX'] = XX
  236. #fit = robjects.r.tryCatch(robjects.r.nls(fmla,start=my_list, control=my_cont)) #, lower=my_lower, algorithm='port')) #, \
  237. fit = robjects.r.tryCatch(robjects.r.nls(fmla, start=my_list, control=my_cont, lower=my_lower, upper=my_upper, algorithm='port')) #, \
  238. report = r.summary(fit)
  239. #print report
  240. #print(r.warnings())
  241. a = r['$'](report,'par')[0]
  242. rT2 = r['$'](report,'par')[1]
  243. nb = r['$'](report,'par')[2]
  244. return a, rT2, nb
  245. #################################################
  246. # Regress for T2 using rpy2 interface
  247. def regressModulus(w, wL, X): #,sigma2=1,intercept=True):
  248. # compute s
  249. s = -1j*w
  250. # TODO, if regression fails, it might be because there is no exponential
  251. # term, maybe do a second regression then on a linear model.
  252. a = 0 # Linear
  253. rT2 = 0.1 # T2 regressed
  254. r = robjects.r
  255. # Variable shared between R and Python
  256. robjects.globalenv['a'] = a
  257. robjects.globalenv['rT2'] = rT2
  258. robjects.globalenv['wL'] = wL
  259. robjects.globalenv['nb'] = 0
  260. s = robjects.ComplexVector(numpy.array(s))
  261. XX = robjects.ComplexVector(X)
  262. Xr = robjects.FloatVector(numpy.real(X))
  263. Xi = robjects.FloatVector(numpy.imag(X))
  264. Xa = robjects.FloatVector(numpy.abs(X))
  265. Xri = robjects.FloatVector(numpy.concatenate((Xr,Xi)))
  266. #my_lower = robjects.r('list(a=.001, rT2=.001, nb=.0001)')
  267. my_lower = robjects.r('list(a=.001, rT2=.001)')
  268. #my_upper = robjects.r('list(a=1.5, rT2=.300, nb =100.)')
  269. my_upper = robjects.r('list(a=1.5, rT2=.300)')
  270. #my_list = robjects.r('list(a=.2, rT2=0.03, nb=.1)')
  271. my_list = robjects.r('list(a=.2, rT2=0.03)')
  272. my_cont = robjects.r('nls.control(maxiter=5000, warnOnly=TRUE, printEval=FALSE)')
  273. #fmla = robjects.Formula('Xri ~ c(a*Re((wL) / (wL^2+(s+1/rT2)^2 )), a*Im((wL)/(wL^2 + (s+1/rT2)^2 )))') # envelope
  274. ##fmla = robjects.Formula('Xri ~ c(a*Re((wL) / (wL^2+(s+1/rT2)^2 )), a*Im((wL)/(wL^2 + (s+1/rT2)^2 )))') # envelope
  275. #fmla = robjects.Formula('XX ~ a*(wL) / (wL^2 + (s+1/rT2)^2 )') # complex
  276. #fmla = robjects.Formula('Xa ~ abs(a*(wL) / (wL^2 + (s+1/rT2)^2 )) + nb') # complex
  277. fmla = robjects.Formula('Xa ~ abs(a*(wL) / (wL^2 + (s+1/rT2)^2 ))') # complex
  278. env = fmla.getenvironment()
  279. env['s'] = s
  280. env['Xr'] = Xr
  281. env['Xa'] = Xa
  282. env['Xi'] = Xi
  283. env['Xri'] = Xri
  284. env['XX'] = XX
  285. #fit = robjects.r.tryCatch(robjects.r.nls(fmla,start=my_list, control=my_cont)) #, lower=my_lower, algorithm='port')) #, \
  286. fit = robjects.r.tryCatch(robjects.r.nls(fmla, start=my_list, control=my_cont, lower=my_lower, upper=my_upper, algorithm='port')) #, \
  287. report = r.summary(fit)
  288. #print report
  289. #print r.warnings()
  290. a = r['$'](report,'par')[0]
  291. rT2 = r['$'](report,'par')[1]
  292. nb = r['$'](report,'par')[2]
  293. return a, rT2
  294. #################################################
  295. # Regress for T2 using rpy2 interface
  296. def regressSpecComplex(w, wL, X, known=True, win=None): #,sigma2=1,intercept=True):
  297. # compute s
  298. s = -1j*w
  299. # TODO, if regression fails, it might be because there is no exponential
  300. # term, maybe do a second regression then on a linear model.
  301. a = 1 # Linear
  302. rT2 = 0.1 # T2 regressed
  303. r = robjects.r
  304. phi2 = 0 # phase
  305. wL2 = wL
  306. # Variable shared between R and Python
  307. robjects.globalenv['a'] = a
  308. robjects.globalenv['rT2'] = rT2
  309. robjects.globalenv['wL'] = wL
  310. robjects.globalenv['wL2'] = 0
  311. robjects.globalenv['nb'] = 0
  312. robjects.globalenv['phi2'] = phi2
  313. s = robjects.ComplexVector(numpy.array(s))
  314. XX = robjects.ComplexVector(X)
  315. Xr = robjects.FloatVector(numpy.real(X))
  316. Xi = robjects.FloatVector(numpy.imag(X))
  317. Xa = robjects.FloatVector(numpy.abs(X))
  318. Xri = robjects.FloatVector(numpy.concatenate((X.real,X.imag)))
  319. robjects.r('''
  320. source('kernel.r')
  321. ''')
  322. #Kw = robjects.globalenv['Kwri']
  323. #print (numpy.shape(X))
  324. #######################################################################
  325. if known:
  326. # known Frequency
  327. my_lower = robjects.r('list(a=.001, rT2=.001, phi2=-3.14)')
  328. my_upper = robjects.r('list(a=3.5, rT2=.300, phi2=3.14)')
  329. my_list = robjects.r('list(a=.2, rT2=0.03, phi2=0)')
  330. else:
  331. # Unknown Frequency
  332. my_lower = robjects.r('list(a=.001, rT2=.001, phi2=-3.14, wL2=wL-5)')
  333. my_upper = robjects.r('list(a=3.5, rT2=.300, phi2=3.14, wL2=wL+5)')
  334. my_list = robjects.r('list(a=.2, rT2=0.03, phi2=0, wL2=wL)')
  335. my_cont = robjects.r('nls.control(maxiter=5000, warnOnly=TRUE, printEval=FALSE)')
  336. #fmla = robjects.Formula('Xri ~ c(a*Re((wL) / (wL^2+(s+1/rT2)^2 )), a*Im((wL)/(wL^2 + (s+1/rT2)^2 )))') # envelope
  337. #fmla = robjects.Formula('Xi ~ Im(a*(sin(phi2)*s + ((1/rT2)*sin(phi2)) + wL*cos(phi2)) / (wL^2+(s+1/rT2)^2 ))') # envelope
  338. #fmla = robjects.Formula('Xri ~ c(Re(a*(sin(phi2)*s + ((1/rT2)*sin(phi2)) + wL*cos(phi2)) / (wL^2+(s+1/rT2)^2 )), Im(a*(sin(phi2)*s + ((1/rT2)*sin(phi2)) + wL*cos(phi2)) / (wL^2+(s+1/rT2)^2 )))') # envelope
  339. #fmlar = robjects.Formula('Xr ~ (Kwr(a, phi2, s, rT2, wL)) ') # envelope
  340. #fmlai = robjects.Formula('Xi ~ (Kwi(a, phi2, s, rT2, wL)) ') # envelope
  341. if known:
  342. ###########################################3
  343. # KNOWN freq
  344. fmla = robjects.Formula('Xri ~ c(Kwr(a, phi2, s, rT2, wL), Kwi(a, phi2, s, rT2, wL) ) ') # envelope
  345. else:
  346. ####################################################################################################3
  347. # unknown frequency
  348. fmla = robjects.Formula('Xri ~ c(Kwr(a, phi2, s, rT2, wL2), Kwi(a, phi2, s, rT2, wL2) ) ') # envelope
  349. #fmla = robjects.Formula('Xri ~ (Kwri(a, phi2, s, rT2, wL)) ') # envelope
  350. #fmla = robjects.Formula('Xa ~ (abs(a*(sin(phi2)*s + ((1/rT2)*sin(phi2)) + wL*cos(phi2)) / (wL^2+(s+1/rT2)^2 )))') # envelope
  351. #fmla = robjects.Formula('XX ~ a*(wL) / (wL^2 + (s+1/rT2)^2 )') # complex
  352. #fmla = robjects.Formula('Xa ~ abs(a*(wL) / (wL^2 + (s+1/rT2)^2 )) + nb') # complex
  353. #fmla = robjects.Formula('Xri ~ c(a*Re((wL) / (wL^2+(s+1/rT2)^2 )), a*Im((wL)/(wL^2 + (s+1/rT2)^2 )))') # envelope
  354. # self.Gw[iw, iT2] = ((np.sin(phi2) * (alpha + 1j*self.w[iw]) + self.wL*np.cos(phi2)) / \
  355. # (self.wL**2 + (alpha+1.j*self.w[iw])**2 ))
  356. # self.Gw[iw, iT2] = ds * self.sc*((np.sin(phi2)*( alpha + 1j*self.w[iw]) + self.wL*np.cos(phi2)) / \
  357. # (self.wL**2 + (alpha+1.j*self.w[iw])**2 ))
  358. # Works Amplitude Only!
  359. #fmla = robjects.Formula('Xa ~ abs(a*(wL) / (wL^2 + (s+1/rT2)^2 ))') # complex
  360. env = fmla.getenvironment()
  361. env['s'] = s
  362. env['Xr'] = Xr
  363. env['Xa'] = Xa
  364. env['Xi'] = Xi
  365. env['Xri'] = Xri
  366. env['XX'] = XX
  367. #fit = robjects.r.tryCatch(robjects.r.nls(fmla,start=my_list, control=my_cont)) #, lower=my_lower, algorithm='port')) #, \
  368. #fit = robjects.r.tryCatch(robjects.r.nls(fmlar, start=my_list, control=my_cont, lower=my_lower, upper=my_upper, algorithm='port')) #, \
  369. fit = robjects.r.tryCatch(robjects.r.nls(fmla, start=my_list, control=my_cont, lower=my_lower, upper=my_upper, algorithm='port')) #, \
  370. #env = fmlai.getenvironment()
  371. #fiti = robjects.r.tryCatch(robjects.r.nls(fmlai, start=my_list, control=my_cont, lower=my_lower, upper=my_upper, algorithm='port')) #, \
  372. #reportr = r.summary(fitr)
  373. #reporti = r.summary(fiti)
  374. report = r.summary(fit)
  375. #print( report )
  376. #exit()
  377. #print( reportr )
  378. #print( reporti )
  379. #exit()
  380. #print ( r.warnings())
  381. #a = (r['$'](reportr,'par')[0] + r['$'](reporti,'par')[0]) / 2.
  382. #rT2 = (r['$'](reportr,'par')[1] + r['$'](reporti,'par')[1]) / 2.
  383. #nb = (r['$'](reportr,'par')[2] + r['$'](reporti,'par')[2]) / 2.
  384. a = r['$'](report,'par')[0]
  385. rT2 = r['$'](report,'par')[1]
  386. nb = r['$'](report,'par')[2] #phi2
  387. #print ("Python wL2", r['$'](report,'par')[3] )
  388. #print ("Python zeta", r['$'](report,'par')[2] )
  389. return a, rT2, nb
  390. ###################################################################
  391. ###################################################################
  392. ###################################################################
  393. if __name__ == "__main__":
  394. dt = .0001
  395. T2 = .1
  396. omega = 2000.*2*numpy.pi
  397. phi = .0
  398. T = 8.*T2
  399. t = numpy.arange(0, T, dt)
  400. # Synthetic data, simple single decaying sinusoid
  401. # with a single decay parameter and gaussian noise added
  402. data = numpy.exp(-t/T2) * numpy.sin(omega * t + phi) + numpy.random.normal(0,.05,len(t)) \
  403. + numpy.random.randint(-1,2,len(t))*numpy.random.exponential(.2,len(t))
  404. cdata = numpy.exp(-t/T2) * numpy.sin(omega * t + phi) #+ numpy.random.normal(0,.25,len(t))
  405. #data = numpy.random.normal(0,.25,len(t))
  406. sigma2 = numpy.std(data[::-len(data)/4])
  407. #sigma2 = numpy.var(data[::-len(data)/4])
  408. print("sigma2", sigma2)
  409. [peaks,times,indices] = peakPicker(data, omega, dt)
  410. [b1,b2,rT2] = regressCurve(peaks,times)
  411. print("rT2 nonweighted", rT2)
  412. [b1,b2,rT2] = regressCurve(peaks,times,sigma2)
  413. print("rT2 weighted", rT2)
  414. envelope = numpy.exp(-t/T2)
  415. renvelope = numpy.exp(-t/rT2)
  416. outf = file('regress.txt','w')
  417. for i in range(len(times)):
  418. outf.write(str(times[i]) + " " + str(peaks[i]) + "\n")
  419. outf.close()
  420. pylab.plot(t,data, 'b')
  421. pylab.plot(t,cdata, 'g', linewidth=1)
  422. pylab.plot(t,envelope, color='violet', linewidth=4)
  423. pylab.plot(t,renvelope, 'r', linewidth=4)
  424. pylab.plot(times, numpy.array(peaks), 'bo', markersize=8, alpha=.25)
  425. pylab.legend(['noisy data','clean data','real envelope','regressed env','picks'])
  426. pylab.savefig("regression.pdf")
  427. # FFT check
  428. fourier = fft(data)
  429. pylab.figure()
  430. freq = fftfreq(len(data), d=dt)
  431. pylab.plot(freq, (fourier.real))
  432. pylab.show()
  433. # TODO do a bunch in batch mode to see if T2 estimate is better with or without
  434. # weighting and which model is best.
  435. # TODO try with real data
  436. # TODO test filters (median, FFT, notch)
  437. # It looks like weighting is good for relatively low sigma, but for noisy data
  438. # it hurts us. Check