# Comparison Analysis Based on the Cubic Spline Wavelet and Daubechies Wavelet of Harmonic Balance Method

Authors: Jing Gao
Journal: Abstract and Applied Analysis
 (2014)

---

## Abstract

This paper develops a theoretical analysis of harmonic balance method, based on the cubic spline wavelet and Daubechies wavelet, for steady state analysis of nonlinear circuits under periodic excitation. The properties of the resulting Jacobian matrix for harmonic balance are analyzed. Numerical experiments illustrate the theoretical analysis.


---

## Body

## 1. Introduction

The rapid growth in integrated circuits has placed new demands on the simulation tools. Many quantities properties of circuits are of interest to circuits designer. Especially, the steady-state analysis of nonlinear circuits represents one of the most computationally challenging problems in microwave design.Harmonic balance (HB) [1–3] is very favorable for periodic or quasiperiodic steady-state analysis of mildly nonlinear circuits using Fourier series expansion. However, the density of the resulting Jacobian matrix seriously affects the efficiency of HB based on Fourier series expansion. More effective simulators are required to study steady-state analysis. Soveiko and Nakhla in [4] have provided the elaborate formulation for HB approach applying Daubechies wavelet series instead of Fourier series and obtain the sparser Jacobian matrix to reduce the whole computational cost. And Steer and Christoffersen in [5] have given the possibility of wavelet expansion for steady-state analysis. One advantage of wavelet bases is a sparse representation matrix of operators or functions which is favorable for solving the nonlinear system by Newton iterative method. But the main disadvantage is the waste of much time in storing Jacobian matrix due to the complex computation of Daubechies wavelet functions. And few studies have been reported on the efficient wavelet matrix transform which is very important in the wavelet HB approach. The cubic spline wavelet in [6, 7] has the explicit form and sparse transform matrix and derivative matrix. In this paper, we provide the theoretical analysis for HB method by using the cubic spline wavelet and Daubechies wavelets.The remainder of this paper is organized as follows. In Section2, we develop the HB method based on the cubic spline wavelet and Daubechies wavelets in [4] for nonlinear circuits simulations, respectively. Section 3 provides the theoretical comparison analysis in the sparsity and computation of Jacobian matrix obtained by the transform. And it is shown that the cubic spline wavelet HB method has sparser Jacobian matrix. Numerical experiments are provided in Section 4. It is concluded in Section 5.


## 2. HB Formulation Based on the Cubic Spline Wavelets and Daubechies Wavelets

### 2.1. Generalized HB Formulation

The harmonic balance (HB) method is a powerful technique for the analysis of high-frequency nonlinear circuits such as mixers, power amplifiers, and oscillators. The basic idea of HB is to expand the unknown state variablex
(
t
) in electrical circuit equations by some series x
(
t
)
=
∑
‍
X
k
v
k
(
t
). Then the problem is transformed into the frequency domain focusing on the coefficients X
k.Let us consider the general approach of HB which assumes obtaining the solutionx
(
t
) of the nonlinear modified nodal analysis (MNA) equation in [8]
(1)
C
x
˙
+
G
x
+
f
(
x
)
+
u
=
0
,
which satisfies the following periodical boundary condition:
(2)
x
(
t
+
L
)
=
x
(
t
)
,
where C and G are N
x
×
N
x matrices, x is a N
x dimensional column vector of unknown circuit variables, and u is a N
x dimensional column vector of independent sources. Let {
v
k
} be the basis; then the unknown function x
(
t
) can be expanded x
(
t
)
=
∑
‍
X
k
v
k
(
t
). To solve (1) with periodic boundary condition (2), assume that the expansion basis is periodic with period τ and [
x
l
] is a discrete vector containing values of x
(
t
) sampled in the time domain at time points [
t
l
], l
=
1
,
…
,
N
t. Then (1) can be written in the transform domain as a nonlinear algebraic equation system:
(3)
Φ
(
X
)
=
(
C
^
D
+
G
^
)
X
+
F
(
X
)
+
U
=
0
,
where
(4)
X
=
T
x
,
x
=
T
-
1
X
,
U
=
T
u
,
C
^
,
D, and G
^ are N
t
N
x
×
N
t
N
x matrices, especially, the matrix D is a representation matrix of the derivative operator d
/
d
t in expansion basis {
v
i
}
(5)
[
D
i
,
j
]
=
〈
d
d
t
v
i
,
v
j
〉
,
and, finally, T and T
-
1 are the matrices associated with the forward and inverse transform arising from the chosen expansion basis. The nonlinear matrix system (3) can be solved by Newton iterative method
(6)
J
(
X
(
i
)
)
(
X
(
i
+
1
)
-
X
(
i
)
)
=
-
Φ
(
X
(
i
)
)
,
where X
(
i
) is the solution of the ith iteration and J
(
X
) is the Jacobian matrix of Φ
(
X
)
(7)
J
(
X
)
=
[
J
k
l
(
X
)
]
=
[
∂
Φ
k
∂
X
l
]
=
∂
Φ
∂
X
=
C
^
D
+
G
^
+
T
[
∂
f
k
∂
x
l
]
T
-
1
,
l
l
l
l
l
l
l
m
k
,
l
=
1
,
…
,
(
N
t
N
x
)
.Hence, the sparsity of this Jacobian matrixJ
(
X
) affects the computational cost of iterative method. Because these matrices C
^ and G
^ have a rather sparse structure due to the MNA formulation and [
(
∂
f
k
)
/
(
∂
x
l
)
] for time-invariant systems is just a block matrix consisting of diagonal blocks, the sparsity of the Jacobian matrix J
(
X
) is determined by three matrices T, T
-
1, and the representation matrix D of the differential operator d
/
d
t.Given the base{
v
k
}
k
=
1
N, the matrices D, T, and T
-
1 are constructed before those iterative methods are used. So the sparsity of the Jacobian matrix based on these different basis functions indicates how to solve the nonlinear algebraic system. Next, we give the formulation for two kinds of wavelet bases.


### 2.2. Description of the Periodic Daubechies Wavelets

Two functionsψ and ϕ are the wavelet function and its corresponding scaling function described by Daubechies [9]. They are defined in the frame of the wavelet theory and can be constructed with finite spatial support under the following conditions:
(8)
ψ
(
t
)
=
2
∑
k
=
0
M
-
1
‍
g
k
+
1
ϕ
(
2
t
-
k
)
,
ϕ
(
t
)
=
2
∑
k
=
0
M
-
1
‍
h
k
+
1
ϕ
(
2
t
-
k
)
,
∫
-
∞
+
∞
‍
ϕ
(
t
)
d
t
=
1
,
where the coefficients {
h
k
}
k
=
0
M
-
1 and {
g
k
}
k
=
0
M
-
1 are the quadrature mirror filters (QMFs) of length L
M. The quadrature mirror filters {
h
k
} and {
g
k
} are defined
(9)
g
k
=
(
-
1
)
k
h
M
-
k
-
1
,
k
=
0,1
,
…
,
M
-
1
.The functionψ has p vanishing moments; that is,
(10)
∫
-
∞
∞
‍
ψ
(
t
)
t
m
d
t
=
0
,
0
≤
m
≤
p
-
1
.
The number M of the filter coefficients is related to the number of vanishing moments p, and M
=
2
p for the wavelets constructed in [9].We observe that once the filter{
h
k
} has been chosen, the functions ϕ and ψ can be confirmed. Moreover, due to the recursive definition of the wavelet bases, via the two-scale equation, all of the manipulations are performed with the quadrature mirror filters {
h
k
} and {
g
k
}. Especially, the wavelet transform matrix T and the derivative matrix D for the differential operator d
/
d
t can be obtained by the filters.In HB method the wavelets on the interval[
0
,
L
] are required. Hence, periodic Daubechies wavelets on the interval [
0
,
L
] are constructed by periodization. Here, we describe the discrete wavelet transform matrix by the periodic Daubechies wavelets. The discrete wavelet transform with the period L
=
2
n can be considered as a linear transformation taking the vector f
J
∈
V
J determined by its sampling data into the vector
(11)
d
=
(
c
0
,
d
0
,
d
1
,
d
2
,
d
3
,
…
,
d
J
-
1
)
T
,
where c
j stands for the scaling coefficients of the function x
(
t
) and d
j for the wavelet coefficients.This linear transform can be represented by theN
=
2
n dimensional matrix T
Daube such that
(12)
T
Daube
f
J
=
d
.
If the level of the DWT is J
≤
n, then the DWT of the sequence has exactly 2
n coefficients. The transform matrix T
Daube is composed of QMFs coefficients {
h
k
} and {
g
k
} as follows:
(13)
T
Daube
=
(
h
0
h
1
⋯
h
M
-
1
0
⋯
0
⋯
0
g
0
g
1
⋯
g
M
-
1
0
⋯
0
⋯
0
0
0
h
0
h
1
⋯
h
M
-
1
0
⋯
0
0
0
g
0
g
1
⋯
g
M
-
1
0
⋯
0
⋯
h
2
h
3
⋯
h
M
-
1
⋯
⋯
0
h
1
h
2
g
2
g
3
⋯
g
M
-
1
⋯
⋯
0
g
1
g
2
)
,
where M is the length of the filters.The periodized Daubechies wavelet HB formulation has been formulated in [4], so we have D
Daube
=
T
Daube
R
T
Daube
-
1, where R is a band limited circulant matrix with its diagonals filled by r
m in [10], where with the following properties:
(14)
r
m
≠
0
,
for
-
M
+
2
≤
m
≤
M
-
2
,
r
0
=
0
,
r
-
m
=
-
r
m
,
∑
m
‍
m
r
m
=
-
1
,
r
m
=
2
[
r
2
m
+
1
2
∑
k
=
1
M
/
2
‍
a
2
k
-
1
(
r
2
m
-
2
k
+
1
+
r
2
m
+
2
k
-
1
)
]
,
in which a
i are autocorrelation coefficients of the QMFs
(15)
a
i
=
2
∑
m
=
0
M
-
i
-
1
‍
h
~
m
h
m
+
1
,
i
=
1
,
…
,
M
-
1
.
And the matrix T
Daube
-
1 is the inverse matrix of the forward transform matrix T
Daube which satisfies T
Daube
-
1
=
T
Daube
T due to the orthogonality of the matrix T
Daube.


### 2.3. The Cubic Spline Wavelet Basis

Consider the cubic spline wavelets as the expansion base in HB technique. The cubic spline wavelets are constructed in [7], which are semiorthogonal wavelets. The high approximation rate and the interpolation property can be inherited from spline functions. Therefore, the cubic spline wavelet transform matrix T
cubic and the differential operator representation matrix D
cubic have the following properties which are suitable for HB method.Due to the periodic conditionx
(
t
+
L
)
=
x
(
t
), we must use the periodization functions of the cubic spline wavelets on the interval [
0
,
L
]
,
L
>
4. For convenience, we still denote by v
i
(
t
) the periodic function. Let us assume that expansion bases are
(16)
{
v
i
}
i
=
1
N
s
=
{
i
2
ϕ
0
,
-
1
,
ϕ
0
,
k
(
0
≤
k
≤
L
-
4
)
,
ϕ
0
,
L
-
3
,
m
ψ
j
,
k
(
0
≤
j
≤
J
-
1
,
-
1
≤
k
≤
n
j
-
2
)
}
,
where n
j
=
2
j
L
,
N
s
=
2
J
L
-
1. Correspondingly, the unknown state variable x
(
t
) is approximated by the bases of these spaces
(17)
V
J
=
V
J
-
1
⊕
W
J
-
1
⋮
=
V
0
⊕
W
0
⊕
⋯
⊕
W
J
-
1
,
where
(18)
V
0
=
span
⁡
{
ϕ
-
1
,
-
1
(
t
)
,
…
,
ϕ
-
1
,
L
-
4
(
t
)
,
ϕ
-
1
,
L
-
3
(
L
-
t
)
}
,
W
i
=
span
⁡
{
ψ
i
,
-
1
(
t
)
,
ψ
i
,
0
(
t
)
,
…
,
ψ
i
,
n
i
-
2
(
t
)
}
,
ψ
i
,
-
1
(
t
)
,
ψ
i
,
0
(
t
)
,
…
iiiiiiiiiiii
,
0
≤
i
≤
J
-
1
.
Based on the interpolation property of the cubic spline wavelets, we have
(19)
P
V
J
x
(
t
)
=
I
V
b
x
(
t
)
+
∑
j
=
0
J
-
1
‍
I
W
j
x
(
t
)
=
x
^
-
1
,
-
3
η
1
(
t
)
+
x
^
-
1
,
-
2
η
2
(
t
)
+
x
^
-
1
,
-
1
ϕ
b
(
t
)
+
∑
k
=
0
L
-
4
‍
x
^
-
1
,
k
ϕ
k
(
t
)
+
x
^
-
1
,
L
-
3
ϕ
b
(
L
-
t
)
+
x
^
-
1
,
L
-
2
η
2
(
L
-
t
)
+
x
^
-
1
,
L
-
1
η
1
(
L
-
t
)
+
∑
j
=
0
J
-
1
‍
[
∑
k
=
-
1
n
j
-
2
‍
x
^
j
,
k
ψ
j
,
k
(
t
)
]
.
Denote the expansion coefficients by a N
s
×
1 dimensional vector x
^
J,
(20)
x
^
J
=
(
x
^
-
1
,
-
3
,
…
,
x
^
-
1
,
L
-
1
,
x
^
0
,
-
1
,
…
,
x
^
0
,
n
0
-
2
,
…
,
x
^
J
-
1
,
-
1
,
m
…
,
x
^
J
-
1
,
k
,
…
,
x
^
J
-
1
,
n
J
-
2
)
T
,
that will be determined by satisfying the collocation conditions, N
s
=
2
J
L
+
3. Interpolate P
V
J at the collocation points
(21)
{
t
1
(
-
1
)
=
0
,
t
2
(
-
1
)
=
1
2
,
m
t
k
(
-
1
)
=
k
-
2
,
k
=
3
,
…
,
L
+
1
;
m
t
L
+
2
(
-
1
)
=
L
-
1
2
,
t
L
+
3
(
-
1
)
=
L
}
,
{
t
-
1
(
j
)
=
1
2
j
+
2
,
t
k
(
j
)
=
k
+
1.5
2
j
,
0
≤
k
≤
n
j
-
3
,
m
t
n
j
-
2
(
j
)
=
L
-
1
2
j
+
2
}
,
as follows:
(22)
P
V
J
x
(
t
k
-
1
)
=
x
(
t
k
-
1
)
,
1
≤
k
≤
L
+
3
,
P
V
J
x
(
t
k
j
)
=
x
(
t
k
j
)
,
j
≥
0
,
-
1
≤
k
≤
n
j
-
2
,
0
≤
j
≤
J
-
1
.
Substituting the expressions into (1), we obtain nonlinear discrete algebraic systems.Denote byT
cubic the cubic spline wavelet transform matrix. We introduce an inverse wavelet transform (IWT) T
cubic
-
1 which maps its wavelet coefficients x
^
J to discrete sample values f
J with length N
s; that is T
cubic
-
1
x
^
J
=
f
J. The inverse transform matrix T
cubic
-
1 is
(23)
T
cubic
-
1
=
(
B
M
0
M
1
⋮
M
J
-
1
)
,
where B denotes a tridiagonal matrix with dimension L
+
2 and M
j is a tridiagonal matrix with dimension 2
j
L.We obtain the derivative matrixD
cubic in [11] as follows:
(24)
D
cubic
=
H
1
-
1
H
2
,
where
(25)
H
1
=
[
λ
1
1
λ
1
2
μ
1
λ
2
2
μ
2
·
·
·
λ
i
2
μ
i
·
·
·
λ
N
s
-
1
2
μ
N
s
-
1
1
μ
N
s
-
1
]
(
N
s
+
1
)
×
(
N
s
+
1
)
,
H
2
=
[
a
1
a
2
a
3
c
1
d
1
e
1
·
·
·
·
·
·
c
i
d
i
e
i
·
·
·
c
N
s
-
1
d
N
s
-
1
e
N
s
-
1
b
3
b
2
b
1
]
(
N
s
+
1
)
×
(
N
s
+
1
)
,
and these constants in these matrices H
1 and H
2 can be referenced from the formulae (2.20a)–(2.20d) in [11].For the whole nonlinear equation system(26)
(
C
^
H
1
-
1
H
2
+
G
^
+
T
cubic
[
∂
f
k
∂
x
l
]
T
cubic
-
1
)
(
X
(
i
+
1
)
-
X
(
i
)
)
=
-
Φ
(
X
(
i
)
)
,
where H
1, H
2, and T
cubic
-
1 are tridiagonal matrices, the triangular decomposition of the tridiagonal matrix can be used to decompose the Jacobian iterative matrix.


## 2.1. Generalized HB Formulation

The harmonic balance (HB) method is a powerful technique for the analysis of high-frequency nonlinear circuits such as mixers, power amplifiers, and oscillators. The basic idea of HB is to expand the unknown state variablex
(
t
) in electrical circuit equations by some series x
(
t
)
=
∑
‍
X
k
v
k
(
t
). Then the problem is transformed into the frequency domain focusing on the coefficients X
k.Let us consider the general approach of HB which assumes obtaining the solutionx
(
t
) of the nonlinear modified nodal analysis (MNA) equation in [8]
(1)
C
x
˙
+
G
x
+
f
(
x
)
+
u
=
0
,
which satisfies the following periodical boundary condition:
(2)
x
(
t
+
L
)
=
x
(
t
)
,
where C and G are N
x
×
N
x matrices, x is a N
x dimensional column vector of unknown circuit variables, and u is a N
x dimensional column vector of independent sources. Let {
v
k
} be the basis; then the unknown function x
(
t
) can be expanded x
(
t
)
=
∑
‍
X
k
v
k
(
t
). To solve (1) with periodic boundary condition (2), assume that the expansion basis is periodic with period τ and [
x
l
] is a discrete vector containing values of x
(
t
) sampled in the time domain at time points [
t
l
], l
=
1
,
…
,
N
t. Then (1) can be written in the transform domain as a nonlinear algebraic equation system:
(3)
Φ
(
X
)
=
(
C
^
D
+
G
^
)
X
+
F
(
X
)
+
U
=
0
,
where
(4)
X
=
T
x
,
x
=
T
-
1
X
,
U
=
T
u
,
C
^
,
D, and G
^ are N
t
N
x
×
N
t
N
x matrices, especially, the matrix D is a representation matrix of the derivative operator d
/
d
t in expansion basis {
v
i
}
(5)
[
D
i
,
j
]
=
〈
d
d
t
v
i
,
v
j
〉
,
and, finally, T and T
-
1 are the matrices associated with the forward and inverse transform arising from the chosen expansion basis. The nonlinear matrix system (3) can be solved by Newton iterative method
(6)
J
(
X
(
i
)
)
(
X
(
i
+
1
)
-
X
(
i
)
)
=
-
Φ
(
X
(
i
)
)
,
where X
(
i
) is the solution of the ith iteration and J
(
X
) is the Jacobian matrix of Φ
(
X
)
(7)
J
(
X
)
=
[
J
k
l
(
X
)
]
=
[
∂
Φ
k
∂
X
l
]
=
∂
Φ
∂
X
=
C
^
D
+
G
^
+
T
[
∂
f
k
∂
x
l
]
T
-
1
,
l
l
l
l
l
l
l
m
k
,
l
=
1
,
…
,
(
N
t
N
x
)
.Hence, the sparsity of this Jacobian matrixJ
(
X
) affects the computational cost of iterative method. Because these matrices C
^ and G
^ have a rather sparse structure due to the MNA formulation and [
(
∂
f
k
)
/
(
∂
x
l
)
] for time-invariant systems is just a block matrix consisting of diagonal blocks, the sparsity of the Jacobian matrix J
(
X
) is determined by three matrices T, T
-
1, and the representation matrix D of the differential operator d
/
d
t.Given the base{
v
k
}
k
=
1
N, the matrices D, T, and T
-
1 are constructed before those iterative methods are used. So the sparsity of the Jacobian matrix based on these different basis functions indicates how to solve the nonlinear algebraic system. Next, we give the formulation for two kinds of wavelet bases.


## 2.2. Description of the Periodic Daubechies Wavelets

Two functionsψ and ϕ are the wavelet function and its corresponding scaling function described by Daubechies [9]. They are defined in the frame of the wavelet theory and can be constructed with finite spatial support under the following conditions:
(8)
ψ
(
t
)
=
2
∑
k
=
0
M
-
1
‍
g
k
+
1
ϕ
(
2
t
-
k
)
,
ϕ
(
t
)
=
2
∑
k
=
0
M
-
1
‍
h
k
+
1
ϕ
(
2
t
-
k
)
,
∫
-
∞
+
∞
‍
ϕ
(
t
)
d
t
=
1
,
where the coefficients {
h
k
}
k
=
0
M
-
1 and {
g
k
}
k
=
0
M
-
1 are the quadrature mirror filters (QMFs) of length L
M. The quadrature mirror filters {
h
k
} and {
g
k
} are defined
(9)
g
k
=
(
-
1
)
k
h
M
-
k
-
1
,
k
=
0,1
,
…
,
M
-
1
.The functionψ has p vanishing moments; that is,
(10)
∫
-
∞
∞
‍
ψ
(
t
)
t
m
d
t
=
0
,
0
≤
m
≤
p
-
1
.
The number M of the filter coefficients is related to the number of vanishing moments p, and M
=
2
p for the wavelets constructed in [9].We observe that once the filter{
h
k
} has been chosen, the functions ϕ and ψ can be confirmed. Moreover, due to the recursive definition of the wavelet bases, via the two-scale equation, all of the manipulations are performed with the quadrature mirror filters {
h
k
} and {
g
k
}. Especially, the wavelet transform matrix T and the derivative matrix D for the differential operator d
/
d
t can be obtained by the filters.In HB method the wavelets on the interval[
0
,
L
] are required. Hence, periodic Daubechies wavelets on the interval [
0
,
L
] are constructed by periodization. Here, we describe the discrete wavelet transform matrix by the periodic Daubechies wavelets. The discrete wavelet transform with the period L
=
2
n can be considered as a linear transformation taking the vector f
J
∈
V
J determined by its sampling data into the vector
(11)
d
=
(
c
0
,
d
0
,
d
1
,
d
2
,
d
3
,
…
,
d
J
-
1
)
T
,
where c
j stands for the scaling coefficients of the function x
(
t
) and d
j for the wavelet coefficients.This linear transform can be represented by theN
=
2
n dimensional matrix T
Daube such that
(12)
T
Daube
f
J
=
d
.
If the level of the DWT is J
≤
n, then the DWT of the sequence has exactly 2
n coefficients. The transform matrix T
Daube is composed of QMFs coefficients {
h
k
} and {
g
k
} as follows:
(13)
T
Daube
=
(
h
0
h
1
⋯
h
M
-
1
0
⋯
0
⋯
0
g
0
g
1
⋯
g
M
-
1
0
⋯
0
⋯
0
0
0
h
0
h
1
⋯
h
M
-
1
0
⋯
0
0
0
g
0
g
1
⋯
g
M
-
1
0
⋯
0
⋯
h
2
h
3
⋯
h
M
-
1
⋯
⋯
0
h
1
h
2
g
2
g
3
⋯
g
M
-
1
⋯
⋯
0
g
1
g
2
)
,
where M is the length of the filters.The periodized Daubechies wavelet HB formulation has been formulated in [4], so we have D
Daube
=
T
Daube
R
T
Daube
-
1, where R is a band limited circulant matrix with its diagonals filled by r
m in [10], where with the following properties:
(14)
r
m
≠
0
,
for
-
M
+
2
≤
m
≤
M
-
2
,
r
0
=
0
,
r
-
m
=
-
r
m
,
∑
m
‍
m
r
m
=
-
1
,
r
m
=
2
[
r
2
m
+
1
2
∑
k
=
1
M
/
2
‍
a
2
k
-
1
(
r
2
m
-
2
k
+
1
+
r
2
m
+
2
k
-
1
)
]
,
in which a
i are autocorrelation coefficients of the QMFs
(15)
a
i
=
2
∑
m
=
0
M
-
i
-
1
‍
h
~
m
h
m
+
1
,
i
=
1
,
…
,
M
-
1
.
And the matrix T
Daube
-
1 is the inverse matrix of the forward transform matrix T
Daube which satisfies T
Daube
-
1
=
T
Daube
T due to the orthogonality of the matrix T
Daube.


## 2.3. The Cubic Spline Wavelet Basis

Consider the cubic spline wavelets as the expansion base in HB technique. The cubic spline wavelets are constructed in [7], which are semiorthogonal wavelets. The high approximation rate and the interpolation property can be inherited from spline functions. Therefore, the cubic spline wavelet transform matrix T
cubic and the differential operator representation matrix D
cubic have the following properties which are suitable for HB method.Due to the periodic conditionx
(
t
+
L
)
=
x
(
t
), we must use the periodization functions of the cubic spline wavelets on the interval [
0
,
L
]
,
L
>
4. For convenience, we still denote by v
i
(
t
) the periodic function. Let us assume that expansion bases are
(16)
{
v
i
}
i
=
1
N
s
=
{
i
2
ϕ
0
,
-
1
,
ϕ
0
,
k
(
0
≤
k
≤
L
-
4
)
,
ϕ
0
,
L
-
3
,
m
ψ
j
,
k
(
0
≤
j
≤
J
-
1
,
-
1
≤
k
≤
n
j
-
2
)
}
,
where n
j
=
2
j
L
,
N
s
=
2
J
L
-
1. Correspondingly, the unknown state variable x
(
t
) is approximated by the bases of these spaces
(17)
V
J
=
V
J
-
1
⊕
W
J
-
1
⋮
=
V
0
⊕
W
0
⊕
⋯
⊕
W
J
-
1
,
where
(18)
V
0
=
span
⁡
{
ϕ
-
1
,
-
1
(
t
)
,
…
,
ϕ
-
1
,
L
-
4
(
t
)
,
ϕ
-
1
,
L
-
3
(
L
-
t
)
}
,
W
i
=
span
⁡
{
ψ
i
,
-
1
(
t
)
,
ψ
i
,
0
(
t
)
,
…
,
ψ
i
,
n
i
-
2
(
t
)
}
,
ψ
i
,
-
1
(
t
)
,
ψ
i
,
0
(
t
)
,
…
iiiiiiiiiiii
,
0
≤
i
≤
J
-
1
.
Based on the interpolation property of the cubic spline wavelets, we have
(19)
P
V
J
x
(
t
)
=
I
V
b
x
(
t
)
+
∑
j
=
0
J
-
1
‍
I
W
j
x
(
t
)
=
x
^
-
1
,
-
3
η
1
(
t
)
+
x
^
-
1
,
-
2
η
2
(
t
)
+
x
^
-
1
,
-
1
ϕ
b
(
t
)
+
∑
k
=
0
L
-
4
‍
x
^
-
1
,
k
ϕ
k
(
t
)
+
x
^
-
1
,
L
-
3
ϕ
b
(
L
-
t
)
+
x
^
-
1
,
L
-
2
η
2
(
L
-
t
)
+
x
^
-
1
,
L
-
1
η
1
(
L
-
t
)
+
∑
j
=
0
J
-
1
‍
[
∑
k
=
-
1
n
j
-
2
‍
x
^
j
,
k
ψ
j
,
k
(
t
)
]
.
Denote the expansion coefficients by a N
s
×
1 dimensional vector x
^
J,
(20)
x
^
J
=
(
x
^
-
1
,
-
3
,
…
,
x
^
-
1
,
L
-
1
,
x
^
0
,
-
1
,
…
,
x
^
0
,
n
0
-
2
,
…
,
x
^
J
-
1
,
-
1
,
m
…
,
x
^
J
-
1
,
k
,
…
,
x
^
J
-
1
,
n
J
-
2
)
T
,
that will be determined by satisfying the collocation conditions, N
s
=
2
J
L
+
3. Interpolate P
V
J at the collocation points
(21)
{
t
1
(
-
1
)
=
0
,
t
2
(
-
1
)
=
1
2
,
m
t
k
(
-
1
)
=
k
-
2
,
k
=
3
,
…
,
L
+
1
;
m
t
L
+
2
(
-
1
)
=
L
-
1
2
,
t
L
+
3
(
-
1
)
=
L
}
,
{
t
-
1
(
j
)
=
1
2
j
+
2
,
t
k
(
j
)
=
k
+
1.5
2
j
,
0
≤
k
≤
n
j
-
3
,
m
t
n
j
-
2
(
j
)
=
L
-
1
2
j
+
2
}
,
as follows:
(22)
P
V
J
x
(
t
k
-
1
)
=
x
(
t
k
-
1
)
,
1
≤
k
≤
L
+
3
,
P
V
J
x
(
t
k
j
)
=
x
(
t
k
j
)
,
j
≥
0
,
-
1
≤
k
≤
n
j
-
2
,
0
≤
j
≤
J
-
1
.
Substituting the expressions into (1), we obtain nonlinear discrete algebraic systems.Denote byT
cubic the cubic spline wavelet transform matrix. We introduce an inverse wavelet transform (IWT) T
cubic
-
1 which maps its wavelet coefficients x
^
J to discrete sample values f
J with length N
s; that is T
cubic
-
1
x
^
J
=
f
J. The inverse transform matrix T
cubic
-
1 is
(23)
T
cubic
-
1
=
(
B
M
0
M
1
⋮
M
J
-
1
)
,
where B denotes a tridiagonal matrix with dimension L
+
2 and M
j is a tridiagonal matrix with dimension 2
j
L.We obtain the derivative matrixD
cubic in [11] as follows:
(24)
D
cubic
=
H
1
-
1
H
2
,
where
(25)
H
1
=
[
λ
1
1
λ
1
2
μ
1
λ
2
2
μ
2
·
·
·
λ
i
2
μ
i
·
·
·
λ
N
s
-
1
2
μ
N
s
-
1
1
μ
N
s
-
1
]
(
N
s
+
1
)
×
(
N
s
+
1
)
,
H
2
=
[
a
1
a
2
a
3
c
1
d
1
e
1
·
·
·
·
·
·
c
i
d
i
e
i
·
·
·
c
N
s
-
1
d
N
s
-
1
e
N
s
-
1
b
3
b
2
b
1
]
(
N
s
+
1
)
×
(
N
s
+
1
)
,
and these constants in these matrices H
1 and H
2 can be referenced from the formulae (2.20a)–(2.20d) in [11].For the whole nonlinear equation system(26)
(
C
^
H
1
-
1
H
2
+
G
^
+
T
cubic
[
∂
f
k
∂
x
l
]
T
cubic
-
1
)
(
X
(
i
+
1
)
-
X
(
i
)
)
=
-
Φ
(
X
(
i
)
)
,
where H
1, H
2, and T
cubic
-
1 are tridiagonal matrices, the triangular decomposition of the tridiagonal matrix can be used to decompose the Jacobian iterative matrix.


## 3. Comparison Analysis

Using HB method to solve nonlinear ODEs, the Newton iterative form is obtained. Here we want to analyze the sparsity of derivative matrixD and wavelet transform matrix T of the Jacobian matrix based on two wavelets.


### 3.1. The Comparison of the Wavelet Transform MatrixesT
D
a
u
b
e and T
c
u
b
i
c

Now we analyze the sparsity of the matrixesT
Daube and T
cubic. The computation cost T
[
∂
f
k
/
∂
x
l
]
T
-
1 of Jacobian matrix J
(
X
) results from the number of nonzero elements of the wavelet transform matrix T. We analyze the nonzero elements (NZ) of matrices T
Daube and T
cubic. Let the maximum level of wavelet decomposition be J, N
=
2
n
,
J
≤
n. From [12], we have the nonzero numbers N
Z
T
Daube of the matrix T
Daube are
(27)
N
Z
T
Daube
≤
J
2
J
-
1
M
+
(
2
J
-
1
)
M
+
2
J
-
1
~
O
(
N
log
⁡
(
N
)
)
.
For the cubic spline interpolation wavelets, the transform matrix T
cubic
-
1 has the following property:
(28)
N
Z
T
cubic
-
1
=
[
3
(
L
-
1
)
-
2
]
+
∑
j
=
0
J
-
1
‍
(
2
j
×
3
×
L
-
2
)
=
3
L
×
2
J
-
2
J
-
5
~
O
(
N
s
)
.


### 3.2. Comparison of the Derivative MatricesD
D
a
u
b
e and D
c
u
b
i
c

The sparsity of the derivative matrixD is an important property of Jacobian matrix J
(
X
). According to the approach in [4], matrix D
Daube is composed of T
Daube
R
T
Daube
-
1, where T
Daube
-
1 is the transpose of the matrix T
Daube. Since both T
Daube and T
Daube
-
1 in this case are band-limited matrices, as well as R, the resulting matrix D
Daube is also a band-limited matrix. Especially, the nonzero element numbers of matrix R are (
M
-
5
)
×
N
+
2
∑
i
=
1
M
-
2
‍
i
~
O
(
N
), so we have
(29)
N
Z
D
Daube
~
O
(
N
log
⁡
(
N
)
)
.By the formulation in Section2.3, the number of nonzero element of matrix D
cubic
=
H
1
-
1
H
2 is O
(
N
s
). It follows that cubic spline wavelets yield a sparser derivative matrix than that of Daubechies wavelets. Thus, the Jacobian matrix of the cubic spline wavelets is much sparser than the periodic Daubechies wavelet.


## 3.1. The Comparison of the Wavelet Transform MatrixesT
D
a
u
b
e and T
c
u
b
i
c

Now we analyze the sparsity of the matrixesT
Daube and T
cubic. The computation cost T
[
∂
f
k
/
∂
x
l
]
T
-
1 of Jacobian matrix J
(
X
) results from the number of nonzero elements of the wavelet transform matrix T. We analyze the nonzero elements (NZ) of matrices T
Daube and T
cubic. Let the maximum level of wavelet decomposition be J, N
=
2
n
,
J
≤
n. From [12], we have the nonzero numbers N
Z
T
Daube of the matrix T
Daube are
(27)
N
Z
T
Daube
≤
J
2
J
-
1
M
+
(
2
J
-
1
)
M
+
2
J
-
1
~
O
(
N
log
⁡
(
N
)
)
.
For the cubic spline interpolation wavelets, the transform matrix T
cubic
-
1 has the following property:
(28)
N
Z
T
cubic
-
1
=
[
3
(
L
-
1
)
-
2
]
+
∑
j
=
0
J
-
1
‍
(
2
j
×
3
×
L
-
2
)
=
3
L
×
2
J
-
2
J
-
5
~
O
(
N
s
)
.


## 3.2. Comparison of the Derivative MatricesD
D
a
u
b
e and D
c
u
b
i
c

The sparsity of the derivative matrixD is an important property of Jacobian matrix J
(
X
). According to the approach in [4], matrix D
Daube is composed of T
Daube
R
T
Daube
-
1, where T
Daube
-
1 is the transpose of the matrix T
Daube. Since both T
Daube and T
Daube
-
1 in this case are band-limited matrices, as well as R, the resulting matrix D
Daube is also a band-limited matrix. Especially, the nonzero element numbers of matrix R are (
M
-
5
)
×
N
+
2
∑
i
=
1
M
-
2
‍
i
~
O
(
N
), so we have
(29)
N
Z
D
Daube
~
O
(
N
log
⁡
(
N
)
)
.By the formulation in Section2.3, the number of nonzero element of matrix D
cubic
=
H
1
-
1
H
2 is O
(
N
s
). It follows that cubic spline wavelets yield a sparser derivative matrix than that of Daubechies wavelets. Thus, the Jacobian matrix of the cubic spline wavelets is much sparser than the periodic Daubechies wavelet.


## 4. Numerical Experiments

In this section, we will give the sparsity figures of the transform matrix and the derivative matrix based on two kinds of wavelets. For simplicity, assume the matricesC
^, G
^, and [
∂
f
/
∂
x
l
] are diagonal. Figure 1 is the sparsity of T
Daube using the periodized D4 Daubechies wavelets.(a): Sparsity pattern of the periodized transform matrixT
Daube by the periodized D4 wavelets; (b): sparsity pattern of the inverse transform matrix T
cubic
-
1 of the cubic spline wavelet.
(a)
(b)In Figure2 we plot the sparsity of the derivative matrix of periodic Daubechies wavelets and the matrix H
2 or H
1 of the derivative matrix D
cubic.(a): Sparsity pattern of the derivative matrixD
Daube by the periodized D4 wavelets; (b): sparsity pattern of H
1 or H
2 of the derivative matrix of the cubic spline wavelet, where D
cubic
=
H
1
-
1
H
2.
(a)
(b)


## 5. Conclusions

In this paper, we formulate the comparison analysis of harmonic balance method based on the cubic spline wavelets and periodic Daubechies wavelets. It is shown that the cubic spline wavelet HB method has the special structure for Jacobian matrix compared to the Daubechies wavelet HB method to solve steady-state analysis of nonlinear circuits.



---

*Source: 634974-2014-04-22.xml*