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.

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