終端短絡線路と終端開放線路の入力インピーダンス周波数特性を計算してみると、 共振点近くで急激に変化し、 その中間点では周波数による変化が少ないことがわかりますので、 open/short 法を共振点の中央、 すなわち 1/8 波長の奇数倍の周波数で使えば誤差が減って、 ある程度高い周波数でも使えるのではないかというアイデアが生まれます。
例えば、次ぎのような測定手順になります。
これは1/8 波長 open/short 法というアイデアの実現には共振法 が必要になるということでもありますから、 共振法から簡単かつ高精度で得られる、 速度係数を位相定数の検証用 あるいは位相定数測定値として使えるという利点が生まれます。
以下、この方針で書いたシミュレーションプログラムの一例ですが、 かなりopen/short 法の欠点をカバーできることがわかります。 どの程度の改善が得られるかは 1/8 波長Open-Short法シミュレーション で実験できるようにしてあります。
プログラムへの誤差要因の追加法は すでに説明しました。
# 1/8 波長 open/short 法
/usr/bin/awk '
BEGIN {
_Z0 = 75.0
_a = 1.373e-6
_b = 8.385e-12
_v = 0.66
_R = 0.0914
_C = 69e-12
l = 1
PI = atan2(0, -1)
Vr = 1/sqrt(2.3)
# 終端短絡試料
# 1/4 波長共振周波数探索
f0s = find_f0s(l, Vr) # 4.94193e+07
print "f0s", f0s >"/dev/stderr"
# n/4 (n = 1, 2, ..) wavelength frequency
K = 0
F0s[K] = f0s # f0
K++
for (i = 1; i < 25; i++) {
f = f0s * (i + 1)
if (1e9 < f)
break
f = find_next(f - f0s * 0.20, f + f0s * 0.20)
F0s[K] = f
K++
}
# 1/8 波長の奇数倍の周波数で入力インピーダンスを測定
N = prev = 0
for (i = 0; i < K; i++) {
Fs[N] = prev + (F0s[i] - prev) / 2
prev = F0s[i]
meas_s(Fs[N])
Xs[N] = x; Ys[N] = y
N++
}
# 終端開放試料
# 1/4 波長共振周波数探索
f0o = find_f0o(l, Vr) # 4.94193e+07
print "f0o", f0o >"/dev/stderr"
K = 0
F0o[K] = f0o # f0
K++
for (i = 1; i < 25; i++) {
f = f0o * (i + 1)
if (1e9 < f)
break
f = find_next(f - f0o * 0.20, f + f0o * 0.20)
F0o[K] = f
K++
}
# Vr
df = (F0o[int(K / 2)] - F0o[0]) / int(K / 2)
print "df =", df >"/dev/stderr"
vr = 2 * 2 * df * l / 2.99792458e8 # velocity ratio
print "vr =", vr >"/dev/stderr"
# 1/8 波長の奇数倍の周波数で入力インピーダンスを測定
N = prev = 0
for (i = 0; i < K; i++) {
Fo[N] = prev + (F0o[i] - prev) / 2
prev = F0o[i]
meas_o(Fo[N])
Xo[N] = x; Yo[N] = y
N++
}
# Z0
R0 = K = 0
for (i = 0; i < N; i++) {
Mul(Xs[i], Ys[i], Xo[i], Yo[i])
Sqrt(x, y)
R0 += x
K++
}
print "Z0 =", R0 / K >"//dev/stderr"
# 伝搬定数 (α + j*β)
offset = b0 = T = K = 0
for (i = 0; i < N; i++) {
f = sqrt(Fs[i] * Fo[i])
if (i == 0)
f0 = f
Mul(Xs[i], Ys[i], Xo[i], Yo[i])
Sqrt(x, y)
T += x # Real(Z0)
Div(Xs[i], Ys[i], Xo[i], Yo[i])
Sqrt(x, y)
Atanh(x, y)
b = y
if (b < 0)
b += PI
if (b < b0)
offset += PI
print f, x, b + offset # 周波数、減衰定数、位相定数
b0 = b
K++
}
z0 = T / K
print "Vr =", Vr, "C =", C, "Z0 =", Z0 >"/dev/stderr"
# 減衰定数を決定
print "x Freq. (Hz)"
print "y Att. (neper/m)"
printf("s Z0 = %.1f Ohm, Vr = %.3f\n", Z0, Vr)
# Alpha
for (i = 0; i < K; i++) {
Alpha[i] = atanh(Alpha[i] / Z0)
print Fr[i], Alpha[i]
}
# 減衰定数回帰式 Aloha(f) = b1 * sqrt(f) の b1 を求める
for (i = 0; i < K; i++) {
Y[i] = Alpha[i]
X[i] = sqrt(Fr[i])
}
reg1(K) # 回帰直線の係数を求める
print "b0 =", b0, "b1 =", b1 >"/dev/stderr"
print "Alpfa(f) =", b1, "*sqrt(f)" >"/dev/stderr"
# exit 0
# 回帰式 Aloha(f) = b0 + b1 * sqrt(f) + b2 がほしい場合
for (i = 0; i < K; i++) {
Y[i] = Alpha[i]
X1[i] = sqrt(Fr[i])
X2[i] = Fr[i]
}
reg2(K)
print "b0 =", b0, "b1 =", b1, "b2 =", b2 >"/dev/stderr"
print "Alpfa(f) =", b0, "+", b1, "* sqrt(f) +", b2, "* f" >"/dev/stderr"
exit 0
}
END {
exit 0
}
# キャパシタンス決定
function cap(Vr, M, i, f1, bl, cm, a, b, c, d, p, y1, y2) {
# 速度係数の決定
# (β*l)/tan(β*l) 補正したキャパシタンスが一定になる β を求める
# 黄金分割法で回帰直線の分散の極小値を見付ける
a = 0.45 # 0.45 <= vr <= 1 の範囲を探索
b = 1.00
p = (sqrt(5) - 1) * 0.5
c = b - p * (b - a)
d = a + p * (b - a)
y1 = cap2(c, M)
y2 = cap2(d, M)
for (i = 0; i < 100; i++) {
if ((d - c) / (c + d) < 1e-3)
break
if (y1 < y2) {
b = d
d = c
y2 = y1
c = b - p * (b - a)
y1 = cap2(c, M)
}
else if (y1 > y2) {
a = c
c = d
y1 = y2
d = a + p * (b - a)
y2 = cap2(d, M)
}
}
vr = (c + d) * 0.5
# 低域の速度係数は高域よりわずかに小さい
# 単位長あたりキャパシタンス
a = b = 0
for (i = 0; i < M; i++) {
w = 2 * PI * fm[i]
cm = -xm[i] / w / (rm[i] * rm[i] + xm[i] * xm[i])
bl = w / vr / 2.99792458e8 * l
c = cm / sin(bl) * cos(bl) * bl
a += c
b++
}
C = a / b
return C
}
# 速度係数 vr を仮定したキャパシタンスの周波数特性
function cap2(vr, M, i, w, c, cm, bl) {
for (i = 0; i < M; i++) {
w = 2 * PI * fm[i]
cm = -xm[i] / w / (rm[i] * rm[i] + xm[i] * xm[i])
bl = w / vr / 2.99792458e8 * l
c = cm / sin(bl) * cos(bl) * bl
X[i] = i * 1e6
Y[i] = c
}
ve = reg1(M)
return ve # 分散
}
# 回帰直線を求める
# Y[i] = b0 + b1 * X[i] (i = 1, 2, .. N)
function reg1(N, \
i, x1, y1, x2, y2, tx, ty, qx, qy, sx, sy, cxy, a, b, ve) {
tx = ty = qx = qy = cxy = offset = 0
for (i = 0; i < N; i++) {
tx += X[i]
ty += Y[i]
qx += X[i] * X[i]
qy += Y[i] * Y[i]
cxy += X[i] * Y[i]
}
sx = qx - tx * tx / N
sy = qy - ty * ty / N
cxy -= tx * ty / N
b1 = cxy / sx
b0 = ty / N - a * (tx / N)
ve = (sy - cxy * cxy / sx) / (N - 2)
return ve # 分散
}
# ２変数の重回帰分析
# y = b0 + b1 * X1[] + b2 * X2[]
function reg2(N, T1, T2, Ty, S11, S22, S12, Syy, S1y, S2y) {
T1 = T2 = Ty = S11 = S22 = S12 = Syy = S1y = S2y = 0
for (i = 0; i < N; i++) {
T1 += X1[i]
T2 += X2[i]
Ty += Y[i]
S11 += X1[i] * X1[i]
S22 += X2[i] * X2[i]
S12 += X1[i] * X2[i]
Syy += Y[i] * Y[i]
S1y += X1[i] * Y[i]
S2y += X2[i] * Y[i]
}
Ty /= N; T1 /= N; T2 /= N # 平均
# T1 = T2 = Ty = 0 # y = b1 * X1 + b2 * X2 (b0 = 0) model
S11 -= N * T1 * T1 # 変動
S22 -= N * T2 * T2
S12 -= N * T1 * T2 # 共変動
Syy -= N * Ty * Ty
S1y -= N * T1 * Ty
S2y -= N * T2 * Ty
b1 = (S1y * S22 - S2y * S12) / (S11 * S22 - S12 * S12)
b2 = (S2y * S11 - S1y * S12) / (S11 * S22 - S12 * S12)
b0 = Ty - b1 * T1 - b2 * T2
return
}
function meas_o(f, a, b, w) {
w = 2 * PI * f
a = _a * sqrt(f) + _b * f
b = w / 2.99792458e8 * sqrt(2.3) + a
Zoc(_Z0, a, b, l * 0.98)
}
function meas_s(f, a, b, w) {
w = 2 * PI * f
a = _a * sqrt(f) + _b * f
b = w / 2.99792458e8 * sqrt(2.3) + a
Zsc(_Z0, a, b, l)
}
function find_next(f1, f2, i, f, w, a, b) {
for (i = 0; i < 50; i++) {
f = (f1 + f2) * 0.5
meas_o(f)
if (y == 0.0 || (f2 - f1) / f < 1e-6)
break
else if (y < 0)
f1 = f
else
f2 = f
}
return f
}
# open circuit 1/4 wave length frequency
function find_f0o(l, Vr, i, f, f1, f2, w, a, b)
{
f0 = 2.99792458e8 * Vr / l * 0.25
f1 = f0 * 0.75; f2 = f0 * 1.25
for (i = 0; i < 50; i++) {
f = (f1 + f2) * 0.5
meas_o(f)
if (y == 0.0 || (f2 - f1) / f < 1e-6)
break
else if (y < 0)
f1 = f
else
f2 = f
}
return f
}
# short circuit 1/4 wave length frequency
function find_f0s(l, Vr, i, f, f1, f2, w, a, b)
{
f0 = 2.99792458e8 * Vr / l * 0.25
f1 = f0 * 0.75; f2 = f0 * 1.25
for (i = 0; i < 50; i++) {
f = (f1 + f2) * 0.5
w = 2 * PI * f
a = _a * sqrt(f) + _b * f
b = w / (2.99792458e8 * Vr) + a
Zsc(_Z0, a, b, l)
if (y == 0.0 || (f2 - f1) / (f1 + f2) < 1e-12)
break
else if (y > 0)
f1 = f
else
f2 = f
}
return f
}
# Zin = (Zt/Z0 + tanh(γ*l) / (1 + (Zt/Z0)*tanh(γ*l)) * Z0
function Zin(z0, x1, y1, l, x2, y2, x3, y3, x4, y4) {
Tanh(x1 * l, y1 * l)
x3 = x
y3 = y
Mul(x, y, x2 / z0, y2 / z0)
x4 = x + 1
y4 = y
Add(x3, y3, x2 / z0, y2 / z0)
Div(x, y, x4, y4)
x *= z0
y *= z0
}
# Zsc = Z0 * tanh((α + j * β) * l)
function Zsc(z0, x1, y1, l) {
Tanh(x1 * l, y1 * l)
Mul(x, y, z0, 0)
}
# Zoc = Z0 * coth((α + j * β) * l)
function Zoc(z0, x1, y1, l, w) {
Tanh(x1 * l, y1 * l)
Div(z0, 0, x, y)
}
function atanh(x) {
return 0.5 * log((1 + x) / (1 - x))
}
function acoth(x) {
return 0.5 * log((x + 1) / (x - 1))
}
# complex functions
function Add(x1, y1, x2, y2) {
x = x1 + x2
y = y1 + y2
}
function Sub(x1, y1, x2, y2) {
x = x1 - x2
y = y1 - y2
}
function Mul(x1, y1, x2, y2) {
x = x1 * x2 - y1 * y2
y = x1 * y2 + y1 * x2
}
function Div(x1, y1, x2, y2, r2) {
r2 = x2 * x2 + y2 * y2
x = (x1 * x2 + y1 * y2) / r2
y = (y1 * x2 - x1 * y2) / r2
}
function Exp(x1, y1, a) {
a = exp(x1)
x = a * cos(y1)
y = a * sin(y1)
}
function Log(x1, y1) {
x = 0.5 * log(x1 * x1 + y1 * y1)
y = atan2(y1, x1)
}
function Sqrt(x1, y1, r, w) {
x = 0.5 * log(x1 * x1 + y1 * y1) * 0.5
y = atan2(y1, x1) * 0.5
a = exp(x)
x = a * cos(y)
y = a * sin(y)
}
function Pow(x1, y1, x2, y2, x3, y3) {
Log(x1, y1)
x3 = x
y3 = y
Mul(x3, y3, x2, y2)
Exp(x, y)
}
function Sin(x1, y1, e, f) {
e = exp(y1)
f = 1 / e
x = 0.5 * cos(x1) * (e - f)
y = 0.5 * sin(x1) * (e + f)
}
function Cos(x1, y1, e, f) {
e = exp(y1)
f = 1 / e
x = 0.5 * cos(x1) * (f + e)
y = 0.5 * sin(x1) * (f - e)
}
function Tan(x1, y1, e, f, d) {
e = exp(2 * y1)
f = 1 / e
d = cos(2 * x1) + 0.5 * (e + f)
x = sin(2 * x1) / d
y = 0.5 * (e - f) / d
}
function Sinh(x1, y1, e, f) {
e = exp(x1)
f = 1 / e
x = 0.5 * (e - f) * cos(y1)
y = 0.5 * (e + f) * sin(y1)
}
function Cosh(x1, y1, e, f) {
e = exp(x1)
f = 1 / e
x = 0.5 * (e + f) * cos(y1)
y = 0.5 * (e - f) * sin(y1)
}
function Tanh(x1, y1, e, f, d) {
e = exp(2 * x1)
f = 1 / e
d = 0.5 * (e + f) + cos(2 * y1)
if (d == 0) {
e = exp(x1)
f = 1 / e
x = 0
y = (e + f) * sin(y1)
}
else {
x = 0.5 * (e - f) / d
y = sin(2 * y1) / d
}
}
function Atanh(x1, y1) {
Div(1 + x1, y1, 1 - x1, -y1)
Log(x, y)
x *= 0.5
y *= 0.5
}
function Abs(x, y) {
return sqrt(x * x + y * y)
}
function Arg(x, y) {
return atan2(y, x)
}
'
平林浩一, 2016-04