2 frequency harmonic modelling working much better

akvo/tressel/harmonic.py

@@ -4,7 +4,6 @@ from scipy.optimize import minimize
 from scipy.linalg import lstsq as sclstsq
 import scipy.linalg as lin
 
-#def harmonicEuler ( f0, sN, fs, nK, t ): 
 def harmonicEuler ( sN, fs, t, f0, k1, kN, ks ): 
     """
     Performs inverse calculation of harmonics contaminating a signal. 
@@ -16,10 +15,9 @@ def harmonicEuler ( sN, fs, t, f0, k1, kN, ks ):
         nK = number of harmonics to calculate 
 
     """
-    #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  )
     KK = np.arange(k1, kN+1, 1/ks )
     nK = len(KK)
-    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  )
+    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)
 
     v = np.linalg.lstsq(A, sN, rcond=None) 
     alpha = np.real(v[0]) 
@@ -29,17 +27,12 @@ def harmonicEuler ( sN, fs, t, f0, k1, kN, ks ):
     phase = np.angle(v[0]) 
 
     h = np.zeros(len(t))
-    #for ik in range(nK):
-    #    h +=  2*amp[ik] * np.cos( 2.*np.pi*(ik+1) * (f0/fs) * np.arange(1, len(t)+1, 1 )  + phase[ik] )
     for ik, k in enumerate(KK):
         h +=  2*amp[ik] * np.cos( 2.*np.pi*(k) * (f0/fs) * np.arange(1, len(t)+1, 1 )  + phase[ik] )
     
     return sN-h
-    
-    res = sN-h # residual 
 
 def harmonicNorm (f0, sN, fs, t, k1, kN, ks): 
-    #print ("norm diff")
     #return np.linalg.norm( harmonicEuler(sN, fs, t, f0, k1, kN, ks)) 
     ii =  sN < (3.* np.std(sN))
     return np.linalg.norm( harmonicEuler(sN, fs, t, f0, k1, kN, ks)[ii] ) 
@@ -50,7 +43,6 @@ def minHarmonic(sN, fs, t, f0, k1, kN, ks):
     print(res)
     return harmonicEuler(sN, fs, t, res.x[0], k1, kN, ks)#[0]
 
-#def harmonicEuler2 ( f0, f1, sN, fs, nK, t ): 
 def harmonicEuler2 ( sN, fs, t, f0, f0k1, f0kN, f0ks, f1, f1k1, f1kN, f1ks ): 
     """
     Performs inverse calculation of harmonics contaminating a signal. 
@@ -62,62 +54,41 @@ def harmonicEuler2 ( sN, fs, t, f0, f0k1, f0kN, f0ks, f1, f1k1, f1kN, f1ks ):
         f0k1 = First harmonic to calulate for f0 
         f0kN = Last base harmonic to calulate for f0
         f0ks = subharmonics to calculate 
+        f1 = second base frequency of the sinusoidal noise 
+        f1k1 = First harmonic to calulate for f1
+        f1kN = Last base harmonic to calulate for f1
+        f1ks = subharmonics to calculate at f1 base frequency
     """
-    
-    #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  )
-    #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  )
-    #A = np.concatenate( (A1, A2), axis=1 )
     KK0 = np.arange(f0k1, f0kN+1, 1/f0ks)
     nK0 = len(KK0)
-    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)
+    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)
 
     KK1 = np.arange(f1k1, f1kN+1, 1/f1ks)
     nK1 = len(KK1)
     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)
-    
 
     A = np.concatenate((A0, A1), axis=1)
-    #A = A0
 
     v = np.linalg.lstsq(A, sN, rcond=None) # rcond=None) #, rcond=1e-8)
+    
     amp0 = np.abs(v[0][0:nK0])     
     phase0 = np.angle(v[0][0:nK0]) 
     amp1 = np.abs(v[0][nK0::])     
     phase1 = np.angle(v[0][nK0::]) 
-    
 
     h = np.zeros(len(t))
-    for ik in range(nK0):
-        h +=  2*amp0[ik] * np.cos( 2.*np.pi*(ik+1) * (f0/fs) * np.arange(1, len(t)+1, 1 )  + phase0[ik] ) 
-    for ik in range(nK1):
-        h +=  2*amp1[ik] * np.cos( 2.*np.pi*(ik+1) * (f1/fs) * np.arange(1, len(t)+1, 1 )  + phase1[ik] ) # + \
+    for ik, k in enumerate(KK0):
+        h +=  2*amp0[ik] * np.cos( 2.*np.pi*(k) * (f0/fs) * np.arange(1, len(t)+1, 1 ) + phase0[ik] )
+    for ik, k in enumerate(KK1):
+        h +=  2*amp1[ik] * np.cos( 2.*np.pi*(k) * (f0/fs) * np.arange(1, len(t)+1, 1 ) + phase1[ik] )
 
     return sN-h
 
 def harmonic2Norm (f0, sN, fs, t, f0k1, f0kN, f0ks, f1k1, f1kN, f1ks): 
-    #def harmonic2Norm ( f0, sN, fs, nK, t ): 
     #return np.linalg.norm(harmonicEuler2(f0[0], f0[1], sN, fs, nK, t))
     ii =  sN < (3.* np.std(sN))
     return np.linalg.norm( harmonicEuler2(sN, fs, t, f0[0], f0k1, f0kN, f0ks, f0[1], f1k1, f1kN, f1ks)[ii] ) 
 
-#def minHarmonic(f0, sN, fs, nK, t):
-#    f02 = guessf0(sN, fs)
-#    print("minHarmonic", f0, fs, nK, " guess=", f02)
-#    # CG, BFGS, Newton-CG, L-BFGS-B, TNC, SLSQP, dogleg, trust-ncg, trust-krylov, trust-exact and trust-constr
-#    res = minimize(harmonicNorm, np.array((f0)), args=(sN, fs, nK, t), jac='2-point', method='BFGS') #, jac=jacEuler) #, hess=None, bounds=None )
-#    print(res)
-#    return harmonicEuler(res.x[0], sN, fs, nK, t)#[0]
-
-
-
-#def minHarmonic2OLD(f1, f2, sN, fs, nK, t):
-    #f02 = guessf0(sN, fs)
-    #print("minHarmonic2", f0, fs, nK, " guess=", f02)
-    #methods with bounds, L-BFGS-B, TNC, SLSQP
-#    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' )
-#    print(res)
-#    return harmonicEuler2(res.x[0], res.x[1], sN, fs, nK, t) 
-
 def minHarmonic2(sN, fs, t, f0, f0k1, f0kN, f0ks, f1, f1k1, f1kN, f1ks):
     # CG, BFGS, Newton-CG, L-BFGS-B, TNC, SLSQP, dogleg, trust-ncg, trust-krylov, trust-exact and trust-constr
     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 )