Surface NMR processing and inversion GUI
Du kan inte välja fler än 25 ämnen Ämnen måste starta med en bokstav eller siffra, kan innehålla bindestreck ('-') och vara max 35 tecken långa.

adapt.py 5.8KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147
  1. import numpy as np
  2. from numpy.linalg import lstsq
  3. from numpy.linalg import norm
  4. from numpy import fft
  5. import pylab
  6. from scipy.signal import correlate
  7. def autocorr(x):
  8. #result = np.correlate(x, x, mode='full')
  9. result = correlate(x, x, mode='full')
  10. return result[result.size/2:]
  11. class AdaptiveFilter:
  12. def __init__(self, mu):
  13. self.mu = mu
  14. def adapt_filt_Ref(self, x, R, M, mu, PCA, lambda2=0.95, H0=0):
  15. """ Taken from .m file
  16. This function is written to allow the user to filter a input signal
  17. with an adaptive filter that utilizes 2 reference signals instead of
  18. the standard method which allows for only 1 reference signal.
  19. Author: Rob Clemens Date: 3/16/06
  20. Modified and ported to Python, now takes arbitray number of reference points
  21. Original public domain source
  22. https://www.mathworks.com/matlabcentral/fileexchange/10447-noise-canceling-adaptive-filter
  23. x = data array
  24. R = reference array
  25. M = number of taps
  26. mu = forgetting factor
  27. PCA = Perform PCA
  28. """
  29. #from akvo.tressel import pca
  30. import akvo.tressel.pca as pca
  31. if np.shape(x) != np.shape(R[0]): # or np.shape(x) != np.shape(rx1):
  32. print ("Error, non aligned")
  33. exit(1)
  34. if PCA == "Yes":
  35. #print("Performing PCA calculation in noise cancellation")
  36. # PCA decomposition on ref channels so signals are less related
  37. R, K, means = pca.pca( R )
  38. # test for in loop reference
  39. #print("Cull nearly zero terms?", np.shape(x), np.shape(R))
  40. #R = R[0:3,:]
  41. #R = R[2:4,:]
  42. #print(" removed zero terms?", np.shape(x), np.shape(R))
  43. #H0 = H0[0:3*np.shape(x)[0]]
  44. #H0 = H0[0:2*np.shape(x)[0]]
  45. if all(H0) == 0:
  46. # corrects for dimensionality issues if a simple 0 is passed
  47. H = np.zeros( (len(R)*M))
  48. else:
  49. H = H0
  50. Rn = np.ones(len(R)*M) / mu
  51. r_ = np.zeros( (len(R), M) )
  52. e = np.zeros(len(x)) # error, in our case the desired output
  53. ilambda = lambda2**-1
  54. for ix in range(0, len(x)):
  55. # Only look forwards, to avoid distorting the lates times
  56. # (run backwards, if opposite and you don't care about distorting very late time.)
  57. for ir in range(len(R)): # number of reference channels
  58. if ix < M:
  59. r_[ir,0:ix] = R[ir][0:ix]
  60. r_[ir,ix:M] = 0
  61. else:
  62. r_[ir,:] = R[ir][ix-M:ix]
  63. # reshape
  64. r_n = np.reshape(r_, -1) # concatenate the ref channels in to a 1D array
  65. K = (Rn* r_n) / (lambda2 + np.dot(r_n*Rn, r_n)) # Create/update K
  66. e[ix] = x[ix] - np.dot(r_n.T, H) # e is the filtered signal, input - r(n) * Filter Coefs
  67. H += K*e[ix]; # Update Filter Coefficients
  68. Rn = ilambda*Rn - ilambda*np.dot(np.dot(K, r_n.T), Rn) # Update R(n)
  69. return e, H
  70. def transferFunctionFFT(self, D, R, reg=1e-2):
  71. from akvo.tressel import pca
  72. """
  73. Computes the transfer function (H) between a Data channel and
  74. a number of Reference channels. The Matrices D and R are
  75. expected to be in the frequency domain on input.
  76. | R1'R1 R1'R2 R1'R3| |h1| |R1'D|
  77. | R2'R1 R2'R2 R2'R3| * |h2| = |R2'D|
  78. | R3'R1 R3'R2 R3'R3| |h3| |R3'D|
  79. Returns the corrected array
  80. """
  81. # PCA decomposition on ref channels so signals are less related
  82. #transMatrix, K, means = pca.pca( np.array([rx0, rx1]))
  83. #RR = np.zeros(( np.shape(R[0])[0]*np.shape(R[0])[1], len(R)))
  84. # RR = np.zeros(( len(R), np.shape(R[0])[0]*np.shape(R[0])[1] ))
  85. # for ir in range(len(R)):
  86. # RR[ir,:] = np.reshape(R[ir], -1)
  87. # transMatrix, K, means = pca.pca(RR)
  88. # #R rx0 = transMatrix[0,:]
  89. # # rx1 = transMatrix[1,:]
  90. # for ir in range(len(R)):
  91. # R[ir] = transMatrix[ir,0]
  92. import scipy.linalg
  93. import akvo.tressel.pca as pca
  94. # Compute as many transfer functions as len(R)
  95. # A*H = B
  96. nref = len(R)
  97. H = np.zeros( (np.shape(D)[1], len(R)), dtype=complex )
  98. for iw in range(np.shape(D)[1]):
  99. A = np.zeros( (nref, nref), dtype=complex )
  100. B = np.zeros( (nref) , dtype=complex)
  101. for ii in range(nref):
  102. for jj in range(nref):
  103. # build A
  104. A[ii,jj] = np.dot(R[ii][:,iw], R[jj][:,iw])
  105. # build B
  106. B[ii] = np.dot( R[ii][:,iw], D[:,iw] )
  107. # compute H(iw)
  108. #linalg.solve(a,b) if a is square
  109. #print "A", A
  110. #print "B", B
  111. # TODO, regularise this solve step? So as to not fit the spurious noise
  112. #print np.shape(B), np.shape(A)
  113. #H[iw, :] = scipy.linalg.solve(A,B)
  114. H[iw, :] = scipy.linalg.lstsq(A,B,cond=reg)[0]
  115. #print "lstqt", np.shape(scipy.linalg.lstsq(A,B))
  116. #print "solve", scipy.linalg.solve(A,B)
  117. #H[iw,:] = scipy.linalg.lstsq(A,B) # otherwise
  118. #H = np.zeros( (np.shape(D)[1], ) )
  119. #print H #A, B
  120. Error = np.zeros(np.shape(D), dtype=complex)
  121. for ir in range(nref):
  122. for q in range( np.shape(D)[0] ):
  123. #print "dimcheck", np.shape(H[:,ir]), np.shape(R[ir][q,:] )
  124. Error[q,:] += H[:,ir]*R[ir][q,:]
  125. return D - Error