コイル(1)

で、対象物のインピーダンスを$Z$とする。対象物に印可されている電圧を$V_D$、対象物を流れる電流を$I_D$とすると、 \[ V_D=V_1-V_2 \] \[ I_D=\frac{V_2}{50} \] なので、 \[ Z=50\left(\frac{V_1}{V_2}-1\right) \]

import matplotlib.pyplot as plt
import numpy as np
from scipy import signal

sig1=sig1-np.mean(sig1)
sig2=sig2-np.mean(sig2)
sig1_hil=signal.hilbert(sig1)
sig2_hil=signal.hilbert(sig2)
Z=50*(sig1_hil/sig2_hil-1)

nop=np.size(t)
analysis_start=int(nop*0.2)
analysis_stop=int(nop*0.8)

t=t[analysis_start:analysis_stop]
sig1=sig1[analysis_start:analysis_stop]
sig2=sig2[analysis_start:analysis_stop]
sig1_hil=sig1_hil[analysis_start:analysis_stop]
sig2_hil=sig2_hil[analysis_start:analysis_stop]
Z=Z[analysis_start:analysis_stop]
Z_re=np.mean(np.real(Z))
Z_im=np.mean(np.imag(Z))
print('Z=',end='\t')
print(Z_re,end='\t')
print(Z_im)

fig=plt.figure()
fig.add_subplot(3,1,1)
plt.plot(t,sig1,t,abs(sig1_hil))
fig.add_subplot(3,1,2)
plt.plot(t,sig2,t,abs(sig2_hil))
ax1=fig.add_subplot(3,1,3)
ax2=ax1.twinx()
ax1.plot(t,np.abs(Z),'C0')
ax2.plot(t,np.angle(Z),'C1')
plt.show()

こんなんで、インピーダンス計算を行う。sig1, sig2, tはオシロからデータを読み出したとする。
計算にあたって$V_1$と$V_2$は複素数じゃないといけないので、測定データをヒルベルト変換しています。またそのためにオシロのオフセットずれを補正しています。

\[ L=\frac{X_{meas}}{\omega} \]

import matplotlib
import matplotlib.pyplot as plt
import numpy as np

import ahkab
from ahkab import circuit

Cp=10e-12
Ls=4.7e-3
Rs=27
Rp=50e3

mycircuit = circuit.Circuit(title="Inductor Equivalent circuit")

gnd = mycircuit.get_ground_node()

mycircuit.add_capacitor("Cp", n1="n1", n2="n3", value=Cp)
mycircuit.add_inductor("Ls", n1="n1", n2="n2", value=Ls)
mycircuit.add_resistor("Rs", n1="n2", n2="n3", value=Rs)
mycircuit.add_resistor("Rp", n1="n1", n2="n3", value=Rp)
mycircuit.add_resistor("R0", n1="n3", n2=gnd, value=50)

mycircuit.add_vsource("V1", n1="n1", n2=gnd, dc_value=0, ac_value=1,)

print(mycircuit)

op_analysis = ahkab.new_op()
ac_analysis = ahkab.new_ac(start=1e3, stop=30e6, points=1001,sweep_type='LOG')

r = ahkab.run(mycircuit, an_list=[op_analysis, ac_analysis])

matplotlib.use('TkAgg') 

f=r['ac']['f']
V1=r['ac']['Vn1']
V2=r['ac']['Vn3']
Z=50*(V1/V2-1)

fig = plt.figure()
ax1=plt.subplot(1,1,1)
ax2=ax1.twinx()
ax1.semilogx(f,np.abs(Z),'C0')
ax2.semilogx(f,np.angle(Z)*180/np.pi,'C1')
plt.show()

for i in range(np.size(f)):
    print(f[i],end=(','))
    print(np.real(Z[i]),end=(','))
    print(np.imag(Z[i]))