ここでは,シングル出力段の設計に役立つRのプログラムを紹介します.
プレート電圧のベクトルと,グリッド電圧のベクトルから, ある演算を施して行列を作成するのに, Rの組み込み関数 outer が使えますが, この関数は2つの引数をとる関数に適用するものであるため, 次のようにして使う必要があります.
Ep <- seq(0, 500, by=1) # プレート電圧:1Vきざみ Eg <- seq(0, -120, by=-10) # グリッド電圧:10Vきざみ ip <- outer(Ep, Eg, function(x, y) Ip(t2A3, x, y))outer の第3の引数には,2つの引数をとる関数を指定します. この関数には,outer の第1の引数と,outer の第2の引数が, 最終的な行列の要素を網羅するような組み合わせとなって渡されます. outer は第3の引数の関数が返してきた値を行列の形に整えて返します.
このように第3の引数として引数の順序を入れ換える関数をいちいち指定するのは面倒なので, touter という関数を作成してあります.
1 touter <- 2 function(func, p, Ep, Eg, ...) 3 { 4 # プレート電圧とグリッド電圧の組み合わせについて、関数を適用する 5 # @param func 適用する関数(Ip など) 6 # @param p 真空管のパラメータ 7 # @param Ep プレート電圧 8 # @param Eg プレート電流 9 # @param ... func に渡す追加の引数 10 # @return Ep x Eg の行列 11 12 z <- outer(Ep, Eg, function(x, y) func(p, x, y, ...)) 13 dimnames(z) <- list(Ep, Eg) 14 z 15 }
先ほどの例を touter を使って求めてみます.
ip <- touter(Ip, t2A3, Ep, Eg)5極管で,一定のスクリーングリッド電圧に対してプレート電流を求める場合, 次のようにします.
Eg2 <- 250 ip <- touter(Ipp, t6L6T, Ep, Eg, Eg2)touter の第5の引数は
...
となっていますが,
これは任意の個数の引数を表しており,
リストの11行目にあるように第一の引数 func で指定された関数の
4個目以降の引数として与えられます.
この例で touter の第一の引数 func として指定された Ipp は,
4つの引数をとり,4つ目の引数はスクリーングリッド電圧ですので,
touter の5番目の引数として Eg2 を指定すると,
Ipp には4番目の引数として渡されます.
ウルトラリニア接続のプレート電流を求めるには, プレート電圧に応じてスクリーングリッド電圧を変化させなければなりません. 例えば,次のようにします.
Ep0 <- 250 # 静止時のプレート電圧 Eg2 <- 250 # 静止時のスクリーングリッド電圧 r <- 0.43 # スクリーングリッドタップ比率 ip <- touter(function(p, ep, eg) { eg2 <- Eg2 + (ep - Ep0) * r Ipp(p, ep, eg, eg2) }, t6L6T, Ep, Eg)touter に指定する関数として,引数が3つの一時的な無名の関数を作成しています. その関数の中でプレート電圧に対応したスクリーングリッド電圧を計算し, Ipp でプレート電流を求めています.
プレート特性を描くには matplot を用いてもよいのですが, 一般に使われないプレート特性図の右上の領域を削除したグラフを描く, g.plate という関数を用意しました.
1 g.plate <- 2 function(Ip, Ebb, RL, Ipmax, drawEg=TRUE, col=2, Ipmin=0, xlim, 3 xname="Ep", yname="Ip", zname="Eg") 4 { 5 # プレート特性のグラフを描く 6 # @param Ip touter で作成したプレート電流の行列(Ep x Eg) 7 # @param Ebb グラフの右上を消すための仮想的なB電圧 8 # @param RL グラフの右上を消すための仮想的な負荷抵抗 9 # @param Ipmax y軸の最大値 10 # @param drawEg グリッド電圧を表示するかどうか 11 # @param col グラフの色 12 # @param Ipmin y軸の最小値 13 # @param xlim x軸の範囲 14 # @param xname x軸の名前(単位は V, kV で自動的に付加) 15 # @param yname y軸の名前(単位は mA, A で自動的に付加) 16 # @param zname グリッド電圧の名前(単位は V で自動的に付加) 17 18 Ep <- as.numeric(dimnames(Ip)[[1]]) 19 Eg <- as.numeric(dimnames(Ip)[[2]]) 20 21 # グラフの右上のデータを消す 22 ilim <- matrix((Ebb - Ep)/RL, nrow=length(Ep), ncol=length(Eg)) 23 Ip[Ip > ilim] <- NA 24 25 # グラフを描く 26 if (missing(xlim)) 27 xlim <- range(Ep) 28 ylim <- c(Ipmin, Ipmax) 29 xlab <- eunit1(pretty(xlim, 8), "V") 30 ylab <- eunit1(pretty(ylim), "A") 31 matplot(Ep, Ip, type="l", col=col, lty=1, ylim=ylim, 32 xaxs="i", yaxs="i", 33 xaxt="n", yaxt="n", 34 xlab=paste(xname, xlab$unit), 35 ylab=paste(yname, ylab$unit)) 36 axis(1, at=xlab$at, lab=xlab$lab) 37 axis(1, at=xlab$at, lab=FALSE, tck=1, lty=2, lwd=0.3) 38 axis(2, at=ylab$at, lab=ylab$lab) 39 axis(2, at=ylab$at, lab=FALSE, tck=1, lty=2, lwd=0.3) 40 41 # Eg を書く 42 if (drawEg) { 43 maxi <- apply(Ip, 2, function(i) { 44 i[i > Ipmax] <- NA 45 max(seq(along=i)[is.finite(i)]) 46 }) 47 x <- pmin(Ep[maxi], (max(Ep) - min(Ep))* 0.97 + min(Ep)) 48 y <- pmin(Ip[cbind(maxi, 1:ncol(Ip))], (Ipmax - Ipmin) * 0.97 + Ipmin) 49 lab <- as.character(Eg) 50 lab[1] <- paste0(zname, "=", lab[1], "V") 51 text(x, y, lab) 52 } 53 invisible() 54 }ip には,touter などで作成した,プレート電流の, プレート電圧 x グリッド電圧 の行列を指定します. Ebb と RL がグラフの右上の部分を消すためのパラメータで, x軸上でプレート電圧が Ebb の点を通り, 傾きが -1/ の直線の右上の部分が表示されないようになります. グラフを描く範囲を狭めたい時は,Ipmax, Ipmin, xlim などを 使います. Ipmax は必ず指定しなければなりません. これ以降の引数はオプションであり, 必要に応じて指定します.
例として2A3のプレート特性図を描いてみましょう.
Ep <- seq(0, 500, by=1) # プレート電圧:1Vきざみ Eg <- seq(0, -120, by=-10) # グリッド電圧:10Vきざみ ip <- touter(Ip, t2A3, Ep, Eg) g.plate(ip, 600, 600/0.2, 0.15)図2.50のような特性図が描かれます.
6L6のプレート特性図を描いてみましょう.
Ep <- seq(0, 500, by=1) # プレート電圧:1Vきざみ Eg <- seq(0, -40, by=-5) # グリッド電圧:5Vきざみ Eg2 <- 250 ip <- touter(Ipp, t6L6T, Ep, Eg, Eg2) g.plate(ip, 600, 600/0.25, 0.2)同様にして,スクリーングリッド電流を求め, グラフに追加します.
ig2 <- touter(Ig2, t6L6T, Ep, Eg, Eg2) matlines(Ep, ig2, col="blue", lty=1)図2.51のような特性図が描かれます.
6L6のUL接続のプレート特性図を描いてみましょう.
Ep <- seq(0, 500, by=1) # プレート電圧:1Vきざみ Eg <- seq(0, -50, by=-5) # グリッド電圧:5Vきざみ Ep0 <- Eg2 <- 250 r <- 0.43 ip <- touter(function(p, ep, eg) { eg2 <- Eg2 + (ep - Ep0) * r Ipp(p, ep, eg, eg2) }, t6L6T, Ep, Eg) g.plate(ip, 600, 600/0.2, 0.15)図2.52のような特性図が描かれます.
1 trans.se <- 2 function(p, ei, Ep0, Eg0, RL, Eg20, rUL, rKF, MaxEp, Zp, Zk, Rbs) 3 { 4 # シングル出力段の伝達特性 5 # @param p 真空管パラメータ 6 # @param ei 入力電圧 7 # @param Ep0 静止時プレート電圧 8 # @param Eg0 バイアス電圧 9 # @param RL (P-K間)負荷インピーダンス 10 # @param Eg20 スクリーングリッド電圧 11 # @param rUL ULタップ比率(P-B間を1とする) 12 # @param rKF KNFタップ比率(P-B間を1とする) 13 # @param MaxEp 探索する最大プレート電圧 14 # @param Zp P-B間インピーダンス 15 # @param Zk K-E間インピーダンス 16 # @param Rbs ブートストラップ負荷抵抗 17 # @return 18 # $ep プレート電圧(対カソード) 19 # $ip プレート電流 20 # $ig グリッド電流 21 # $ek カソード電圧(対グラウンド) 22 # $ig2 スクリーン電流 23 # $bUL P-Kに対するULタップ比率 24 # $bKF P-Kに対するKNFタップ比率 25 # $ibs ブートストラップによる出力電流(P-K間相当) 26 27 triode <- missing(rUL) && missing(Eg20) 28 if (missing(RL) && missing(Zp)) 29 stop("RL or Zp must be specified.") 30 if (missing(RL)) { 31 if (!missing(Zk)) { 32 rKF <- sqrt(Zk / Zp) 33 } else if (missing(rKF)) { 34 rKF <- 0 35 } 36 RL <- Zp * (1 + rKF)^2 37 } else { 38 if (missing(rKF)) 39 rKF <- 0 40 } 41 bKF <- rKF / (1 + rKF) 42 43 # 動作点のプレート電流を計算し、交点を求めるための関数を準備する 44 if (triode) { 45 ip0 <- Ip(p, Ep0, Eg0) # 動作点のプレート電流 46 if (missing(Rbs)) { 47 f <- function(ep) { 48 dEp <- Ep0 - ep # プレート電圧の変化(*(-1)) 49 ek <- dEp * bKF # カソード電位(グラウンド基準) 50 ip <- Ip(p, ep, eii-ek) # プレート特性によるプレート電流 51 irl <- ip0 + dEp/RL # ロードラインによるプレート電流 52 ip - irl # 0になったところが交点 53 } 54 } else { 55 f <- function(ep) { 56 dEp <- Ep0 - ep # プレート電圧の変化(*(-1)) 57 ek <- dEp * bKF # カソード電位(グラウンド基準) 58 ip <- Ip(p, ep, eii-ek) + # プレート特性によるプレート電流 59 (eii-ek-Eg0) / Rbs * bKF 60 irl <- ip0 + dEp/RL # ロードラインによるプレート電流 61 ip - irl # 0になったところが交点 62 } 63 } 64 } else { 65 if (missing(Eg20)) 66 Eg20 <- Ep0 67 if (missing(rUL)) 68 rUL <- 0 69 bUL <- (rUL + rKF) / (1 + rKF) # P-Kに対するUL比率 70 z <- Ipp.sub(p, Ep0, Eg0, Eg20) 71 ip0 <- z$ip + z$ig2 * bUL # P-K 間に換算した信号電流 72 if (missing(Rbs)) { 73 f <- function(ep) { 74 dEp <- Ep0 - ep # プレート電圧の変化(*(-1)) 75 ek <- dEp * bKF # カソード電位(グラウンド基準) 76 eg2 <- Eg20 - dEp * bUL # スクリーングリッド電圧(カソード基準) 77 z <- Ipp.sub(p, ep, eii-ek, eg2) 78 ip <- z$ip + z$ig2 * bUL # プレート特性による信号電流 79 irl <- ip0 + dEp/RL # ロードラインによる信号電流 80 ip - irl # 0になったところが交点 81 } 82 } else { 83 f <- function(ep) { 84 dEp <- Ep0 - ep # プレート電圧の変化(*(-1)) 85 ek <- dEp * bKF # カソード電位(グラウンド基準) 86 eg2 <- Eg20 - dEp * bUL # スクリーングリッド電圧(カソード基準) 87 z <- Ipp.sub(p, ep, eii-ek, eg2) 88 ip <- z$ip + z$ig2 * bUL + # プレート特性による信号電流 89 (eii-ek-Eg0) / Rbs * bKF 90 irl <- ip0 + dEp/RL # ロードラインによる信号電流 91 ip - irl # 0になったところが交点 92 } 93 } 94 } 95 96 if (missing(MaxEp)) 97 MaxEp <- (Ep0 + ip0 * RL) * 1.1 # 探索範囲の既定値 98 99 # 入力信号に対するプレート電圧、プレート電流を求める 100 ig <- ek <- ep <- ip <- rep(NA, length(ei)) 101 ibs <- rep(0, length(ei)) 102 if (!triode) 103 ig2 <- ip 104 for (i in seq(along=ei)) { 105 eii <- Eg0 + ei[i] # 入力電圧を加えたグリッド電圧 106 ep[i] <- uniroot(f, c(1, MaxEp))$root # 入力電圧に対するプレート電圧 107 dEp <- Ep0 - ep[i] 108 ek[i] <- dEp * bKF # カソード電位(グラウンド基準) 109 if (triode) { 110 z <- Ip.sub(p, ep[i], eii-ek[i]) 111 ip[i] <- z$ip 112 ig[i] <- z$ig 113 } else { 114 eg2 <- Eg20 - dEp * bUL # スクリーングリッド電圧(カソード基準) 115 z <- Ipp.sub(p, ep[i], eii-ek[i], eg2) 116 ip[i] <- z$ip 117 ig[i] <- z$ig 118 ig2[i] <- z$ig2 119 } 120 if (!missing(Rbs)) 121 ibs[i] <- (eii-ek[i]-Eg0) / Rbs * bKF 122 } 123 124 # 結果を返す 125 if (triode) 126 list(ep=ep, ip=ip, ig=ig, ek=ek, ibs=ibs, bKF=bKF) 127 else 128 list(ep=ep, ip=ip, ig=ig, ek=ek, ig2=ig2, ibs=ibs, bUL=bUL, bKF=bKF) 129 }p には,真空管のパラメータを指定します. ei には,入力電圧(信号成分)を与えます. Ep0 には,静止時のプレート電圧を指定します. Eg0 には,静止時のグリッド電圧を指定します. これら2つの値は,自己バイアスであってもカソードを基準とした電圧で指定します. RL には,負荷となるインピーダンスを指定します.
これ以降の引数はオプションであり, 必要に応じて指定します. Eg2 は,5極管のスクリーングリッド電圧です. これが指定されていない場合は,3極管接続として扱われます. rUL は,スクリーングリッドのタップ比率です. これを指定して Eg2 を指定しないばあいは, Eg2 の値として Ep0 が使われます. rKF は,カソード帰還のタップ比率です. 5極管の場合,スクリーンが電源に直結している場合は, rUL にも同じ値を指定してください. カソード帰還とULを併用している場合は, > となるはずです. また,rKF を使わずに,Zp と Zk で 一次インピーダンスとカソード帰還巻線のインピーダンスを指定することもできます. カソード帰還をかける場合の RL を計算するのは面倒なので, こちらの指定のほうが簡単です. Rbs, ibs については,5.5節で説明します.
例として,2A3のシングルで,フルスイングした場合の ロードラインと伝達特性を描きます.
Ep0 <- 250 Eg0 <- -43.5 ei <- seq(-Eg0, Eg0, len=21) # 加えうる最大の振幅を20分割 RL <- 2500 # 負荷インピーダンス tr <- trans.se(t2A3, ei, Ep0, Eg0, RL) points(tr$ep, tr$ip) # プレート特性図上に表示 lines(tr$ep, tr$ip, col="blue") # ロードライン plot(ei+Eg0, tr$ip, type="l", col="red", ylim=c(0, 0.15), yaxs="i") # 伝達特性のグラフ
シングルの場合は,ロードラインが直線になるので, throughline を使って,
Ip0 <- Ip(t2A3, Ep0, Eg0) throughline(Ep0, Ip0, -1/RL, col="blue")のように引いても構いませんが, 伝達特性を求めれば, ロードライン上のどの部分が使われているか,明確にわかります.
正弦波を入力したときの出力波形とその歪率を求めます. gen.sig という関数で正弦波を生成します.
1 gen.sig <- 2 function(a=1, n=128) 3 { 4 # 正弦波を生成する. 5 # @param a 振幅(尖頭値) 6 # @param n 1周期の分割数 7 # @return 各時点の値 8 # $uniq 重複する点を除いた値 9 # $idx 重複する点を除いた値を各時点の値に戻すための添字 10 # $sin 正弦波の値 11 12 if (n %% 4 != 0) 13 stop("n must be a multiple of 4.") 14 n2 <- n / 2 15 n4 <- n / 4 16 half <- a * sin((-n4:n4) / n * 2 * pi) 17 idx <- c((n4+1):(n2+1), n2:2, 1:n4) 18 list(uniq=half, idx=idx, sin=half[idx]) 19 }a は,振幅で,波高値で指定します. n は,1周期の分割数で,デフォルトは128となっています. 伝達特性の計算時間がかかる場合は,64や32に減らすとよいのですが, その分精度が悪くなります.
ロードラインが直線の場合,正弦波の1周期には同じ電圧の点が2箇所あり, すべての点について伝達特性を求めるのは無駄です. このため,gen.sin が返す値はリストとなっていて, $sin には正弦波の波形が、 $uniq には重複した値を除いたベクトルが, $idx には元の波形に戻すための添え字が格納されています.
> sig <- gen.sig(-Eg0) > sig$sin [1] 0.000000 2.134444 4.263746 6.382776 8.486429 10.569638 [7] 12.627383 14.654709 16.646729 18.598647 ... ... [127] -4.263746 -2.134444 > sig$idx [1] 33 34 35 36 37 38 39 40 41 ... ... [126] 30 31 32 > sig$uniq [1] -43.500000 -43.447602 -43.290536 -43.029178 ... ... [61] 42.664160 43.029178 43.290536 43.447602 43.500000 > plot(sig$sin) # 正弦波が描かれる > plot(sig$uniq) # 正弦波の一意の部分が描かれる > plot(sig$uniq[sig$idx]) # 元通りの正弦波が描かれる
2A3シングルの例で,出力波形を求めます.
tr <- trans.se(t2A3, sig$uniq, Ep0, Eg0, RL) ip <- tr$ip[sig$idx] # プレート電流の波形 ep <- tr$ep[sig$idx] # プレート電圧の波形 x <- 0:length(ip) # x軸の値 plot(x, c(ip, ip[1]), type="l", col="red") # プレート電流の波形
FFTを使って歪率を求めるには,distortion を使います.
1 distortion <- 2 function(wave, component=FALSE) 3 { 4 # FFTにより歪率を求める 5 # @param wave 波形 6 # @param component 歪み波形を求めるか? 7 # @return 8 # $wave: 与えられた波形 9 # $power: パワースペクトル 10 # $distortion: 歪率 11 # $wave.component: 歪み波形 12 13 n <- length(wave) 14 wave.fft <- fft(wave)/n 15 16 # パワースペクトル 17 power.c <- wave.fft[1:(n/2)] * c(1, rep(2, n/2-1)) 18 power <- abs(power.c) 19 names(power) <- c("DC", paste("H", 1:(n/2-1), sep="")) 20 21 # 歪率 22 dist <- c(sqrt(sum(power[-c(1,2)]^2)), power[-c(1,2)]) / power[2] 23 names(dist) <- c("THD", paste("H", 2:(n/2-1), sep="")) 24 25 # 歪成分(波形) 26 if (component) { 27 wave.comp <- matrix(NA, n, 4) 28 wave.comp[, 1] <- Re(fft(wave.fft * c(0, 1, rep(0, n-3), 1), inv=TRUE)) 29 wave.comp[, 2] <- Re(fft(wave.fft * c(0, 0, 1, rep(0, n-5), 1, 0), inv=TRUE)) 30 wave.comp[, 3] <- Re(fft(wave.fft * c(0, 0, 0, 1, rep(0, n-7), 1, 0, 0), inv=TRUE)) 31 wave.fft[c(1, 2, n)] <- 0 # DC と基本波を消す 32 wave.comp[, 4] <- Re(fft(wave.fft, inv=TRUE)) 33 wave.comp <- rbind(wave.comp, wave.comp[1, ]) 34 dimnames(wave.comp) <- list(NULL, c(paste("H", 1:3, sep=""), "Total")) 35 } else 36 wave.comp <- NULL 37 list(wave=wave, power=power, power.c=power.c, distortion=dist, 38 wave.component=wave.comp) 39 }wave には,プレート電流,プレート電圧などの波形を与えます. component を TRUE にすると, 3次までの高調波の波形を求めます. 歪率は,返値の distortion 成分に, 各高調波成分は,返値の wave.component 成分に格納されています.
2A3 シングルの歪率を求め, 歪み成分の波形を描いてみます.
> d <- distortion(ip, comp=TRUE) > d$dist * 100 # 歪率を%で表示 THD H2 H3 H4 H5 H6 5.369067e+00 5.363715e+00 9.312125e-03 2.384032e-01 1.095285e-02 1.969539e-02 ... > plot(2:10, d$dist[2:10]*100, type="h", col="red", xlab="Harmonics", ylab="Distortion (%)") > matplot(x, d$wave.comp, type="l", col=2:5, lty=1)全高調波歪率は5.4%であり, グラフから,3次高調波はほとんど存在しないことがわかります. ayumi