Surface NMR processing and inversion GUI
You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.

rotate.py 9.9KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264
  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. T2D = the time series to gate integrate, complex
  115. T2T = the abscissa values
  116. gpd = gates per decade
  117. sigma = estimate of standard deviation for theoretical gate noise
  118. stackEfficiency = exponential in theoretical gate noise, 2 represents ideal stacking
  119. """
  120. # use artificial time gates so that early times are fully captured
  121. T2T0 = T2T[0]
  122. T2TD = T2T[0] - (T2T[1]-T2T[0])
  123. T2T -= T2TD
  124. #####################################
  125. # calculate total number of decades #
  126. # windows edges are approximate until binning but will be adjusted to reflect data timing, this
  127. # primarily impacts bins with a few samples
  128. nd = np.log10(T2T[-1]/T2T[0])
  129. tdd = np.logspace( np.log10(T2T[0]), np.log10(T2T[-1]), (int)(gpd*nd)+1, base=10, endpoint=True)
  130. tdl = tdd[0:-1] # approximate window left edges
  131. tdr = tdd[1::] # approximate window right edges
  132. td = (tdl+tdr) / 2. # approximate window centres
  133. Vars = np.zeros( len(td) )
  134. htd = np.zeros( len(td), dtype=complex )
  135. isum = np.zeros( len(td), dtype=int )
  136. ii = 0
  137. for itd in range(len(T2T)):
  138. if ( T2T[itd] > tdr[ii] ):
  139. ii += 1
  140. # correct window edges to centre about data
  141. tdr[ii-1] = (T2T[itd-1]+T2T[itd])*.5
  142. tdl[ii ] = (T2T[itd-1]+T2T[itd])*.5
  143. isum[ii] += 1
  144. htd[ii] += T2D[ itd ]
  145. Vars[ii] += sigma**2
  146. td = (tdl+tdr) / 2. # actual window centres
  147. sigma2 = np.sqrt( Vars * ((1/(isum))**stackEfficiency) )
  148. # Reset abscissa where isum == 1
  149. # when there is no windowing going on
  150. td[isum==1] = T2T[0:len(td)][isum==1]
  151. tdd = np.append(tdl, tdr[-1])
  152. htd /= isum # average
  153. T2T += T2TD
  154. return td+T2TD, htd, tdd+T2TD, sigma2, isum # centre abscissa, data, window edges, error
  155. if __name__ == "__main__":
  156. dt = 1e-4
  157. TT = 1.5
  158. t = np.arange(0, TT, dt)
  159. vL = 2057.
  160. wL = 2.*np.pi*vL
  161. wL2 = 2.*np.pi*(vL-2.5) #-2) #-2.2) # 3 Hz off
  162. zeta = -np.pi/6. #4.234
  163. t2 = .150
  164. xs = np.exp(-t/t2) * np.cos(wL2*t + zeta)
  165. xe = np.exp(-t/t2)
  166. xn = xs + np.random.normal(0,.1,len(xs))# + (np.sign(xs)
  167. # np.random.random_integers(-1,1,len(xs))*0.6*np.random.lognormal(0, .35, len(xs)) + \
  168. # np.random.random_integers(-1,1,len(xs))*.004*np.random.weibull(.25, len(xs)), 60)))
  169. # quadrature detection downsampling
  170. T = 50 # sampling period, grab every T'th oscilation
  171. DT = .002 #85 # dead time ms
  172. #[irsamp, iisamp, abse, b0, rt20, times, b, rt2, phase, tdec, he, dsamp] = quadDetect(T, vL, wL, dt, xn, DT)
  173. [Q, I, tt] = quadrature(T, vL, wL, dt, xn, DT, t)
  174. [E0,df,phi,T2] = quadratureDetect(Q, I, tt)
  175. print("df", df)
  176. D = RotateAmplitude(I, Q, phi, df, tt)
  177. fig = plt.figure(figsize=[pc2in(20), pc2in(14)]) #
  178. ax1 = fig.add_axes([.125,.2,.8,.7])
  179. #ax1.plot(tt*1e3, np.exp(-tt/t2), linewidth=2, color='black', label="actual")
  180. ax1.plot(tt*1e3, D.imag, label="CA", color='red')
  181. ax1.plot(t*1e3, xn, color='blue', alpha=.25)
  182. ax1.plot(tt*1e3, I, label="inphase", color='blue')
  183. ax1.plot(tt*1e3, Q, label="quadrature", color='green')
  184. #ax1.plot(tt*1e3, np.angle( Q + 1j*I), label="angle", color='purple')
  185. GT, GD = gateIntegrate( D.imag, tt, 10 )
  186. GT, GDR = gateIntegrate( D.real, tt, 10 )
  187. GT, GQ = gateIntegrate( Q, tt, 10 )
  188. GT, GI = gateIntegrate( I, tt, 10 )
  189. #ax1.plot(tt*1e3, np.arctan( Q/I), label="angle", color='purple')
  190. #ax1.plot(GT*1e3, np.real(GD), 'o', label="GATE", color='purple')
  191. #ax1.plot(GT*1e3, np.real(GDR), 'o', label="GATE Real", color='red')
  192. #ax1.plot(GT*1e3, np.arctan( np.real(GQ)/np.real(GI)), 'o',label="GATE ANGLE", color='magenta')
  193. ax1.set_xlabel(r"time [ms]")
  194. ax1.set_ylim( [-1.25,1.65] )
  195. #light_grey = np.array([float(248)/float(255)]*3)
  196. legend = plt.legend( frameon=True, scatterpoints=1, numpoints=1, labelspacing=0.2 )
  197. #rect = legend.get_frame()
  198. fixLeg(legend)
  199. #rect.set_color('None')
  200. #rect.set_facecolor(light_grey)
  201. #rect.set_linewidth(0.0)
  202. #rect.set_alpha(0.5)
  203. # Remove top and right axes lines ("spines")
  204. spines_to_remove = ['top', 'right']
  205. for spine in spines_to_remove:
  206. ax1.spines[spine].set_visible(False)
  207. #ax1.xaxis.set_ticks_position('none')
  208. #ax1.yaxis.set_ticks_position('none')
  209. ax1.get_xaxis().tick_bottom()
  210. ax1.get_yaxis().tick_left()
  211. plt.savefig('rotatetime.pdf',dpi=600)
  212. plt.savefig('rotatetime.eps',dpi=600)
  213. # phase part
  214. plt.figure()
  215. plt.plot( tt*1e3, D.real, label="CA", color='red' )
  216. plt.show()
  217. exit()