-
Notifications
You must be signed in to change notification settings - Fork 0
/
tube_peak3.py
411 lines (346 loc) · 17.7 KB
/
tube_peak3.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
#coding:utf-8
# compute peak and drop-peak frequency detail of the tube
# by scipy.optimize.minimize_scalar
#
# This version is using frequency ratio.
import sys
import argparse
import numpy as np
from scipy import signal
from scipy import optimize
import matplotlib.pyplot as plt
# Check version
# Python 3.10.4, 64bit on Win32 (Windows 10)
# numpy 1.22.3
# matplotlib 3.5.2
# scipy 1.8.0
class compute_tube_peak(object):
def __init__(self, rg0=0.95, rl0=0.9 ,NUM_TUBE=3, disp=False):
self.rg0=rg0
self.rl0=0.9
self.C0=35000.0 # speed of sound in air, round 35000 cm/second
self.NUM_TUBE=NUM_TUBE
self.Delta_Freq=2.5 # 5
self.f_min=200
self.f_max=5000
self.f_out=100000 # 候補がないときに代入する値
self.f=np.arange(self.f_min, self.f_max, self.Delta_Freq)
self.xw= 2.0 * np.pi * self.f
self.sign0=1.0 # normal
self.disp=disp
self.counter=0
#
def __call__(self, X, xw_input=None):
# X[0,1]= L1,L2
# X[2,3]= A1,A2
# they should be same as get_ft5
if xw_input is not None:
xw= xw_input
else:
xw=self.xw
if (len(X) == 8) or (len(X) == 7) : # X=[L1,L2,L3,L4,A1,A2,A3,A4] or X=[L1,L2,L3,L4,r1,r2,r3] when four tube model
tu1= X[0] / self.C0 # delay time in 1st tube
tu2= X[1] / self.C0 # delay time in 2nd tube
tu3= X[2] / self.C0 # delay time in 3rd tube
tu4= X[3] / self.C0 # delay time in 4th tube
if len(X) == 8:
r1=( X[5] - X[4]) / ( X[5] + X[4]) # reflection coefficient between 1st tube and 2nd tube
r2=( X[6] - X[5]) / ( X[6] + X[5]) # reflection coefficient between 2nd tube and 3rd tube
r3=( X[7] - X[6]) / ( X[7] + X[6]) # reflection coefficient between 3rd tube and 4th tube
else:
r1=X[4]
r2=X[5]
r3=X[6]
func1= self.func_yb_t4
args1=(tu1,tu2,tu3,tu4,r1,r2,r3)
# abs(yi) = abs( const * (cos wv + j sin wv)) becomes constant. So, max/min(abs(val)) depends on only yb
self.yi= 0.5 * ( 1.0 + self.rg0 ) * ( 1.0 + r1) * ( 1.0 + r2) * ( 1.0 + r3) * ( 1.0 + self.rl0 ) * \
np.exp( -1.0j * ( tu1 + tu2 + tu3 + tu4 ) * xw)
# yb
yb1= 1.0 + r2 * r1 * np.exp( -2.0j * tu2 * xw ) # 2.3
yb1= yb1 + r3 * r2 * np.exp( -2.0j * tu3 * xw ) # 1.2
yb1= yb1 + self.rl0 * r3 * np.exp( -2.0j * tu4 * xw ) # 0
yb2= r3 * r1 * np.exp( -2.0j * ( tu2 + tu3) * xw ) # 2.2
yb2= yb2 + self.rl0 * r2 * np.exp( -2.0j * ( tu3 + tu4) * xw ) # 1.1
yb3= self.rl0 * r3 * r2 * r1 * np.exp( -2.0j * ( tu2 + tu4) * xw ) # 2.4
yb4= self.rl0 * r1 * np.exp( -2.0j * ( tu2 + tu3 + tu4) * xw ) # 2.1
yb31= self.rg0 * self.rl0 * np.exp( -2.0j * ( tu1 + tu2 + tu3 + tu4) * xw ) # 3.1
yb32= self.rg0 * r3 * np.exp( -2.0j * ( tu1 + tu2 + tu3 ) * xw ) # 3.2
yb41= self.rg0 * r2 * np.exp( -2.0j * ( tu1 + tu2 ) * xw ) # 4.1
yb42= self.rg0 * r2 * r3 * self.rl0 * np.exp( -2.0j * ( tu1 + tu2 + tu4) * xw ) # 4.2
yb51= self.rg0 * r1 * np.exp( -2.0j * tu1 * xw ) # 5.1
yb52= self.rg0 * r1 * r3 * self.rl0 * np.exp( -2.0j * ( tu1 + tu4 ) * xw ) # 5.2
yb53= self.rg0 * r1 * r2 * self.rl0 * np.exp( -2.0j * ( tu1 + tu3 + tu4 ) * xw ) # 5.3
yb54= self.rg0 * r1 * r2 * r3 * np.exp( -2.0j * ( tu1 + tu3 ) * xw ) # 5.4
self.yb= yb1 + yb2 + yb3 + yb4 + yb31 + yb32 + yb41 + yb42 + yb51 + yb52 + yb53 + yb54
elif (len(X) == 6) or (len(X) == 5) : # X=[L1,L2,L3,A1,A2,A3] or X=[L1,L2,L3,r1,r2] when three tube model
tu1= X[0] / self.C0 # delay time in 1st tube
tu2= X[1] / self.C0 # delay time in 2nd tube
tu3= X[2] / self.C0 # delay time in 3rd tube
if len(X) == 6:
r1=( X[4] - X[3]) / ( X[4] + X[3]) # reflection coefficient between 1st tube and 2nd tube
r2=( X[5] - X[4]) / ( X[5] + X[4]) # reflection coefficient between 2nd tube and 3rd tube
else:
r1=X[3]
r2=X[4]
func1= self.func_yb_t3
args1=(tu1,tu2,tu3,r1,r2)
# abs(yi) = abs( const * (cos wv + j sin wv)) becomes constant. So, max/min(abs(val)) depends on only yb
self.yi= 0.5 * ( 1.0 + self.rg0 ) * ( 1.0 + r1) * ( 1.0 + r2) * ( 1.0 + self.rl0 ) * \
np.exp( -1.0j * ( tu1 + tu2 + tu3 ) * xw)
# yb
yb1= 1.0 + r1 * self.rg0 * np.exp( -2.0j * tu1 * xw )
yb1= yb1 + r2 * r1 * np.exp( -2.0j * tu2 * xw )
yb1= yb1 + self.rl0 * r2 * np.exp( -2.0j * tu3 * xw )
yb2= r2 * self.rg0 * np.exp( -2.0j * ( tu1 + tu2) * xw )
yb2= yb2 + self.rl0 * r1 * np.exp( -2.0j * ( tu2 + tu3) * xw )
yb3= self.rl0 * r2 * r1 * self.rg0 * np.exp( -2.0j * ( tu1 + tu3) * xw )
yb4= self.rl0 * self.rg0 * np.exp( -2.0j * ( tu1 + tu2 + tu3) * xw )
self.yb= yb1 + yb2 + yb3 + yb4
elif (len(X) == 4) or (len(X) == 3): # else X=[L1,L2,A1,A2] or X=[L1,L2,r1] two tube model
tu1= X[0] / self.C0 # delay time in 1st tube
tu2= X[1] / self.C0 # delay time in 2nd tube
if len(X) == 4:
r1=( X[3] - X[2]) / ( X[3] + X[2]) # reflection coefficient between 1st tube and 2nd tube
else:
r1= X[2]
func1= self.func_yb_t2
args1=(tu1,tu2,r1)
# compute frequency response
# abs(yi) = abs( const * (cos wv + j sin wv)) becomes constant. So, max/min(abs(val)) depends on only yb
self.yi= 0.5 * ( 1.0 + self.rg0 ) * ( 1.0 + r1) * ( 1.0 + self.rl0 ) * \
np.exp( -1.0j * ( tu1 + tu2 ) * xw)
# yb
self.yb= 1.0 + r1 * self.rg0 * np.exp( -2.0j * tu1 * xw ) + \
self.rl0 * r1 * np.exp( -2.0j * tu2 * xw ) + \
self.rl0 * self.rg0 * np.exp( -2.0j * ( tu1 + tu2) * xw )
else:
print ('error: len(X) is not expected value.', len(X))
val= self.yi / self.yb
if xw_input is not None: # return response if there is xw_input
return np.sqrt(val.real ** 2 + val.imag ** 2)
else:
self.response=np.sqrt(val.real ** 2 + val.imag ** 2)
# get peak and drop-peak list
self.peaks_list=signal.argrelmax(self.response)[0] # signal.argrelmax output is triple
peaks= self.f[ self.peaks_list ]
self.drop_peaks_list=signal.argrelmin(self.response)[0]
drop_peaks= self.f[ self.drop_peaks_list ]
# 候補点がNUM_TUBEより少ないときは f_outを入れておく
if len(peaks) < self.NUM_TUBE:
peaks= np.concatenate( ( peaks, np.ones( self.NUM_TUBE - len(peaks)) * self.f_out ) )
elif len(peaks) > self.NUM_TUBE:
peaks=peaks[0: self.NUM_TUBE]
self.peaks_list=self.peaks_list[0: self.NUM_TUBE]
if len(drop_peaks) < self.NUM_TUBE:
drop_peaks= np.concatenate( ( drop_peaks, np.ones( self.NUM_TUBE - len(drop_peaks)) * self.f_out ) )
elif len(drop_peaks) > self.NUM_TUBE:
drop_peaks=drop_peaks[0: self.NUM_TUBE]
self.drop_peaks_list=self.drop_peaks_list[0: self.NUM_TUBE]
# より詳細に探索する
peaks_detail=np.zeros( len(peaks) )
drop_peaks_detail=np.zeros( len(drop_peaks) )
## peak
self.sign0=1.0 # normal
for l, xinit in enumerate( peaks ):
if xinit >= self.f_max:
peaks_detail[l]= xinit
else:
# Use brent method: 囲い込み戦略と二次近似を組み合わせ
b_xinit=[ xinit - self.Delta_Freq , xinit + self.Delta_Freq ]
res = optimize.minimize_scalar(func1, bracket=b_xinit, args=args1)
peaks_detail[l]= res.x
if self.disp:
print ('b_xinit', b_xinit)
print ('result x', res.x)
## drop-peak
self.sign0=-1.0 # turn upside down
for l, xinit in enumerate( drop_peaks ):
if xinit >= self.f_max:
drop_peaks_detail[l]= xinit
else:
b_xinit=[ xinit - self.Delta_Freq , xinit + self.Delta_Freq ]
res = optimize.minimize_scalar(func1, bracket=b_xinit, args=args1)
drop_peaks_detail[l]= res.x
if self.disp:
print ('b_xinit', b_xinit)
print ('result x', res.x)
return peaks_detail, drop_peaks_detail
def func_yb_t2(self, x, *args): # two tube *は可変長の引数
x = x
tu1,tu2,r1= args
xw= x * 2.0 * np.pi
yb= 1.0 + r1 * self.rg0 * np.exp( -2.0j * tu1 * xw ) + self.rl0 * r1 * np.exp( -2.0j * tu2 * xw ) + \
self.rl0 * self.rg0 * np.exp( -2.0j * ( tu1 + tu2) * xw )
return (yb.real**2 + yb.imag**2) * self.sign0
def func_yb_t3(self, x, *args): # three tube *は可変長の引数
tu1,tu2,tu3,r1,r2= args
xw= x * 2.0 * np.pi
yb1= 1.0 + r1 * self.rg0 * np.exp( -2.0j * tu1 * xw )
yb1= yb1 + r2 * r1 * np.exp( -2.0j * tu2 * xw )
yb1= yb1 + self.rl0 * r2 * np.exp( -2.0j * tu3 * xw )
yb2= r2 * self.rg0 * np.exp( -2.0j * ( tu1 + tu2) * xw )
yb2= yb2 + self.rl0 * r1 * np.exp( -2.0j * ( tu2 + tu3) * xw )
yb3= self.rl0 * r2 * r1 * self.rg0 * np.exp( -2.0j * ( tu1 + tu3) * xw )
yb4= self.rl0 * self.rg0 * np.exp( -2.0j * ( tu1 + tu2 + tu3) * xw )
yb= yb1 + yb2 + yb3 + yb4
return (yb.real**2 + yb.imag**2) * self.sign0
def func_yb_t4(self, x, *args): # four tube *は可変長の引数
tu1,tu2,tu3,tu4,r1,r2,r3= args
xw= x * 2.0 * np.pi
yb1= 1.0 + r2 * r1 * np.exp( -2.0j * tu2 * xw ) # 2.3
yb1= yb1 + r3 * r2 * np.exp( -2.0j * tu3 * xw ) # 1.2
yb1= yb1 + self.rl0 * r3 * np.exp( -2.0j * tu4 * xw ) # 0
yb2= r3 * r1 * np.exp( -2.0j * ( tu2 + tu3) * xw ) # 2.2
yb2= yb2 + self.rl0 * r2 * np.exp( -2.0j * ( tu3 + tu4) * xw ) # 1.1
yb3= self.rl0 * r3 * r2 * r1 * np.exp( -2.0j * ( tu2 + tu4) * xw ) # 2.4
yb4= self.rl0 * r1 * np.exp( -2.0j * ( tu2 + tu3 + tu4) * xw ) # 2.1
yb31= self.rg0 * self.rl0 * np.exp( -2.0j * ( tu1 + tu2 + tu3 + tu4) * xw ) # 3.1
yb32= self.rg0 * r3 * np.exp( -2.0j * ( tu1 + tu2 + tu3 ) * xw ) # 3.2
yb41= self.rg0 * r2 * np.exp( -2.0j * ( tu1 + tu2 ) * xw ) # 4.1
yb42= self.rg0 * r2 * r3 * self.rl0 * np.exp( -2.0j * ( tu1 + tu2 + tu4) * xw ) # 4.2
yb51= self.rg0 * r1 * np.exp( -2.0j * tu1 * xw ) # 5.1
yb52= self.rg0 * r1 * r3 * self.rl0 * np.exp( -2.0j * ( tu1 + tu4 ) * xw ) # 5.2
yb53= self.rg0 * r1 * r2 * self.rl0 * np.exp( -2.0j * ( tu1 + tu3 + tu4 ) * xw ) # 5.3
yb54= self.rg0 * r1 * r2 * r3 * np.exp( -2.0j * ( tu1 + tu3 ) * xw ) # 5.4
yb= yb1 + yb2 + yb3 + yb4 + yb31 + yb32 + yb41 + yb42 + yb51 + yb52 + yb53 + yb54
return (yb.real**2 + yb.imag**2) * self.sign0
def show_freq(self,):
# show rough(accuracy=Delta_Freq) result
fig = plt.figure()
ax1 = fig.add_subplot(211)
plt.title('frequency response')
plt.xlabel('Frequency [Hz]')
plt.ylabel('Amplitude [dB]')
if 1: # show peak and drop peak
ax1.semilogy(self.f, self.response, 'b', ms=2)
ax1.semilogy(self.f[self.peaks_list] , self.response[self.peaks_list], 'ro', ms=3)
ax1.semilogy(self.f[self.drop_peaks_list] , self.response[self.drop_peaks_list], 'co', ms=3)
if 1: # show yi and yb
ax1.plot( self.f , np.abs(self.yi), 'g')
ax1.plot( self.f , np.abs(self.yb), 'y')
plt.grid()
plt.axis('tight')
plt.show()
def reset_counter(self,):
self.counter=0
def cost_0(self, peaks2, drop_peaks2, peaks, drop_peaks, USE_COST_RATIO=False):
# lower cost function
if USE_COST_RATIO:
# Use only ratio portion, that is [1:]
# 周波数の比である2番目以降の要素を使って計算する。
if drop_peaks is None:
return abs(peaks[1:] - peaks2[1:]).mean()
else:
return (abs(peaks[1:] - peaks2[1:]).mean() + abs(drop_peaks[1:] - drop_peaks2[1:]).mean()) / 2.0
else:
if drop_peaks is None:
return abs(peaks - peaks2).mean()
else:
return (abs(peaks - peaks2).mean() + abs(drop_peaks - drop_peaks2).mean()) / 2.0
def calc_cost(self, X , peaks, drop_peaks, display_count=100, disp=False, OVER_VALUE=0.95):
# get mean of difference between target and new computed ones
peaks2, drop_peaks2= self.__call__(X)
cost0= self.cost_0( peaks2, drop_peaks2, peaks, drop_peaks, USE_COST_RATIO=False )
# add penalty if reflection coefficient abs is over than OVER_VALUE
if len(X) == 3 and abs( X[2]) > OVER_VALUE:
cost0 += 1000.0
elif len(X) == 5 and ( abs( X[3]) > OVER_VALUE or abs( X[4] ) > OVER_VALUE ):
cost0 += 1000.0
elif len(X) == 7 and ( abs( X[4]) > OVER_VALUE or abs( X[5] ) > OVER_VALUE or abs( X[6] ) > OVER_VALUE ):
cost0 += 1000.0
if disp :
print (X,cost0, peaks2, drop_peaks2)
self.counter +=1
# show present counter value, don't show if display_count is negative
if display_count > 0 and self.counter % display_count == 0:
sys.stdout.write("\r%d" % self.counter)
sys.stdout.flush()
return cost0
# helper functions
def get_r1( X ):
return ( X[1] - X[0]) / ( X[1] + X[0]) # return reflection coefficient between 1st tube and 2nd tube
def get_A2( r1, A1 ):
if abs(r1) >= 1.0:
print ('error: abs(r1) > 1.0')
return (( 1.0 + r1) / ( 1 - r1)) * A1 # return cross-section area of 2nd tube
def get_A1( r1, A2 ):
if abs(r1) >= 1.0:
print ('error: abs(r1) > 1.0')
return (( 1.0 - r1) / ( 1 + r1)) * A2 # return cross-section area of 1st tube
def get_A1A2( r1, A_min=1.0):
# return cross-section area A1 and A2 under the condition of
# minimum cross-section is fixed as A_min
if r1 >= 0.0:
return A_min, get_A2(r1, A_min)
else:
return get_A1(r1, A_min), A_min
def get_A1A2A3( r1, r2, A_min=1.0):
# return cross-section area A1 A2 and A3 under the condition of
# minimum cross-section is fixed as A_min
A1=1.0
A2=get_A2( r1, A1 )
A3=get_A2( r2, A2 )
min_index= np.argmin( [A1,A2,A3] )
if min_index == 0:
A1= A_min
A2= get_A2( r1, A1 )
A3= get_A2( r2, A2 )
elif min_index == 1:
A2= A_min
A1= get_A1( r1, A2)
A3= get_A2( r2, A2)
elif min_index == 2:
A3= A_min
A2= get_A1( r2, A3)
A1= get_A1( r1, A2)
return A1, A2, A3
def get_A1A2A3A4( r1, r2, r3, A_min=1.0):
# return cross-section area A1 A2 A3 and A4 under the condition of
# minimum cross-section is fixed as A_min
A1=1.0
A2=get_A2( r1, A1 )
A3=get_A2( r2, A2 )
A4=get_A2( r3, A3 )
min_index= np.argmin( [A1,A2,A3,A4] )
if min_index == 0:
A1= A_min
A2= get_A2( r1, A1 )
A3= get_A2( r2, A2 )
A4= get_A2( r3, A3 )
elif min_index == 1:
A2= A_min
A1= get_A1( r1, A2)
A3= get_A2( r2, A2)
A4= get_A2( r3, A3)
elif min_index == 2:
A3= A_min
A4= get_A2( r3, A3)
A2= get_A1( r2, A3)
A1= get_A1( r1, A2)
elif min_index == 3:
A4= A_min
A3= get_A1( r3, A4)
A2= get_A1( r2, A3)
A1= get_A1( r1, A2)
return A1, A2, A3, A4
if __name__ == '__main__':
# shape symmetry check
# 全長が同じ長さの２管モデルは2種類存在するが、その多くは ピークとドロップピークの位置がかなり近くなる対称性がある。
# 大きく差がでるケースもある。完全な対称性はない。
# X[ L1, L2, r1] and X2[ L2, L1, r1] are same total lenght LT=L1+L2 (tu1+tu2).
# instance
tube = compute_tube_peak(rg0=0.95, rl0=0.9)
#LA_ranges = ( slice(0.5,13,0.5), slice(0.5,13,0.5), slice(-0.9, 0.9, 0.1) ) # specify X value range
total_lenght=26.0
L1= np.arange( 0.5,total_lenght,0.5)
L2= total_lenght - L1
r1= -0.8
for l1 in L1:
l2= total_lenght - l1
X_init= [ l1, l2, r1]
X_init2=[ l2, l1, r1]
peaks, drop_peaks= tube(X_init)
peaks2, drop_peaks2= tube(X_init2)
cost0= tube.cost_0( peaks2, drop_peaks2, peaks, drop_peaks)
print ( cost0)