社団法人 電子情報通信学会
THE INSTITUTE OF ELECTRONICS,
INFORMATION AND COMMUNICATION ENGINEERS

信学技報
TECHNICAL REPORT OF IEICE.

MATLAB と SPICE を併用した変調回路の周波数領域での解法
ボルテラ級数展開法と調波平衡法
河田 淳治†山上喜広††

西尾 芳文††

牛田 明夫†† 徳島文理大学 工学部
香川県さぬき市志度 1314-1†† 徳島大学 工学部
徳島県徳島市南常三島町 2-1
E-mail: †{kawata,ushida}@fe.bunri-u.ac.jp, ††{yamagami,nishio}@ee.tokushima-u.ac.jp
あらまし ボルテラ級数展開法は非線形回路の周波数領域での求解法であり，解が解析的な関数で与えられることから，RF 回路や通信回路の解析に広く用いられている。この手法は準線形理論に基づいており，弱非線形系に対してのみ有効である。しかしながら，その解の精度について定量的に議論した文献は見当たらない。一方，調波平衡法（HB法) は強非線形系にも適用できるが，一般的な電子回路を解析する場合 決定方程式が煩雑になるために 不向きであった。本稿では，HB 法に Spice と MATLAB を併用することによって，電子回路の周波数領域での解析を可能とした。
変調回路の解析に対してボルテラ級数展開法と HB 法を適用し，解の精度について議論すると共に，応用例として 差動型ミクサ回路を解析する。
キーワード ボルテラ級数展開法，調波平衡法，MATLAB，Spice，変調回路，ミクサ回路

Spice-oriented frequency-domain intermodulation analysis combining
with MATLAB
Comparisons between Volterra series and harmonic balance methods
Junji KAWATA† , Yoshihiro YAMAGAMI†† , Yoshifumi NISHIO†† , and Akio USHIDA†† Faculty of Eng., Tokushima Bunri University
1314-1 Shido, Sanuki, 769-2193 JAPAN†† Faculty of Eng., Tokushima University
2-1 Minami-josanjima, Tokushima, 770-8506 JAPANE-mail: †{kawata,ushida}@fe.bunri-u.ac.jp, ††{yamagami,nishio}@ee.tokushima-u.ac.jp
Abstract It is very important to analyze RF circuits, mixers and modulators driven by multiple input frequencies.
Volterra series methods are widely used for the frequency-domain analysis of nonlinear circuits, because they give thesolutions in analytical forms. On the other hand, HB (harmonic balance) method is well-known that can be appliedeven for relatively strong nonlinear circuits. In this article, we propose a Spice-oriented HB method combining with
MATLAB. Firstly, for the nonlinear devices modeled by the special functions, their characteristics are approximatedby the Taylor series, and the Fourier coefficients to the input and output relations can be calculated by MATLAB
in the symbolic forms. Thus, the determining equation of HB method can be formulated by the equivalent circuit
and/or net-list. It can be efficiently solved by the DC analysis of Spice. Thus, the frequency-domain solutions
such as frequency response curves can be easily obtained by the DC sweep. We found from the examples that,
although Volterra series method can be efficiently applied to weakly nonlinear circuits, it becomes erroneous forstrong nonlinear circuits. On the other hand, our HB method can be stably applied to relatively strong nonlinear
circuits.
Key words Volterra series method, harmonic balance method, MATLAB, Spice，modulator and mixer circuit—1—1. ま え が き
非線形電子回路の周波数領域での求解において，周波数応

id = IS exp vd /VT , IS = 10−12 [A], VT = 0.025.

(1)

まず最初に，vd を次のように仮定したとき
ただし Vd0 = 0.5822 （注 2）

vd = Vd0 + Vd cos ωt,

答曲線を求めることは IC や通信回路を設計する上で重要であ(2)る。ボルテラ級数展開法は解が解析的な関数で得られるためこ

応答 id を求めフーリエ級数に展開した直接法と，3 次のオーのような目的に適しており，欧米を中心に広く用いられている

ダのテーラ級数で近似したときの応答の DC, cos ωt, cos 2ωt,[1-3]。この解析手法は準線形理論に基づいている。この手法でcos 3ωt 成分を比較してみよう。ただしは，非線形特性がべき級数で与えられる必要があり，弱非線形

id = k0 + k1 vd + k2 vd2 + k3 vd3 ,回路の解析に適していると云われている。ボルテラ核の計算は

k0 = 0.0130, k1 = 0.5199, k2 = 10.39, k3 = 138.6

高次になると煩雑になるため，通常，3 次までのオーダが用い

振幅 Vd をパラメータとして変化させた場合の結果が表 1 に示されている。

られている（注 1）。一方，調波平衡法（HB) は古くから非線形振表1動の解析に広く用いられており，多くの興味ある非線形現象を(3)

3 次のテーラ級数と厳密解
3 次テーラ近似

厳密解

解明してきた [5-8]。 この解法は規模の大きい電子回路の解析にVd

対しては，°
1 決定方程式が煩雑になることと，°
2 トランジスタ

0.02

0.0151 0.0113 0.0022 0.0003

0.0151 0.0112 0.0021 0.0003などのモデルが特殊関数を含んでいるために解析的にフーリエ0.05

0.0295 0.0413 0.0179 0.0055

0.0259 0.0390 0.0130 0.0043

0.07

0.0538 0.0859 0.0469 0.0191

0.0383 0.0720 0.0255 0.0119

0.10

0.1460 0.2542 0.1624 0.0870

0.0647 0.1560 0.0520 0.0347

級数展開ができない等の理由で，一般的な電子回路の解析には

直流

基本

2倍

3倍

直流

基本

2倍

3倍

用いられなかった [9,10]。そのため反復解法に基づく数値的な手法が多く開発されている [11-14]。これらはある指定した入力に対して，正確な定常振動波形を求めることはできるが，周波数応答曲線のような大域的な特性解析には有効でない。そこで，
我々は古典的な HB 法が任意の非線形素子をもつ回路に適用で

表 1 より，入力振幅 Vd が小さい場合には誤差は小さいが，振幅の増加と共に急速に誤差が大きくなることが分かった。たとえば，Vd = 0.1 に対する 3 次近似の基本調波の誤差は 63 パー
セントである。そこで，精度を上げるために，テーラ近似の次数を 5 次，7 次と上げたときの結果が表 2 に示されている。
表2

きるように，フーリエ変換器なるものを Spice の ABM（アナテーラ級数の次数と厳密解 (Vd = 0.1)
直流成分

基本調波

2 倍調波

3 倍調波

解析を可能にした [15-16]。しかしながらこの方法では，2 つの厳密解

0.14590.25380.16700.08703 次近似

0.0646

0.15590.05190.0347異なる周波数成分を持つ回路に対しては，フーリエ変換器の構

5 次近似

0.11620.22500.1216

0.0696

成が煩雑となり，演算時間等の問題があることが分かった。そ7 次近似

0.13970.24830.15590.0832

ログビヘービア モデル）を用いて開発し，調波平衡法に基づく

こで，本稿では，非線形特性をテーラ級数で近似し，MATLAB

テーラ級数の次数を上げることにより厳密解に近づいており，を用いて入出力関係をシンボリックな形でフーリエ級数に展開

Vd = 0.1 に対しては 7 次まで考慮すれば十分であることが分し [17]，ネットリストとして調波成分ごとに Spice に取り込み，かる。

HB 法を用いて解析する手法を提案する。この HB 法では決定

2. 2 ボルテラ級数展開法

方程式がネットリストとして与えられ，Spice の直流解析で周ボルテラ級数展開法を理解するために，2 つの異なった周波

波数応答曲線を容易に得ることができるため，記号解析と同様

数成分をもつ入力で駆動される図 1 の回路について考えてみよな性質を持っている。この場合，トランジスタなどをあらかじう [3]。

めモジュール化し，ライブラリに登録しておくことによって，iG

容易に汎用化できる。
本解法の有効性を調べるために，ボルテラ級数法との精度のvin1(ω 1 t)

Linearvin2 (ω 2 t)

Network

比較をベース変調回路について行った。応用例として，差動型

図1のミクサ回路を解析した。いずれの回路に対しても周波数応答

+
vG
-

v(t)

2 つの入力を持つ非線形回路

曲線が安定に追跡できることが分かった。なお，上述の各種解

ただし，非線形特性は (3) の 3 次のべき級数で与えられるもの

法の特長を比較したものが，付録にまとめられている。

とする。ここで，複素周波数 “s” の領域でボルテラ核を求め

2. ボルテラ級数展開法と調波平衡法の比較ることにする（注 3）。1 次のボルテラ核は線形増分回路である図2(a) における “r = 1/k1 ” の出力電圧に等しい [3]。すなわち，

2. 1 非線形素子のテーラ級数近似

2 入力電源に対してはステップ入力 U1 (s), U2 (s) と置き，これ

ここでは，非線形素子をテーラ級数によって近似するとき，らの応答 H11 (s), H12 (s) を用いて 1 次のボルテラ核は次のよ

何次のオーダまで採用すれば良いかについて考える。バイポー

うに書かれる。

ラトランジスタは 2 個のダイオードと従属電流源からなるエH1k (s) = H11 (s) + H12 (s)(4)

バース・モルモデルによってよくモデリングされるため，次のような特性 (id , vd ) を持つダイオードについて議論する。

（注 2）：これはあとに出てくるベース変調回路の動作点電圧である。
（注 3）：ボルテラ級数展開は次数を幾らでも高くすることは可能であるが，次数（注 1）：ボルテラ級数展開法の結果は漸近的解法に等しいことが証明されているを高くするとボルテラ核が複雑化する。このため通常，3 次までが用いられてい

[4]

る。

—2—

ここで，H11 (s), H12 (s) は それぞれ入力 vin1 (ω1 t), vin2 (ωt)
に対応していることに注意されたい。

I2N (±jω1 , ±jω2 ) = k2 ((H11 (±jω1 ) + H12 (±jω2 ))× ((H11 (±jω1 ) + H12 (±jω2 ))2 次のボルテラ核は (3) と (4) から求められる 2 次のオーダの電流 I2N を求めることから始まる。すなわち，

I2N (s1 , s2 ) = k2 H1k (s1 )H1k (s2 )(12)= k2 [H11 (jω1 )H11 (jω1 ) + H11 (jω1 )H11 (−jω1 )+H12 (jω2 )H11 (jω2 ) + H12 (jω2 )H11 (−jω2 )(5)+2H11 (jω1 )H12 (jω2 ) + 2H11 (jω1 )H12 (−jω2 )]

このとき，線形抵抗 “r” から左側をみたアドミッタンスを Y (s)ここで，第 1 項が 2ω1 成分，第 3 項が 2ω2 成分，第 2，4 項がとすると，2 次のボルテラ核 H2k (s1 , s2 ) は直流成分，第 5 項が ω1 + ω2 成分，第 6 項が |ω1 − ω2 | 成分にY (s1 + s2 )H2k (s1 , s2 ) = −I2N (s1 , s2 )(6)の解として求まる。
3 次のボルテラ核も同様にして，3 次のオーダの電流を求め
I3N (s1 , s2 , s3 ) = k3 H1k (s1 )H1k (s2 )H1k (s3 ) + 32 k2 [H1k (s1 )×H2k (s2 , s3 ) + H1k (s2 )H2k (s1 , s3 ) + H1k (s3 )H2k (s1 , s2 )](7)

これを用いて 3 次のボルテラ核 H3k (s1 , s2 , s3 ) が次式から求められる。

Y (s1 + s2 + s3 )H3k (s1 , s2 , s3 ) = −I3N (s1 , s2 , s3 ).
U1(s)Linearr
Network+
H11(s)−

Linearr
Network(b) 2nd order(8)+
H2k(s1 ,s2 )−I 2N (s 1 ,s2 )対応している。この電流を (6) の右辺に代入することで，対応する周波数成分の核が求まる。3 次のボルテラ核についても同様である。このようにして，表 3 が得られる。これらの計算はMATLAB を用いることにより容易に計算される [17]。
2. 3 Spice による非線形電子回路の調波平衡法
ここでは，Spice と MATLAB を併用した調波平衡法（HB)の解析手法について述べる。まず，代表的な非線形素子の一つであるバイポーラトランジスタの HB モジュールを導出してみよう。エバース・モル モデルでは，ベース・エミッタ間電圧
VBE ，ベース・コレクタ間電圧 VBC によって，コレクタ電流
iC ，ベース電流 iB ，エミッタ電流 iE はそれぞれ次のように与
えられる。




iC = IS (αF exp(VBE /VT ) − exp(VBC /VT ))iB = IS ((1−αF ) exp(VBE /VT )+(1−αR ) exp(VBC /VT ))

iE = −IS (exp(VBE /VT ) + αR exp(VBC /VT ))Linearr
Network

U2(s)+
H12(s)−

LinearNetwork(c) 3rd order(a) 1st order
図2+
H3k(s 1 ,s 2 ,s3 )−I3N (s1 ,s2 ,s 3 )r

以下では次のような値を用いる。

IS = 10−12 [A], VT = 0.026[V], αF = 0.99, αR = 0.3
ここで，与えられた回路が 2 つの入力 (10) で駆動されるもの

ボルテラ核の計算とし，2 つの周波数の干渉によって生ずる周波数成分の内，次このようにして，複素周波数領域での解は次のように書ける。
Hout (s1 , s2 , s3 ) = H1k (s1 ) + H2k (s1 , s2 ) + H3k (s1 , s2 , s3 ). (9)のような 13 個の周波数成分を考慮する。

DC, ω1 , ω2 , 2ω1 , 2ω2 , 3ω1 , 3ω2
ω1 + ω2 , |ω1 − ω2 |, 2ω1 + ω2 , |2ω1 − ω2 |ここで，次の 2 入力に対する周波数領域での解を求めると，周
波数成分ごとに表 3 のようになる。

vin1 (t) = A1 sin ω1 t,
表3

vin2 (t) = A2 sin ω2 t

(13)

ω1 + 2ω2 , |ω1 − 2ω2 |







.

(14)

エバース・モル モデルは (13) のようにダイオード特性を利用(10)しているので，バイポーラトランジスタの HB モジュールはダイオードの HB モジュールを用いて構成することができる。そ周波数領域での出力電圧 [3]

周波数

振幅

こで，ダイオード HB モジュールを導出する。ダイオード特性

DC| 12 A21 H2k (jω1 , −jω1 ) + 21 A22 H2k (jω2 , −jω2 )|を形式的に動作点 Vd0 における 3 次のテーラ級数で近似する。

ω1|A1 H1k (jω1 ) + 32 A1 A22 H3k (jω1 , jω2 , −jω2 )id = k0 (Vd0 ) + k1 (Vd0 )vd + k2 (Vd0 )vd2 + k3 (Vd0 )vd3

+ 34 A31 H3k (jω1 , jω1 , −jω1 )|ω2|A2 H1k (jω2 ) + 32 A21 A2 H3k (jω1 , −jω1 , jω2 )式 (14) の干渉周波数に対して，ダイオードの入出力関係を次+ 34 A32 H3k (jω2 , jω2 , −jω2 )|のように仮定する。

ω1 + ω2|A1 A2 H2k (jω1 , jω2 )||ω1 − ω2 ||A1 A2 H2k (jω1 , −jω2 )|2ω1| 12 A21 H2k (jω1 , jω1 )|2ω2| 12 A22 H2k (jω2 , jω2 )|2ω1 + ω2
|2ω1 − ω2 |
ω1 + 2ω2|ω1 − 2ω2 |
3ω1
3ω2

(15)

vd (t) = Vd0 +k=1

id (t) = Id0 +| 34 A21 A2 H3k (jω1 , jω1 , jω2 )|| 34 A21 A2 H3k (jω1 , jω1 , −jω2 )|| 34 A1 A22 H3k (jω1 , jω2 , jω2 )|| 34 A1 A22 H3k (jω1 , −jω2 , −jω2 )|| 14 A31 H3k (jω1 , jω1 , jω1 )|| 14 A32 H3k (jω2 , jω2 , jω2 )|

K
X




{Id,2k−1 cos νk t+Id,2k sin νk t) 

, (16)

k=1

次に，式 (16) の第 1 式を (15) に代入し，フーリエ級数 id を

MATLAB によりシンボリック形式で求める。このとき生ずる
係数の掛け算の総数は

1 + 25 + 252 + 253 = 16276.

このとき，核の計算は

si = {±jω1 , ±jω2 }, i = 1, 2, 3(Vd,2k−1 cos νk t+Vd,2k sin νk t) 



K
X

(17)

と膨大な数となる。これらの中から (14) に相当する周波数成(11)分を選び出して Spice のネットリストを作成し，モジュール化と置くことによって求められる。ただし，式 (5) の計算には (4)すると図 3 が得られる。これを用いて，バイポーラトランジ

を用いてスタの HB モジュールを容易に作成することができる。このモ—3—

ジュールをライブラリに登録しておくと便利である。
C

Vd0
Vd1 cosν1 t
Vd2 sinν 1t

id

Diode HB

R

L

C0

Vout

Q

E

e (t)
1

model

e (t)

α F IdC0

R1

Id0
I d1 cosν1 t
I d2 sinν 1t

R2

V1 (jω 1 )Id

V2 (jω 2)

g

+

d

C

R

L Vout
-

2

+

-

vd

I d2K sinνKt

Vd2K sinνKt

(a)

(b)

図5

(a)
図3

(a) ベース変調回路， (b) １次のボルテラ級数の回路
R1 = 50[kΩ], R2 = 10[kΩ], R = 0.2[kΩ]C = 1[nF], C0 = 1[nF], C0 = 1[nF], L = 10[µF], E = 10[V]v1 = 0.02 sin ω1 t, v2 = 0.02 sin ω2 t, for ω2 = 8.2 × 106 [rad/sec](b)(a) ダイオードモデル，(b) ダイオードの HB モジュール

また，LCR 素子に対しては

2.5
[V]





iC = C vdtC ,L
vL = L idt,

iL = I1 cos νk t + I2 sin νk t

vR = RiR ,iR = IR1 cos νk t + IR2 sin νk t

vC = V1 cos νk t + V2 sin νk t

ω1 component(18)

IC1 = νk CV2 ,VL1 = νk LI2 ,VL2 = −νk LI1

VR1 = RIR1 ,VR2 = RIR2

|ω1 −ω2 |component

Cosine component
Inductor
Capacitor

vL

0.0




5

10

15

|2ω1 −ω2 |component

20
(c)

25
6
[x10 rad/s]

30(c) ボルテラ級数展開法によるベース変調回路特性

Sine component

Intermodulation phenomena

1.6

ω 2 component

1.2

ω1 component

I2

V L1

V L2

0.8

V L1 = νk LI 2

V L2 = −νkLI 1

0.4

2ω1 component|ω1 −ω2 |componentω1 +ω component
2

iC
0

vC

I C1

C

V1

V2

I C2

0

5

10

| 2ω1 −ω2 | component

15

20

(d)

I C1= νk CV2

Resistor

0

(19)
[V]

I1

L

2ω1 component ω1 +ω2 component

0.5となり，電流制御電圧源，電圧制御電流源を用いて実現される。

iL

ω 2 component

1.5

1.0

なる関係から Sine-Cosine 変換 [15,16] をすると


IC2 = −νk CV1 Intermodulation phenomena(Volterra)

2.0

iR

vR

R

R

(d) 調波平衡法によるベース変調回路特性

I C2 = − νk CV1

[V]
4.0

IR2

I R1
VR1

VR2

25
30
6
[x10 rad/s]

Steady-state waveformR
2.0

VR1 =RI R1

図4

VR2 =RI R2

0.0

LCR 素子の Sine-Cosine 変換-2.0

非線形素子のモジュールと共に上のような LCR の Sine-Cosine
変換によって構成される等価回路のネットリストは調波平衡法の決定方程式に等しく，Spice の DC 解析により解析される。

2. 4 ボルテラ級数展開法と調波平衡法の比較
ここでは，図 5(a) に示された簡単なベース変調回路を考え
よう [18]。バイポーラトランジスタにはエバース・モルモデル
を用い，ダイオードの特性は 3 次のテーラ級数で近似する。また，ボルテラ級数展開法の次数も 3 次のオーダとする。

2.4.1 ボルテラの場合

-4.0
10

(e) 定常振動波形

12

14

16

18

[µs]

20

ω1 = 107 [rad/sec], ω2 = 8.2 × 106 [rad/sec]

2.4.2 HB 法の場合 トランジスタをエバース・モルでモデリ
ングし，3 次のテーラ級数
iBC = IBC,0 (VBC,0 )+IBC,1 (VBC,0 )vBC
2
3
+IBC,2 (VBC,0 )vBC
+IBC,3 (VBC,0 )vBC

iBE = IBE,0 (VBE,0 )+IBE,1 (VBE,0 )vBE2
3
+IBE,2 (VBE,0 )vBE
+IBE,3 (VBE,0 )vBE







(22)

まず 動作点を求め，テーラ級数に展開して，この入出力関係を HB モジュール化し，これを用に 展 開 す る 必 要 が あ る 。図 5(a) の 回 路 に 対 し て は 動 作 点

いて同一パラメータに対する周波数応答を求めたものが図 5(d)vd0 = 0.5822[V] が得られ，このときの 3 次のテーラ級数はに示されている。全体的な特性はよく似ているが，HB 法の方

id = 0.0130 + 0.52vd + 10.40vd2 + 138.60vd3 .

が振幅が小さくなっていることと，ω2 成分はボルテラ法ではほ(20)

ぼ一定であるのに対して，HB 法では変化している。このとき（注 4）

このようにして，係数が次のように与えられる。

k1 = 0.52(= gd ), k2 = 10.40, k3 = 138.60(21)

2.3 節の手法を用いて 3 次までのボルテラ核を求め，表 3 の結果を利用し，MATLAB によって周波数領域での応答を求めると図 5(c) の周波数応答曲線が得られた。
（注 4）：この回路の場合，ベース・コレクタ間電圧は VBC < 0 であるから，

IS exp(VBC /VT ) = 0 とした。

の，演算時間は 78.38[sec] であった。次に，トランジスタに直接エバース・モル モデルを適用し，過渡解析により LCR タン
ク回路の電位差波形を計算したものが図 5(e) に示されている。
2.4.3 誤差比較 入力信号の振幅を変化させた場合の ω1 と ω2
成分を，ボルテラ級数展開法と HB 法により求めて誤差評価を行う（表 4）。ここで，ボルテラ法は直流動作点において，ダイオード特性を 3 次のテーラ級数で近似した。一方，HB 法では 全体の回路から求まった動作点を用いた。すなわち，この直—4—

流動作点の値は調波解析における他の調波成分にも関係することに注意しよう。過渡解析では，トランジスタモデルとしてエ
バース・モル モデルを用いている。

+
v

表4
ω1 と ω2 の出力振幅の比較(ω1 = 107 [rad/s], ω2 = 8.2 × 106 [rad/s])

HB 法

ボルテラ

V
out C

R

CC

L

Q

Q

1

2

過渡解

V1 = V2

ω1

ω2

ω1

ω2

ω1

ω2

0.01

0.950

0.693

0.929

0.691

0.913

0.698

0.02

2.20

1.41

1.72

1.31

1.68

1.29

0.03

3.91

2.32

2.50

1.92

2.17

1.72

0.04

6.25

3.35

3.35

2.59

2.50

2.11

0.05

9.80

4.55

4.31

3.34

2.75

2.38

v1 (t)Q3v 2 (t)

R2
V

V

1

図6

ここで，過渡解析の結果は (13) の特性を用いたので真値と見なすことができる。振幅の小さい V1 = V2 = 0.01 では，いずれの結果もほぼ等しいが，入力振幅が大きくなると，真値から

1[V]

ω 2 component

0.8

大きく異なっている。特に，ボルテラ級数展開法の誤差の増加0.6が目立つ。HB 法は真値に近いが，これも入力振幅の増加と共0.4に誤差が増加している。これらには 2 つの要因が考えられる。

0.2（ 1 ） 2.1 節で示したように，振幅が大きくなると，テーラ

ω component
1

Intermodulation phenomena[ω2 =57x106 rad/s]
ω1 +ω2 component[x20]|ω1 −ω2 | component[x20]03級数を高次まで採用する必要がある。しかし，ボルテラ級数展開法では 3 次以上は難しいので精度を上げることは期待でき

2(a) ミクサ回路 R2 = 10[kΩ], R = 100[Ω], L = 10[µH],
VCC = 10[V], V1 = −10[V], V2 = −5[V]

4

5

7

8

9
10
ω1 [x107 rad/s]

C = 1[pF], v1 = 0.002 sin ω1 t, v2 = 0.2 sin 57×106 t.

(b) 周波数応答
3.0
[V]

ない。

6
(b)

Output waveform

2.0（ 2 ） 解法の違いから生ずるもので，準線形法に基づくボル1.0

テラ法と，考慮する周波数成分を総合的に取り扱う HB 法の違0.0

いが考えられ，HB 法がより厳密解に近いのは当然である。

-1.0

3. 応 用 例-2.0

ここでは，3 つトランジスタから構成された図 6(a) の差動型ミクサ回路 [18] を解析してみよう。ただし，変調回路として使用する場合には， v1 (t) が搬送信号，v2 (t) が入力信号に相当
し，さらに LC 共振周波数は搬送信号に等しく選ばれる。さて，
トランジスタをエバース・モルモデル (13) を用いてモデリングを行い，
i d = IS e

λvd,

for IS = 10

−12[A], λ = 40, α = 0.99

(23)

形式的なテーラ級数で次のように近似する。
2

-3.0
20

21

22

24

23

(c)

[µsec]

25

v1 (t) = 0.002 sin 50 × 106 t, v2 (t) = 0.2 sin 57 × 106 t(c) 定常振動波形
4[V]Intermodulation phenomena6ω2=0.2[x10 rad/s]

3

ω 1 component2ω1 +ω2 component[x20]1| ω1−ω2 | component[x20]3id ' k0 (vd0 ) + k1 (vd0 )v + k2 (vd0 )v + k3 (vd0 )v , vd = vd0 + v(24)
ここで，式 (14) で示された，13 個の干渉周波数を考慮し，MATLAB を用いてシンボリック形式でフーリエ級数に展開し，ト
ランジスタの HB モジュール化を行う。このモジュールと LC
素子の Sine-Cosine 変換（図 4）を用いて，調波成分ごとに等価回路を作成する。このとき，これら部分回路は従属電源を用いて互いに結合されていることに注意しよう。これらは Spiceの直流解析によって解析される。図 6(b) に求めた周波数応答曲線を示す。ここで，ω2 = 57 × 106 [rad/s] に固定し，ω1 を連
続的に変化させた。次に ω1 = 50 × 106 [rad/s] における，出力波形を図 6(c) に示す。エバース・モル モデルを使用して過渡解析したものと，HB 法を比較したところ，次のような結果が得られた。
Vout,1 = 0.7919[V] (HB) V1 = 0.8311[V] (Tran. + DFT)
Vout,2 = 0.9603[V] (HB) V2 = 0.9585[V] (Trans. + DFT)
ω1 = 50 × 106 [rad/s], ω2 = 57 × 106 [rad/s]

両結果はほぼ一致している。一方，干渉波成分 (ω1 + ω2 ) と(|ω1 − ω2 |) は十分小さく，ω1 = ω = 57 × 106 [rad/s] のとこω2 component[x20]

0
4

6

8

10

12

14

(d)

16

18
20
ω1 [x10 6 rad/s](d) 変調回路の周波数応答
C = 1[nF], v1 = 0.001 sin ω1 t, v2 = 0.1 sin 0.2×106 t.
[V]
4.0
2.0
0.0

-2.0

-4.0
10

12

14
(e)

16

18

20
[µsec]

7(e) 定常波形 v1 (t) = 0.001 sin 10 t, v2 (t) = 0.1 sin 0.2 × 106 t

次に，図 6(a) を変調回路として考えることにする。このとき，RLC の共振周波数は 107 [rad/s] であり，搬送周波数 ω1 または，ω2 を 107 [rad/s] と設定する必要がある。他方は入力信号であり，低い周波数に設定される。この周波数応答が図 6(d)に示されている。過渡解析からの結果と我々の HB 法を比較す
ると次のようになる。

ろで，|ω1 − ω| 成分が零となっている。この場合の演算時間は102.05[s] であった。
—5—

Vout,1 = 3.987[V] (HB) V1 = 4.218[V] (Tran. + DFT)
Vout,2 = 0.307[V] (HB) V2 = 0.337[V] (Tran. + DFT)
ω1 = 107 [rad/s], ω2 = 0.2 × 106 [rad/s]




[10]

A.S.Sedra and K.C.Smith, Microelectronic Circuits, Oxford



[11]

R.Telichevesky, K.S.Kundert and J.K.White, “Efficient steady-Univ. Press, 2004.
state analysis based on matrix-free Krylov-subspace methods,”

いずれの結果もほぼ一致している（注 5）。計算時間は 209.30[s]ACM, pp.480-485, 1995.
[12]

であった。

K.S.Kundert, J.K.White and A.Sangiovanni-Vincentelle, SteadyState methods for simulating Analog and Microwave Circuits,

4. む す び

Kluwe Academic, Pub. 1990.

本稿では，ボルテラ級数展開法と Spice と MATLAB を併用[13]“Analysis of communication circuits based on multidimensional

した HB 法の比較を変調回路とミクサ回路について行った。そFourier transformation,” IEEE Trans. Comput.-Aided Des.

れらをまとめると次のようになる。
（ 1 ） ボルテラ級数展開法は記述関数的解析手法であるが，Interg. Circuits Syst., vol.18, pp.1165-1177, 1999.
[14]

HB 法もフーリエ級数展開が MATLAB によって行われ，このシンボリック形式を用いて，決定方程式が等価回路で与えられると云う意味で数値解析的と云うより記述関数法に近い。また，

A.Ushida,T.Adachi and L.O.Chua, “Steady-state response of
nonlinear circuitsbased on hybrid method,” IEEE Trans. on
Circuits and Systems-I, vol.39, pp.649-661, 1992.

[15]

A.Ushida, Y.Yamagami and Y.Nishio, “Frequency responses
of nonlinear networks using curve tracing algorithm,”ISCAS

いずれの方法も非線形特性はべき級数で表す必要がある。ただし，ボルテラの場合には動作点があらかじめ与えられる必要が

Y.Yamagami,Y.Nishio, A.Ushida, M.Takahashi and K.Ogawa,

2002, vol.I, pp.641-644, 2002.
[16]

J.Kawata, Y.Taniguchi, M.Oda,Y.Yamagami and Y.Nishio and
A.Ushida “ Spice-oriented frequencydomain analysis of nonlin-

あるが，HB 法では形式的なテーラ級数展開法を用いるため，

ear electronic circuits,” IEICE Trans. Fundamentals, vol.E90-

その必要はない。

A, pp.406-410, 2007.

（ 2 ） ボルテラ級数展開法は MATLAB を用いて，効率よく[17]

大石進一、MATLAB による数値計算、培風館、2004.

計算されるが，オーダが 3 次以上になると計算式が煩雑になり，[18]

K.K.Clarke and D.T.Hess, Communication Circuits: Analysis

これ以上のオーダで計算した文献は見当たらない。

[19]

R.Telichevesky and K.Kundert, SpectreRF Primer, Cadence（ 3 ） HB 法でのフーリエ級数への展開は MATLAB によってシンボリックに行われ，この結果をネットリストとしてライブラリに登録しておけば便利である。しかしながら，登録可能なフーリエ級数の項数は用いる Spice に依存するようである。
（ 4 ） 過渡解析（厳密解）との比較ではボルテラ級数展開法
よりも，HB 法の方が誤差の点では優れているようだ。

2.1 節に述べたように，非線形素子のテーラ級数展開の次数が解の精度に大きく影響しているようである。そこで，今後の課題として，高次オーダのテーラ級数にフーリエ級数展開を適用する HB 法を開発したい。
文
[1]

献

M.Schetzen, The Volterra and Wiener Theorems of Nonlinear
Systems, John Wiley and Sons, 1978.

[2]

J.Wood and D.E.Root, Fundamentals of Nonlinear Behavioral
Modeling for RF and Microwave Design, Artech House, 2005.

[3]
[4]

P.Wambacq and W.Sansen, Distortion Analysis of Analog In-付

録

1. 各種周波数領域での求解法の比較
ここでは，次の項目について比較してみる。

°
1 非線形素子特性の関数，°
2 非線形性の強さ，°
3 解析原理，°4
解の精度，°
5 多入力回路への対応，°
6 複数解が存在する場合の求解の可能性，°7 その他特長・欠点（ 1 ） ボルテラ級数展開法 [1-3]°
1 べき級数，°2 弱非線形系，°
3 記述関数，°4 強非線形系では誤差が大（3 次のオーダまで），°
5 簡単，°
6 不可能，°
7 MATLAB を用いて容易に解析できる。
（ 2 ） HB0 法 [15-16] °1 任意，°
2 非線形性が強くなると，等価回路の規模が大きくなる，°3 Spice，°
4 よい，°
5 煩雑化する，°
6 可能，°
7 非線形性が強くても解が安定に求められる。
（ 3 ） 本 HB 法 °1 べき級数，°2 非線形性が強くなると，等価回路の規模が大きくなる，°3 Spice，°
4 よい，°
5 簡単，°
6可

B.J.Leon and D.J.Shaefer, “Volterra series and Picard iteration

能，°
7 非線形性が強くても解が安定に求められる。

Systems, vol.25, pp.789-793, 1978.
C.Hayashi,

Nonlinear Oscillations in Physical Systems,

McGraw-Hill, 1964.
[6]

Y.Ueda, The Road to Chaos-II, Aerial Press. Inc., 2001.

[7]

R.J.Gilmore and M.B.Steer, “Nonlinear circuit analysis using（ 4 ） 数値解析的手法 [11-13] °1 任意，°
2 任意，°
3 FFT，ニュートン反復法，°
4 高精度，°
5 煩雑化し，演算効率悪い，°
6不可能，°
7 周波数応答曲線の求解には適さない。
（ 5 ） 緩和法 [13] °1 任意，°
2 任意，°
3 FFT，回路分割，緩the method of harmonic balance-A review of the Art. Part I. In-

和反復法，°
4 高精度，°
5 煩雑化し，演算効率悪い，°
6 不可能，

troductory concepts,” Int. Jour. of Microwave and Millimeter-°
7 周波数応答曲線の求解には適さない。

Wave Computer-Aided Eng. vol.1, pp.22-37, 1991.
[8]

Design Systems, San Jose, California, 1996.

tegrated Circuits, Kluwer Academic Pub., 1998.
for nonlinear circuits and systems,” IEEE Trans. Circuits and[5]and Design, Addison-Wesley Pub. Co., 1971.

R.J.Gilmore and M.B.Steer, “Nonlinear circuit analysis using（ 6 ） SpectreRF[19] °1 任意，°
2 任意，°
3 時間領域での解the method of harmonic balance-A review of the Art. Part II.

法，°
4 高精度，°
5 煩雑化し，演算効率悪い，°
6 不可能，°
7 周波

Advanced concepts,” Int. Jour. of Microwave and Millimeter-

数応答曲線の求解には適さない。

Wave Computer-Aided Eng. vol.1, pp.159-180, 1991.
[9]

牛田明夫、田中衛、電子回路シミュレーション、コロナ社、2003.

（注 5）：これらの誤差はトランジスタのテーラ級数近似から生じているものと考えられる。

—6—