|
|
@@ -106,10 +106,21 @@ def harmonicEuler ( f0, sN, fs, nK, t ):
|
|
106
|
106
|
nK = number of harmonics to calculate
|
|
107
|
107
|
t = time samples
|
|
108
|
108
|
"""
|
|
109
|
|
- print("building Euler matrix ")
|
|
110
|
|
- A = np.zeros( (len(t), nK), dtype=np.complex64)
|
|
111
|
|
- for irow, tt in enumerate(t):
|
|
112
|
|
- A[irow,:] = np.exp(1j* np.arange(1,nK+1) * 2*np.pi* (f0/fs) * irow)
|
|
|
109
|
+
|
|
|
110
|
+ #print("building Euler matrix ")
|
|
|
111
|
+ #A = np.zeros( (len(t), nK), dtype=np.complex64)
|
|
|
112
|
+ #for irow, tt in enumerate(t):
|
|
|
113
|
+ # A[irow,:] = np.exp(1j* np.arange(1,nK+1) * 2*np.pi* (f0/fs) * irow)
|
|
|
114
|
+
|
|
|
115
|
+ #AA = np.zeros( (len(t), nK), dtype=np.complex64)
|
|
|
116
|
+ A = np.exp(1j* np.tile( np.arange(1,nK+1),(len(t), 1)) * 2*np.pi* (f0/fs) * np.tile(np.arange(len(t)),(nK,1)).T )
|
|
|
117
|
+ #AA = np.tile(np.arange(len(t)), (nK,1)).T
|
|
|
118
|
+
|
|
|
119
|
+ #plt.matshow( np.imag(A), aspect='auto' )
|
|
|
120
|
+ #plt.matshow( np.real(AA), aspect='auto' )
|
|
|
121
|
+ #plt.show()
|
|
|
122
|
+ #exit()
|
|
|
123
|
+ #print ("A norm", np.linalg.norm(A - AA))
|
|
113
|
124
|
|
|
114
|
125
|
v = np.linalg.lstsq(A, sN, rcond=None) # rcond=None) #, rcond=1e-8)
|
|
115
|
126
|
alpha = np.real(v[0]) #[0::2]
|
|
|
@@ -175,7 +186,7 @@ def minHarmonic(f0, sN, fs, nK, t):
|
|
175
|
186
|
f02 = guessf0(sN, fs)
|
|
176
|
187
|
print("minHarmonic", f0, fs, nK, " guess=", f02)
|
|
177
|
188
|
# CG, BFGS, Newton-CG, L-BFGS-B, TNC, SLSQP, dogleg, trust-ncg, trust-krylov, trust-exact and trust-constr
|
|
178
|
|
- res = minimize( harmonicNorm, np.array((f0)), args=(sN, fs, nK, t)) #, method='CG', jac=jacEuler) #, hess=None, bounds=None )
|
|
|
189
|
+ res = minimize(harmonicNorm, np.array((f0)), args=(sN, fs, nK, t)) #, method='CG', jac=jacEuler) #, hess=None, bounds=None )
|
|
179
|
190
|
print(res)
|
|
180
|
191
|
return harmonicEuler(res.x[0], sN, fs, nK, t)
|
|
181
|
192
|
|
|
|
@@ -183,7 +194,7 @@ def minHarmonic2(f1, f2, sN, fs, nK, t):
|
|
183
|
194
|
#f02 = guessf0(sN, fs)
|
|
184
|
195
|
#print("minHarmonic2", f0, fs, nK, " guess=", f02)
|
|
185
|
196
|
#methods with bounds, L-BFGS-B, TNC, SLSQP
|
|
186
|
|
- 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='TNC' )
|
|
|
197
|
+ res = minimize( harmonic2Norm, np.array((f1,f2)), args=(sN, fs, nK, t), jac='2-point') #, bounds=((f1-1.,f1+1.0),(f2-1.0,f2+1.0)), method='TNC' )
|
|
187
|
198
|
print(res)
|
|
188
|
199
|
return harmonicEuler2(res.x[0], res.x[1], sN, fs, nK, t)
|
|
189
|
200
|
|
|
|
@@ -202,7 +213,7 @@ if __name__ == "__main__":
|
|
202
|
213
|
import matplotlib.pyplot as plt
|
|
203
|
214
|
|
|
204
|
215
|
f0 = 60 # Hz
|
|
205
|
|
- f1 = 62 # Hz
|
|
|
216
|
+ f1 = 61 # Hz
|
|
206
|
217
|
delta = np.random.rand()
|
|
207
|
218
|
delta2 = np.random.rand()
|
|
208
|
219
|
print("delta", delta)
|
|
|
@@ -211,7 +222,7 @@ if __name__ == "__main__":
|
|
211
|
222
|
phi = np.random.rand()
|
|
212
|
223
|
phi2 = np.random.rand()
|
|
213
|
224
|
A = 1.0
|
|
214
|
|
- A2 = 0.25
|
|
|
225
|
+ A2 = 0.0
|
|
215
|
226
|
A3 = 1.0
|
|
216
|
227
|
nK = 35
|
|
217
|
228
|
T2 = .200
|
|
|
@@ -229,10 +240,10 @@ if __name__ == "__main__":
|
|
229
|
240
|
|
|
230
|
241
|
# single freq
|
|
231
|
242
|
#h = harmonicEuler( f0, sN, fs, nK, t)
|
|
232
|
|
- #h = minHarmonic( f0, sN, fs, nK, t)
|
|
|
243
|
+ h = minHarmonic( f0, sN, fs, nK, t)
|
|
233
|
244
|
|
|
234
|
245
|
# two freqs
|
|
235
|
|
- h = minHarmonic2( f0, f1, sN, fs, nK, t)
|
|
|
246
|
+ #h = minHarmonic2( f0, f1, sN, fs, nK, t)
|
|
236
|
247
|
#h = harmonic2( f0, f1, sN, fs, nK, t)
|
|
237
|
248
|
#h = harmonicEuler2( f0, f1, sN, fs, nK, t)
|
|
238
|
249
|
|
