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

rotate.py 11KB

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