Surface NMR processing and inversion GUI
Vous ne pouvez pas sélectionner plus de 25 sujets Les noms de sujets doivent commencer par une lettre ou un nombre, peuvent contenir des tirets ('-') et peuvent comporter jusqu'à 35 caractères.

harmonic.py 6.8KB

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