Subsections

2.6 Rのプログラム

ここでは,シングル出力段の設計に役立つRのプログラムを紹介します.

2.6.1 プレート特性を描く

プレート特性を描くには,プレート電圧とグリッド電圧を変化させて, プレート電流を求めておきます. Rの場合,グラフ(散布図)を描くときには,列ごとに線が描かれるので, プレート電流の値を プレート電圧 x グリッド電圧 の 形の行列にすると便利です.

プレート電圧のベクトルと,グリッド電圧のベクトルから, ある演算を施して行列を作成するのに, 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 という関数を作成してあります.

2.6.1.0.1 プレート電圧,グリッド電圧を組み合わせた値を求める 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 という関数を用意しました.

2.6.1.0.2 プレート特性図を描く 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 グリッド電圧 の行列を指定します. EbbRL がグラフの右上の部分を消すためのパラメータで, x軸上でプレート電圧が Ebb の点を通り, 傾きが -1/$ \tt RL$ の直線の右上の部分が表示されないようになります. グラフを描く範囲を狭めたい時は,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のような特性図が描かれます.
図 2.50: 2A3のプレート特性図を描く
\includegraphics{figs/2A3_plate.ps}

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のような特性図が描かれます.
図 2.51: 6L6のプレート特性図を描く
\includegraphics{figs/6L6_plate.ps}

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のような特性図が描かれます.
図 2.52: 6L6のUL接続プレート特性図を描く
\includegraphics{figs/6L6_UL_plate.ps}

2.6.2 ロードラインと伝達特性

2.6.2.0.1 シングル出力段の伝達特性 trans.se

    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を併用している場合は, $ \tt rUL$ > $ \tt rKF$ となるはずです. また,rKF を使わずに,ZpZk で 一次インピーダンスとカソード帰還巻線のインピーダンスを指定することもできます. カソード帰還をかける場合の 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")    # 伝達特性のグラフ
図 2.53: 2A3のロードラインと伝達特性
\includegraphics{figs/2A3_ll_trans.ps}

シングルの場合は,ロードラインが直線になるので, throughline を使って,

Ip0 <- Ip(t2A3, Ep0, Eg0)
throughline(Ep0, Ip0, -1/RL, col="blue")
のように引いても構いませんが, 伝達特性を求めれば, ロードライン上のどの部分が使われているか,明確にわかります.

2.6.3 出力波形

正弦波を入力したときの出力波形とその歪率を求めます. gen.sig という関数で正弦波を生成します.

2.6.3.0.1 正弦波を生成する 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])        # 元通りの正弦波が描かれる
図 2.54: gen.sig により生成された正弦波
\includegraphics[width=\linewidth]{figs/gen_sig.ps}

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")  # プレート電流の波形
図 2.55: プレート電流の波形を描く
\includegraphics{figs/2A3_wf.ps}

FFTを使って歪率を求めるには,distortion を使います.

2.6.3.0.2 歪率を求める 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 には,プレート電流,プレート電圧などの波形を与えます. componentTRUE にすると, 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次高調波はほとんど存在しないことがわかります.
図 2.56: 歪率と歪み成分のグラフを描く
\includegraphics{figs/2A3_dist.ps}
ayumi
2016-03-07