Переглянути джерело

Fixed bug in harmonic with two frequencies

tags/1.6.1
Trevor Irons 5 роки тому
джерело
коміт
eadd0e4ef8
1 змінених файлів з 16 додано та 11 видалено
  1. 16
    11
      akvo/tressel/harmonic.py

+ 16
- 11
akvo/tressel/harmonic.py Переглянути файл

@@ -63,29 +63,34 @@ def harmonicEuler2 ( sN, fs, t, f0, f0k1, f0kN, f0ks, f1, f1k1, f1kN, f1ks ):
63 63
         f0kN = Last base harmonic to calulate for f0
64 64
         f0ks = subharmonics to calculate 
65 65
     """
66
+    
66 67
     #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  )
67 68
     #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  )
68 69
     #A = np.concatenate( (A1, A2), axis=1 )
69
-    KK0 = np.arange(f0k1, f0kN+1, 1/f0ks )
70
+    KK0 = np.arange(f0k1, f0kN+1, 1/f0ks)
70 71
     nK0 = len(KK0)
71
-    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  )
72
-    KK1 = np.arange(f1k1, f1kN+1, 1/f1ks )
72
+    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)
73
+
74
+    KK1 = np.arange(f1k1, f1kN+1, 1/f1ks)
73 75
     nK1 = len(KK1)
74
-    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  )
75
-    A = np.concatenate( (A0, A1), axis=1 )
76
+    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)
77
+    
78
+
79
+    A = np.concatenate((A0, A1), axis=1)
80
+    #A = A0
76 81
 
77 82
     v = np.linalg.lstsq(A, sN, rcond=None) # rcond=None) #, rcond=1e-8)
78
-    amp = np.abs(v[0][0:nK0])     
79
-    phase = np.angle(v[0][0:nK0]) 
83
+    amp0 = np.abs(v[0][0:nK0])     
84
+    phase0 = np.angle(v[0][0:nK0]) 
80 85
     amp1 = np.abs(v[0][nK0::])     
81 86
     phase1 = np.angle(v[0][nK0::]) 
87
+    
82 88
 
83 89
     h = np.zeros(len(t))
84 90
     for ik in range(nK0):
85
-        h +=  2*amp[ik]  * np.cos( 2.*np.pi*(ik+1) * (f0/fs) * np.arange(1, len(t)+1, 1 )  + phase[ik] ) 
91
+        h +=  2*amp0[ik] * np.cos( 2.*np.pi*(ik+1) * (f0/fs) * np.arange(1, len(t)+1, 1 )  + phase0[ik] ) 
86 92
     for ik in range(nK1):
87
-        h +=  2*amp1[ik]  * np.cos( 2.*np.pi*(ik+1) * (f1/fs) * np.arange(1, len(t)+1, 1 )  + phase1[ik] ) # + \
88
-        #      2*amp1[ik] * np.cos( 2.*np.pi*(ik+1) * (f1/fs) * np.arange(1, len(t)+1, 1 )  + phase1[ik] )
93
+        h +=  2*amp1[ik] * np.cos( 2.*np.pi*(ik+1) * (f1/fs) * np.arange(1, len(t)+1, 1 )  + phase1[ik] ) # + \
89 94
 
90 95
     return sN-h
91 96
 
@@ -117,7 +122,7 @@ def minHarmonic2(sN, fs, t, f0, f0k1, f0kN, f0ks, f1, f1k1, f1kN, f1ks):
117 122
     # CG, BFGS, Newton-CG, L-BFGS-B, TNC, SLSQP, dogleg, trust-ncg, trust-krylov, trust-exact and trust-constr
118 123
     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 )
119 124
     print(res)
120
-    return harmonicEuler2(sN, fs, t, res.x[0], f0k1, f0kN, f0ks, res.x[1], f1kN, f1kN, f1ks)#[0]
125
+    return harmonicEuler2(sN, fs, t, res.x[0], f0k1, f0kN, f0ks, res.x[1], f1k1, f1kN, f1ks)#[0]
121 126
 
122 127
 def guessf0( sN, fs ):
123 128
     S = np.fft.fft(sN)

Завантаження…
Відмінити
Зберегти