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.

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183
  1. import numpy as np
  2. from scipy.optimize import least_squares
  3. from scipy.optimize import minimize
  4. from scipy.linalg import lstsq as sclstsq
  5. import scipy.linalg as lin
  6. def harmonicEuler ( f0, sN, fs, nK, t ):
  7. """
  8. Performs inverse calculation of harmonics contaminating a signal.
  9. Args:
  10. f0 = base frequency of the sinusoidal noise
  11. sN = signal containing noise
  12. fs = sampling frequency
  13. nK = number of harmonics to calculate
  14. t = time samples
  15. """
  16. A = np.exp(1j* np.tile( np.arange(1,nK+1),(len(t), 1)) * 2*np.pi* (f0/fs) * np.tile(np.arange(len(t)),(nK,1)).T )
  17. v = np.linalg.lstsq(A, sN, rcond=None)
  18. alpha = np.real(v[0])
  19. beta = np.imag(v[0])
  20. amp = np.abs(v[0])
  21. phase = np.angle(v[0])
  22. h = np.zeros(len(t))
  23. for ik in range(nK):
  24. h += 2*amp[ik] * np.cos( 2.*np.pi*(ik+1) * (f0/fs) * np.arange(0, len(t), 1 ) + phase[ik] )
  25. return sN-h
  26. res = sN-h # residual
  27. # calculate jacobian
  28. #J = jacEuler(f0, sN, fs, nK, t)
  29. #plt.matshow(np.real(J), aspect='auto')
  30. #Jv = J * np.tile(v[0], (len(t) ,1))
  31. #plt.matshow(np.real(Jv), aspect='auto')
  32. #print ("shape J", np.shape(J))
  33. #print ("shape v", np.shape(v[0]))
  34. #print ("shape Jv", np.shape(Jv))
  35. #plt.figure()
  36. #plt.plot(v[0])
  37. #plt.figure()
  38. #plt.plot(Jv)
  39. #plt.show()
  40. #exit()
  41. #jac = np.linalg.norm( np.dot(Jv.T, res))
  42. #print("norm ||", jac , "||", " at freq", f0)
  43. #return res, jac
  44. def jacEuler( f0, sN, fs, nK, t):
  45. print("building Jacobian matrix ")
  46. J = np.zeros( (len(t), nK), dtype=np.complex64 )
  47. for it, tt in enumerate(t):
  48. for ik, k in enumerate( np.arange(0, nK) ):
  49. c = 1j*2.*np.pi*(ik+1.)
  50. J[it, ik] = (c*np.exp(it*c*f0/fs)) / fs
  51. return J
  52. def harmonicEuler2 ( f0, f1, sN, fs, nK, t ):
  53. """
  54. Performs inverse calculation of harmonics contaminating a signal.
  55. Args:
  56. f0 = base frequency of the sinusoidal noise
  57. sN = signal containing noise
  58. fs = sampling frequency
  59. nK = number of harmonics to calculate
  60. t = time samples
  61. """
  62. A1 = np.exp(1j* np.tile( np.arange(1,nK+1),(len(t), 1)) * 2*np.pi* (f0/fs) * np.tile(np.arange(len(t)),(nK,1)).T )
  63. A2 = np.exp(1j* np.tile( np.arange(1,nK+1),(len(t), 1)) * 2*np.pi* (f1/fs) * np.tile(np.arange(len(t)),(nK,1)).T )
  64. A = np.concatenate( (A1, A2), axis=1 )
  65. v = np.linalg.lstsq(A, sN, rcond=None) # rcond=None) #, rcond=1e-8)
  66. amp = np.abs(v[0][0:nK])
  67. phase = np.angle(v[0][0:nK])
  68. amp1 = np.abs(v[0][nK:2*nK])
  69. phase1 = np.angle(v[0][nK:2*nK])
  70. h = np.zeros(len(t))
  71. for ik in range(nK):
  72. h += 2*amp[ik] * np.cos( 2.*np.pi*(ik+1) * (f0/fs) * np.arange(0, len(t), 1 ) + phase[ik] ) + \
  73. 2*amp1[ik] * np.cos( 2.*np.pi*(ik+1) * (f1/fs) * np.arange(0, len(t), 1 ) + phase1[ik] )
  74. return sN-h
  75. def harmonicNorm ( f0, sN, fs, nK, t ):
  76. return np.linalg.norm( harmonicEuler(f0, sN, fs, nK, t) )
  77. def harmonic2Norm ( f0, sN, fs, nK, t ):
  78. return np.linalg.norm(harmonicEuler2(f0[0], f0[1], sN, fs, nK, t))
  79. def minHarmonic(f0, sN, fs, nK, t):
  80. f02 = guessf0(sN, fs)
  81. print("minHarmonic", f0, fs, nK, " guess=", f02)
  82. # CG, BFGS, Newton-CG, L-BFGS-B, TNC, SLSQP, dogleg, trust-ncg, trust-krylov, trust-exact and trust-constr
  83. res = minimize(harmonicNorm, np.array((f02)), args=(sN, fs, nK, t), jac='2-point', method='BFGS') #, jac=jacEuler) #, hess=None, bounds=None )
  84. print(res)
  85. return harmonicEuler(res.x[0], sN, fs, nK, t)#[0]
  86. def minHarmonic2(f1, f2, sN, fs, nK, t):
  87. #f02 = guessf0(sN, fs)
  88. #print("minHarmonic2", f0, fs, nK, " guess=", f02)
  89. #methods with bounds, L-BFGS-B, TNC, SLSQP
  90. res = minimize( harmonic2Norm, np.array((f1,f2)), args=(sN, fs, nK, t), jac='2-point', method='BFGS') #, bounds=((f1-1.,f1+1.0),(f2-1.0,f2+1.0)), method='TNC' )
  91. print(res)
  92. return harmonicEuler2(res.x[0], res.x[1], sN, fs, nK, t)
  93. def guessf0( sN, fs ):
  94. S = np.fft.fft(sN)
  95. w = np.fft.fftfreq( len(sN), 1/fs )
  96. imax = np.argmax( np.abs(S) )
  97. #plt.plot( w, np.abs(S) )
  98. #plt.show()
  99. #print(w)
  100. #print ( w[imax], w[imax+1] )
  101. return abs(w[imax])
  102. if __name__ == "__main__":
  103. import matplotlib.pyplot as plt
  104. f0 = 60 # Hz
  105. f1 = 60 # Hz
  106. delta = np.random.rand() - .5
  107. delta2 = np.random.rand() - .5
  108. print("delta", delta)
  109. print("delta2", delta2)
  110. fs = 10000 # GMR
  111. t = np.arange(0, 1, 1/fs)
  112. phi = 2.*np.pi*np.random.rand() - np.pi
  113. phi2 = 2.*np.pi*np.random.rand() - np.pi
  114. print("phi", phi, phi2)
  115. A = 1.0
  116. A2 = 1.0
  117. A3 = 1.0
  118. nK = 10
  119. T2 = .200
  120. sN = A *np.sin( ( 1*(delta +f0))*2*np.pi*t + phi ) + \
  121. A2*np.sin( ( 1*(delta2 +f1))*2*np.pi*t + phi2 ) + \
  122. np.random.normal(0,.1,len(t)) + \
  123. + A3*np.exp( -t/T2 )
  124. sNc = A *np.sin( (1*(delta +f0))*2*np.pi*t + phi ) + \
  125. A2*np.sin( (1*(delta2+f1))*2*np.pi*t + phi2 ) + \
  126. + A3*np.exp( -t/T2 )
  127. guessf0(sN, fs)
  128. # single freq
  129. #h = harmonicEuler( f0, sN, fs, nK, t)
  130. #h = minHarmonic( f0, sN, fs, nK, t)
  131. # two freqs
  132. h = minHarmonic2( f0+1e-2, f1-1e-2, sN, fs, nK, t)
  133. #h = harmonicEuler2( f0, f1, sN, fs, nK, t)
  134. plt.figure()
  135. plt.plot(t, sN, label="sN")
  136. #plt.plot(t, sN-h, label="sN-h")
  137. plt.plot(t, h, label='h')
  138. plt.title("harmonic")
  139. plt.legend()
  140. plt.figure()
  141. plt.plot(t, sN-sNc, label='true noise')
  142. plt.plot(t, h, label='harmonic removal')
  143. plt.plot(t, np.exp(-t/T2), label="nmr")
  144. plt.legend()
  145. plt.title("true noise")
  146. plt.show()