Surface NMR processing and inversion GUI
Você não pode selecionar mais de 25 tópicos Os tópicos devem começar com uma letra ou um número, podem incluir traços ('-') e podem ter até 35 caracteres.

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263
  1. # #################################################################################
  2. # # GJI final pub specs #
  3. # import matplotlib #
  4. # from matplotlib import rc #
  5. # matplotlib.rcParams['text.latex.preamble']=[r"\usepackage{timet,amsmath}"] #
  6. # rc('font',**{'family':'serif','serif':['timet']}) #
  7. # rc('font',**{'size':8}) #
  8. # rc('text', usetex=True) #
  9. # # converts pc that GJI is defined in to inches #
  10. # # In GJI \textwidth = 42pc #
  11. # # \columnwidth = 20pc #
  12. # def pc2in(pc): #
  13. # return pc*12/72.27 #
  14. # #################################################################################
  15. #from GJIPlot import *
  16. import numpy as np
  17. import matplotlib.pyplot as plt
  18. #from invertColours import *
  19. from akvo.tressel.decay import *
  20. from scipy import signal
  21. def quadrature(T, vL, wL, dt, xn, DT, t):
  22. # decimate
  23. # blind decimation
  24. # 1 instead of T
  25. irsamp = (T) * int( (1./vL) / dt) # real
  26. iisamp = int( ((1./vL)/ dt) * ( .5*np.pi / (2.*np.pi) ) ) # imaginary
  27. dsamp = int( DT / dt) # real
  28. iisamp += dsamp
  29. ############################################################
  30. # simple quadrature-detection via sampling
  31. xr = xn[dsamp::irsamp]
  32. xi = xn[iisamp::irsamp]
  33. phase = np.angle( xr + 1j*xi )
  34. abse = np.abs( xr + 1j*xi )
  35. # times
  36. #ta = np.arange(0, TT, dt)
  37. #te = np.arange(DT, TT, TT/len(abse))
  38. #############################################################
  39. # hilbert transform
  40. ht = signal.hilbert(xn) #, 100))
  41. he = np.abs(ht) #, 100))
  42. hp = ((np.angle(ht[dsamp::irsamp])))
  43. #############################################################
  44. # Resample ht
  45. #htd = signal.decimate(he, 100, ftype='fir')
  46. #td = signal.decimate(t, 100, ftype='fir')
  47. #[htd, td] = signal.resample(he, 21, t)
  48. #toss first, and use every third
  49. #htd = htd[1::3]
  50. #td = td[1::3]
  51. #############################################################
  52. # Pre-envelope
  53. #gplus = xn + 1j*ht
  54. #############################################################
  55. # Complex envelope
  56. #gc = gplus / np.exp(1j*wL*t)
  57. #############################################################
  58. ## Design a low-pass filter
  59. FS = 1./dt # sampling rate
  60. FC = 10.05/(0.5*FS) # cutoff frequency at 0.05 Hz
  61. N = 11 # number of filter taps
  62. a = [1] # filter denominator
  63. b = signal.firwin(N, cutoff=FC, window='hamming') # filter numerator
  64. #############################################################
  65. ## In-phase
  66. #2*np.cos(wL*t)
  67. dw = -2.*np.pi*2
  68. Q = signal.filtfilt(b, a, xn*2*np.cos((wL+dw)*t)) # X
  69. I = signal.filtfilt(b, a, xn*2*np.sin((wL+dw)*t)) # Y
  70. ###############################################
  71. # Plots
  72. #plt.plot(ht.real)
  73. #plt.plot(ht.imag)
  74. #plt.plot(np.abs(ht))
  75. #plt.plot(gc.real)
  76. #plt.plot(gc.imag)
  77. #plt.plot(xn)
  78. #plt.plot(xn)
  79. #plt.plot(ta, xn)
  80. #plt.plot(te, abse, '-.', linewidth=2, markersize=10)
  81. #plt.plot(ta, he, '.', markersize=10 )
  82. #plt.plot(td, htd, color='green', linewidth=2)
  83. # Phase Plots
  84. #ax2 = plt.twinx()
  85. #ax2.plot(te, hp, '.', markersize=10, color='green' )
  86. #ax2.plot(te, phase, '-.', linewidth=2, markersize=10, color='green')
  87. return Q[N:-N], I[N:-N], t[N:-N]
  88. # #####################################################################
  89. # # regress raw signal
  90. #
  91. # #[peaks, times, ind] = peakPicker(xn, wL, dt)
  92. # #[a0,b0,rt20] = regressCurve(peaks, times) #,sigma2=1,intercept=True):
  93. #
  94. # dsamp = int( DT / dt) # real
  95. # # regress analytic signal
  96. # [a0,b0,rt20] = regressCurve(he[dsamp::], t[dsamp::], intercept=True) #,sigma2=1,intercept=True):
  97. # #[b0,rt20] = regressCurve(he[dsamp::], t[dsamp::], intercept=False) #,sigma2=1,intercept=True):
  98. # #[a0,b0,rt20] = regressCurve(he, t) #,sigma2=1,intercept=True):
  99. #
  100. # # regress downsampled
  101. # [a,b,rt2] = regressCurve(abse, t[dsamp::irsamp], intercept=True) #,sigma2=1,intercept=True):
  102. # #[b,rt2] = regressCurve(htd, td, intercept=False) #,sigma2=1,intercept=True):
  103. #
  104. # return irsamp, iisamp, htd, b0, rt20, ta, b, rt2, phase, td, he, dsamp
  105. # #return irsamp, iisamp, abse, a0, b0, rt20, times, a, b, rt2, phase
  106. def RotateAmplitude(X, Y, zeta, df, t):
  107. V = X + 1j*Y
  108. return np.abs(V) * np.exp( 1j * ( np.angle(V) - zeta - 2.*np.pi*df*t ) )
  109. #return np.abs(V) * np.exp( 1j * ( np.angle(V) - zeta - df*t ) )
  110. #def Gate(x, t):
  111. # pass
  112. def gateIntegrate(T2D, T2T, gpd, sigma, stackEfficiency=2.):
  113. """ Gate integrate the signal to gpd, gates per decade
  114. """
  115. # use artificial time gates so that early times are fully captured
  116. T2T0 = T2T[0]
  117. T2TD = T2T[0] - (T2T[1]-T2T[0])
  118. T2T -= T2TD
  119. # calculate total number of decades
  120. nd = np.log10(T2T[-1]/T2T[0]) #np.log10(self.T2T[-1]) - np.log10(self.T2T[-1])
  121. tdd = np.logspace( np.log10(T2T[0]), np.log10(T2T[-1]), (int)(gpd*nd)+1, base=10, endpoint=True)
  122. tdl = tdd[0:-1] # these are left edges
  123. tdr = tdd[1::] # these are left edges
  124. td = (tdl+tdr) / 2. # window centres
  125. Vars = np.ones( len(td) ) * sigma**2 #* .15**2
  126. htd = np.zeros( len(td), dtype=complex )
  127. isum = np.zeros( len(td) )
  128. ii = 0
  129. SIGSTACK = {}
  130. SIGSTACK[ii]= []
  131. for itd in range(len(T2T)):
  132. if ( T2T[itd] > tdr[ii] ):
  133. if (ii < len(td)-1):
  134. ii += 1
  135. SIGSTACK[ii] = []
  136. else:
  137. #print "overshoot??", ii
  138. break
  139. isum[ii] += 1
  140. htd[ii] += T2D[ itd ]
  141. SIGSTACK[ii].append(T2D[itd])
  142. Vars[ii] += sigma**2
  143. sigma = np.sqrt( Vars * (1/isum)**stackEfficiency )
  144. for ii in range(len(td)):
  145. if len(SIGSTACK[ii]) > 30:
  146. sigma[ii] = np.var(SIGSTACK[ii]) / ((len(SIGSTACK[ii])-1)**(1/stackEfficiency))
  147. #sigma[ii] = np.std(SIGSTACK[ii]) * ( 1./(len(SIGSTACK[ii]))**stackEfficiency)
  148. # RESET times where isum == 1
  149. ii = 0
  150. while (isum[ii] == 1):
  151. td[ii] = T2T[ii]
  152. ii += 1
  153. htd /= isum
  154. T2T += T2TD
  155. return td+T2TD, htd, tdd+T2TD, sigma, isum # centre abscissa, data, window edges, error
  156. if __name__ == "__main__":
  157. dt = 1e-4
  158. TT = 1.5
  159. t = np.arange(0, TT, dt)
  160. vL = 2057.
  161. wL = 2.*np.pi*vL
  162. wL2 = 2.*np.pi*(vL-2.5) #-2) #-2.2) # 3 Hz off
  163. zeta = -np.pi/6. #4.234
  164. t2 = .150
  165. xs = np.exp(-t/t2) * np.cos(wL2*t + zeta)
  166. xe = np.exp(-t/t2)
  167. xn = xs + np.random.normal(0,.1,len(xs))# + (np.sign(xs)
  168. # np.random.random_integers(-1,1,len(xs))*0.6*np.random.lognormal(0, .35, len(xs)) + \
  169. # np.random.random_integers(-1,1,len(xs))*.004*np.random.weibull(.25, len(xs)), 60)))
  170. # quadrature detection downsampling
  171. T = 50 # sampling period, grab every T'th oscilation
  172. DT = .002 #85 # dead time ms
  173. #[irsamp, iisamp, abse, b0, rt20, times, b, rt2, phase, tdec, he, dsamp] = quadDetect(T, vL, wL, dt, xn, DT)
  174. [Q, I, tt] = quadrature(T, vL, wL, dt, xn, DT, t)
  175. [E0,df,phi,T2] = quadratureDetect(Q, I, tt)
  176. print("df", df)
  177. D = RotateAmplitude(I, Q, phi, df, tt)
  178. fig = plt.figure(figsize=[pc2in(20), pc2in(14)]) #
  179. ax1 = fig.add_axes([.125,.2,.8,.7])
  180. #ax1.plot(tt*1e3, np.exp(-tt/t2), linewidth=2, color='black', label="actual")
  181. ax1.plot(tt*1e3, D.imag, label="CA", color='red')
  182. ax1.plot(t*1e3, xn, color='blue', alpha=.25)
  183. ax1.plot(tt*1e3, I, label="inphase", color='blue')
  184. ax1.plot(tt*1e3, Q, label="quadrature", color='green')
  185. #ax1.plot(tt*1e3, np.angle( Q + 1j*I), label="angle", color='purple')
  186. GT, GD = gateIntegrate( D.imag, tt, 10 )
  187. GT, GDR = gateIntegrate( D.real, tt, 10 )
  188. GT, GQ = gateIntegrate( Q, tt, 10 )
  189. GT, GI = gateIntegrate( I, tt, 10 )
  190. #ax1.plot(tt*1e3, np.arctan( Q/I), label="angle", color='purple')
  191. #ax1.plot(GT*1e3, np.real(GD), 'o', label="GATE", color='purple')
  192. #ax1.plot(GT*1e3, np.real(GDR), 'o', label="GATE Real", color='red')
  193. #ax1.plot(GT*1e3, np.arctan( np.real(GQ)/np.real(GI)), 'o',label="GATE ANGLE", color='magenta')
  194. ax1.set_xlabel(r"time [ms]")
  195. ax1.set_ylim( [-1.25,1.65] )
  196. #light_grey = np.array([float(248)/float(255)]*3)
  197. legend = plt.legend( frameon=True, scatterpoints=1, numpoints=1, labelspacing=0.2 )
  198. #rect = legend.get_frame()
  199. fixLeg(legend)
  200. #rect.set_color('None')
  201. #rect.set_facecolor(light_grey)
  202. #rect.set_linewidth(0.0)
  203. #rect.set_alpha(0.5)
  204. # Remove top and right axes lines ("spines")
  205. spines_to_remove = ['top', 'right']
  206. for spine in spines_to_remove:
  207. ax1.spines[spine].set_visible(False)
  208. #ax1.xaxis.set_ticks_position('none')
  209. #ax1.yaxis.set_ticks_position('none')
  210. ax1.get_xaxis().tick_bottom()
  211. ax1.get_yaxis().tick_left()
  212. plt.savefig('rotatetime.pdf',dpi=600)
  213. plt.savefig('rotatetime.eps',dpi=600)
  214. # phase part
  215. plt.figure()
  216. plt.plot( tt*1e3, D.real, label="CA", color='red' )
  217. plt.show()
  218. exit()