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 5.6KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163
  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):
  33. # CG, BFGS, Newton-CG, L-BFGS-B, TNC, SLSQP, dogleg, trust-ncg, trust-krylov, trust-exact and trust-constr
  34. res = minimize(harmonicNorm, np.array((f0)), args=(sN, fs, t, k1, kN, ks), jac='2-point', method='BFGS') # hess=None, bounds=None )
  35. #print(res)
  36. return harmonicEuler(sN, fs, t, res.x[0], k1, kN, ks)#[0]
  37. def harmonicEuler2 ( sN, fs, t, f0, f0k1, f0kN, f0ks, f1, f1k1, f1kN, f1ks ):
  38. """
  39. Performs inverse calculation of harmonics contaminating a signal.
  40. Args:
  41. sN = signal containing noise
  42. fs = sampling frequency
  43. t = time samples
  44. f0 = first base frequency of the sinusoidal noise
  45. f0k1 = First harmonic to calulate for f0
  46. f0kN = Last base harmonic to calulate for f0
  47. f0ks = subharmonics to calculate
  48. f1 = second base frequency of the sinusoidal noise
  49. f1k1 = First harmonic to calulate for f1
  50. f1kN = Last base harmonic to calulate for f1
  51. f1ks = subharmonics to calculate at f1 base frequency
  52. """
  53. KK0 = np.arange(f0k1, f0kN+1, 1/f0ks)
  54. nK0 = len(KK0)
  55. 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)
  56. KK1 = np.arange(f1k1, f1kN+1, 1/f1ks)
  57. nK1 = len(KK1)
  58. 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)
  59. A = np.concatenate((A0, A1), axis=1)
  60. v = np.linalg.lstsq(A, sN, rcond=None) # rcond=None) #, rcond=1e-8)
  61. amp0 = np.abs(v[0][0:nK0])
  62. phase0 = np.angle(v[0][0:nK0])
  63. amp1 = np.abs(v[0][nK0::])
  64. phase1 = np.angle(v[0][nK0::])
  65. h = np.zeros(len(t))
  66. for ik, k in enumerate(KK0):
  67. h += 2*amp0[ik] * np.cos( 2.*np.pi*(k) * (f0/fs) * np.arange(1, len(t)+1, 1 ) + phase0[ik] )
  68. for ik, k in enumerate(KK1):
  69. h += 2*amp1[ik] * np.cos( 2.*np.pi*(k) * (f0/fs) * np.arange(1, len(t)+1, 1 ) + phase1[ik] )
  70. return sN-h
  71. def harmonic2Norm (f0, sN, fs, t, f0k1, f0kN, f0ks, f1k1, f1kN, f1ks):
  72. #return np.linalg.norm(harmonicEuler2(f0[0], f0[1], sN, fs, nK, t))
  73. ii = sN < (3.* np.std(sN))
  74. return np.linalg.norm( harmonicEuler2(sN, fs, t, f0[0], f0k1, f0kN, f0ks, f0[1], f1k1, f1kN, f1ks)[ii] )
  75. def minHarmonic2(sN, fs, t, f0, f0k1, f0kN, f0ks, f1, f1k1, f1kN, f1ks):
  76. # CG, BFGS, Newton-CG, L-BFGS-B, TNC, SLSQP, dogleg, trust-ncg, trust-krylov, trust-exact and trust-constr
  77. 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 )
  78. print(res)
  79. return harmonicEuler2(sN, fs, t, res.x[0], f0k1, f0kN, f0ks, res.x[1], f1k1, f1kN, f1ks)#[0]
  80. def guessf0( sN, fs ):
  81. S = np.fft.fft(sN)
  82. w = np.fft.fftfreq( len(sN), 1/fs )
  83. imax = np.argmax( np.abs(S) )
  84. #plt.plot( w, np.abs(S) )
  85. #plt.show()
  86. #print(w)
  87. #print ( w[imax], w[imax+1] )
  88. return abs(w[imax])
  89. if __name__ == "__main__":
  90. import matplotlib.pyplot as plt
  91. f0 = 60 # Hz
  92. f1 = 60 # Hz
  93. delta = np.random.rand() - .5
  94. delta2 = np.random.rand() - .5
  95. print("delta", delta)
  96. print("delta2", delta2)
  97. fs = 10000 # GMR
  98. t = np.arange(0, 1, 1/fs)
  99. phi = 2.*np.pi*np.random.rand() - np.pi
  100. phi2 = 2.*np.pi*np.random.rand() - np.pi
  101. print("phi", phi, phi2)
  102. A = 1.0
  103. A2 = 0.0
  104. A3 = 1.0
  105. nK = 10
  106. T2 = .200
  107. sN = A *np.sin( ( 1*(delta +f0))*2*np.pi*t + phi ) + \
  108. A2*np.sin( ( 1*(delta2 +f1))*2*np.pi*t + phi2 ) + \
  109. np.random.normal(0,.1,len(t)) + \
  110. + A3*np.exp( -t/T2 )
  111. sNc = A *np.sin( (1*(delta +f0))*2*np.pi*t + phi ) + \
  112. A2*np.sin( (1*(delta2+f1))*2*np.pi*t + phi2 ) + \
  113. + A3*np.exp( -t/T2 )
  114. guessf0(sN, fs)
  115. # single freq
  116. #h = harmonicEuler( f0, sN, fs, nK, t)
  117. h = minHarmonic( f0, sN, fs, nK, t)
  118. # two freqs
  119. #h = minHarmonic2( f0+1e-2, f1-1e-2, sN, fs, nK, t)
  120. #h = harmonicEuler2( f0, f1, sN, fs, nK, t)
  121. plt.figure()
  122. plt.plot(t, sN, label="sN")
  123. #plt.plot(t, sN-h, label="sN-h")
  124. plt.plot(t, h, label='h')
  125. plt.title("harmonic")
  126. plt.legend()
  127. plt.figure()
  128. plt.plot(t, sN-sNc, label='true noise')
  129. plt.plot(t, h, label='harmonic removal')
  130. plt.plot(t, np.exp(-t/T2), label="nmr")
  131. plt.legend()
  132. plt.title("true noise")
  133. plt.show()