Surface NMR processing and inversion GUI
Nelze vybrat více než 25 témat Téma musí začínat písmenem nebo číslem, může obsahovat pomlčky („-“) a může být dlouhé až 35 znaků.

harmonic.py 5.8KB

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. def harmonic2 ( f1, f2, sN, fs, nK, t ):
  6. """
  7. Performs inverse calculation of two harmonics contaminating a signal.
  8. Args:
  9. f01 = base frequency of the first sinusoidal noise
  10. f02 = base frequency of the second sinusoidal noise
  11. sN = signal containing noise
  12. fs = sampling frequency
  13. nK = number of harmonics to calculate
  14. t = time samples
  15. """
  16. print("building matrix ")
  17. A = np.zeros( (len(t), 4*nK) )
  18. #f1 = f1MHz * 1e-3
  19. #f2 = f2MHz * 1e-3
  20. for irow, tt in enumerate(t):
  21. A[irow, 0:2*nK:2] = np.cos( np.arange(nK)*2*np.pi*(f1/fs)*irow )
  22. A[irow, 1:2*nK:2] = np.sin( np.arange(nK)*2*np.pi*(f1/fs)*irow )
  23. A[irow, 2*nK::2] = np.cos( np.arange(nK)*2*np.pi*(f2/fs)*irow )
  24. A[irow, 2*nK+1::2] = np.sin( np.arange(nK)*2*np.pi*(f2/fs)*irow )
  25. v = np.linalg.lstsq(A, sN, rcond=1e-8)
  26. #v = sclstsq(A, sN) #, rcond=1e-6)
  27. alpha = v[0][0:2*nK:2]
  28. beta = v[0][1:2*nK:2]
  29. amp = np.sqrt( alpha**2 + beta**2 )
  30. phase = np.arctan(- beta/alpha)
  31. alpha2 = v[0][2*nK::2]
  32. beta2 = v[0][2*nK+1::2]
  33. amp2 = np.sqrt( alpha2**2 + beta2**2 )
  34. phase2 = np.arctan(- beta2/alpha2)
  35. h = np.zeros(len(t))
  36. for ik in range(nK):
  37. h += np.sqrt(alpha[ik]**2 + beta[ik]**2) * np.cos( 2.*np.pi*ik * (f1/fs) * np.arange(0, len(t), 1 ) + phase[ik] ) \
  38. + np.sqrt(alpha2[ik]**2 + beta2[ik]**2) * np.cos( 2.*np.pi*ik * (f2/fs) * np.arange(0, len(t), 1 ) + phase2[ik] )
  39. return sN-h
  40. def harmonic ( f0, sN, fs, nK, t ):
  41. """
  42. Performs inverse calculation of harmonics contaminating a signal.
  43. Args:
  44. f0 = base frequency of the sinusoidal noise
  45. sN = signal containing noise
  46. fs = sampling frequency
  47. nK = number of harmonics to calculate
  48. t = time samples
  49. """
  50. print("building matrix ")
  51. A = np.zeros( (len(t), 2*nK) )
  52. for irow, tt in enumerate(t):
  53. A[irow, 0::2] = np.cos( np.arange(nK)*2*np.pi*(f0/fs)*irow )
  54. A[irow, 1::2] = np.sin( np.arange(nK)*2*np.pi*(f0/fs)*irow )
  55. v = np.linalg.lstsq(A, sN, rcond=None) #, rcond=1e-8)
  56. alpha = v[0][0::2]
  57. beta = v[0][1::2]
  58. amp = np.sqrt( alpha**2 + beta**2 )
  59. phase = np.arctan(- beta/alpha)
  60. #print("amp:", amp, " phase", phase)
  61. h = np.zeros(len(t))
  62. for ik in range(nK):
  63. h += np.sqrt(alpha[ik]**2 + beta[ik]**2) * np.cos( 2.*np.pi*ik * (f0/fs) * np.arange(0, len(t), 1 ) + phase[ik] )
  64. #plt.matshow(A, aspect='auto')
  65. #plt.colorbar()
  66. #plt.figure()
  67. #plt.plot(alpha)
  68. #plt.plot(beta)
  69. #plt.plot(amp)
  70. #plt.figure()
  71. #plt.plot(h)
  72. #plt.title("modelled noise")
  73. return sN-h
  74. def jacobian( f0, sN, fs, nK, t):
  75. print("building Jacobian matrix ")
  76. A = np.zeros( (len(t), 2*nK) )
  77. for irow, tt in enumerate(t):
  78. #A[irow, 0::2] = np.cos( np.arange(nK)*2*np.pi*(f0/fs)*irow )
  79. #A[irow, 1::2] = np.sin( np.arange(nK)*2*np.pi*(f0/fs)*irow )
  80. # brutal
  81. for k, ik in enumerate( np.arange(0, 2*nK, 2) ):
  82. #A[irow, ik ] = np.cos( k*2*np.pi*(f0/fs)*irow )
  83. #A[irow, ik+1] = np.sin( k*2*np.pi*(f0/fs)*irow )
  84. A[irow, ik ] = - (2.*np.pi*k*irow * sin((2.*np.pi*irow*f0)/fs)) / fs
  85. A[irow, ik+1] = (2.*np.pi*k*irow * cos((2.*np.pi*irow*f0)/fs)) / fs
  86. def harmonicNorm ( f0, sN, fs, nK, t ):
  87. return np.linalg.norm( harmonic(f0, sN, fs, nK, t))
  88. def harmonic2Norm ( f0, sN, fs, nK, t ):
  89. return np.linalg.norm(harmonic2(f0[0], f0[1], sN, fs, nK, t))
  90. def minHarmonic(f0, sN, fs, nK, t):
  91. f02 = guessf0(sN, fs)
  92. print("minHarmonic", f0, fs, nK, " guess=", f02)
  93. res = minimize( harmonicNorm, np.array((f0)), args=(sN, fs, nK, t)) #, method='Nelder-Mead' )# jac=None, hess=None, bounds=None )
  94. print(res)
  95. return harmonic(res.x[0], sN, fs, nK, t)
  96. def minHarmonic2(f1, f2, sN, fs, nK, t):
  97. #f02 = guessf0(sN, fs)
  98. #print("minHarmonic2", f0, fs, nK, " guess=", f02)
  99. #methods with bounds, L-BFGS-B, TNC, SLSQP
  100. res = minimize( harmonic2Norm, np.array((f1,f2)), args=(sN, fs, nK, t)) #, bounds=((f1-1.,f1+1.0),(f2-1.0,f2+1.0)), method='SLSQP' )
  101. print(res)
  102. return harmonic2(res.x[0], res.x[1], sN, fs, nK, t)
  103. def guessf0( sN, fs ):
  104. S = np.fft.fft(sN)
  105. w = np.fft.fftfreq( len(sN), 1/fs )
  106. imax = np.argmax( np.abs(S) )
  107. #plt.plot( w, np.abs(S) )
  108. #plt.show()
  109. #print(w)
  110. #print ( w[imax], w[imax+1] )
  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 = 0 #np.random.rand()
  117. delta2 = 0 #np.random.rand()
  118. print("delta", delta)
  119. fs = 10000 # GMR
  120. t = np.arange(0, 1, 1/fs)
  121. phi = 0 #np.random.rand()
  122. phi2 = 0 # np.random.rand()
  123. A = 1.0
  124. A2 = 0.0
  125. nK = 10
  126. T2 = .200
  127. sN = A * np.sin( ( 1*(delta +f0))*2*np.pi*t + phi ) + \
  128. A2* np.sin( ( 1*(delta2 +f1))*2*np.pi*t + phi2 ) + \
  129. np.random.normal(0,.1,len(t)) + \
  130. + np.exp( -t/T2 )
  131. sNc = A * np.sin( (1*(delta +f0))*2*np.pi*t + phi ) + \
  132. A2* np.sin( (1*(delta2+f1))*2*np.pi*t + phi2 ) + \
  133. + np.exp( -t/T2 )
  134. guessf0(sN, fs)
  135. #h = harmonic( f0, sN, fs, nK, t)
  136. #h = minHarmonic2( f0, f1, sN, fs, nK, t)
  137. h = harmonic2( f0, f1, sN, fs, nK, t)
  138. plt.figure()
  139. plt.plot(t, sN, label="sN")
  140. #plt.plot(t, sN-h, label="sN-h")
  141. plt.plot(t, h, label='h')
  142. plt.title("harmonic")
  143. plt.legend()
  144. plt.figure()
  145. plt.plot(t, sN-sNc, label='true noise')
  146. plt.plot(t, h, label='harmonic removal')
  147. plt.plot(t, np.exp(-t/T2), label="nmr")
  148. plt.legend()
  149. plt.title("true noise")
  150. plt.show()