# Recent Developments on Hybrid Time-Frequency Numerical Simulation Techniques for RF and Microwave Applications

Authors: Jorge F. Oliveira; José C. Pedro
Journal: Journal of Function Spaces and Applications
 (2013)

---

## Abstract

This paper reviews some of the promising doors that functional analysis techniques have recently opened in the field of electronic circuit simulation. Because of the modulated nature of radio frequency (RF) signals, the corresponding electronic circuits seem to operate in a slow time scale for the aperiodic information and another, much faster, time scale for the periodic carrier. This apparent multirate behavior can be appropriately described using partial differential equations (PDEs) within a bivariate framework, which can be solved in an efficient way using hybrid time-frequency techniques. With these techniques, the aperiodic information dimension is treated in the discrete time domain, while the periodic carrier dimension is processed in the frequency domain, in which the solution is evaluated within a space of harmonically related sinusoidal functions. The objective of this paper is thus to provide a general overview on the most important hybrid time-frequency techniques, as the ones found in commercial tools or the ones recently published in the literature.


---

## Body

## 1. Introduction

Numerical simulation plays an important role in electronics, helping engineers to verify correctness and debug circuits during their design, and so avoiding breadboarding and physical prototyping. The advantages of numerical simulation are especially significant in integrated circuits design, where manufacturing is expensive and probing internal nodes is difficult or prohibitive.Circuit simulation has emerged in the early 1970’s, and many numerical techniques have been developed and improved along the years. Radio frequency (RF) and microwave system design is a field that was an important driver for numerical simulation development, and continues to be so nowadays. Indeed, computing the solution of some current electronic circuits, as is the case of modern wireless communication systems, is still today a hot topic. In effect, serious difficulties arise when these nonlinear systems are highly heterogeneous circuits operating in multiple time scales. Current examples of these are wireless RF integrated circuits (RFICs), or systems-on-a-chip (SoC), combining RF, baseband analog, and digital blocks in the same the circuit.Signals handed by wireless communication systems can usually be described by a high frequency RF carrier modulated by some kind of slowly varying baseband information signal. Hence, the analysis of any statistically relevant information time frame requires the processing of thousands or millions of time points of the composite modulated signal, turning any conventional numerical integration of the circuit’s system of ordinary differential equations (ODEs) highly inefficient, or even impractical. However, if the waveforms produced by the circuit are not excessively demanding on the number of harmonics for a convenient frequency-domain representation, this class of problems can be efficiently simulated with hybrid time-frequency techniques. Handling the response to the slowly varying baseband information signal in the conventional time step by time-step basis, but representing the reaction to the periodic RF carrier as a small set of Fourier components (a harmonic balance algorithm for computing the steady-state response to the carrier) new circuit simulators are taking an enormous profit from functional analysis techniques. But, beyond overcoming the signals’ time-scale disparity, one of the recently proposed hybrid time-frequency techniques is also able to deal with highly heterogeneous RF circuits in an efficient way, by applying different numerical strategies to state variables in different parts (blocks) of the circuits.


## 2. Theoretical Background Material

### 2.1. Mathematical Model of an Electronic Circuit

The behavior of an electronic circuit can be described with a system of equations involving voltages, currents, charges, and fluxes. This system of equations can be constructed from a circuit description using, for example, nodal analysis, which involves applying the Kirchhoff current law to each node in the circuit, and applying the constitutive or branch equations to each circuit element. Systems generated this way have, in general, the following form(1)p(y(t))+dq(y(t))dt=x(t),
where x(t)∈Rn and y(t)∈Rn stand for the excitation (independent voltage or current sources) and state variable (node voltages and branch currents) vectors, respectively. p:ℝn→ℝn stands for all linear or nonlinear elements, as resistors, nonlinear voltage-controlled current sources, and so forth, while q:ℝn→ℝn models dynamic linear or nonlinear elements, as capacitors (represented as linear or nonlinear voltage-dependent electric charges), or inductors (represented as linear or nonlinear current-dependent magnetic fluxes).The system of (1) is, in general, a differential algebraic equations’ (DAE) system, which represents the general mathematical formulation of lumped problems. However, as reviewed in [1], this DAE circuit model formulation could even be extended to include linear distributed elements. For that, these are substituted, one-by-one, by their lumped-element equivalent circuit  models or are replaced, as whole sub-circuits, by reduced order models derived from their frequency-domain characteristics whenever larger distributed linear networks are dealt with.The substitution of distributed devices by lumped-equivalent models is especially reasonable when the size of the circuit elements is small in comparison to the wavelengths, as is the case of most emerging RF technologies (e.g., new systems on chip (SoCs), or systems in package (SiPs), integrating digital high-speed CMOS baseband processing and RFCMOS hardware).


### 2.2. Steady-State Simulation

The most natural way of simulating an electronic circuit is to numerically time-step integrate, in time domain, the ordinary differential system describing its operation. This straightforward technique was used in the first digital computer programs of circuit analysis and is still widely used nowadays. It is the core of all SPICE (which means simulation program with integrated circuit emphasis) [2] or SPICE-like computer programs. The dilemma is that these tools focus on transient analysis, and sometimes electronics designers, as is the case of RF and microwave designers,  are not interested in the circuits’ transient response, but, instead, in their steady-state regimes. This is because certain aspects of circuits’ performance are better characterized, or simply only defined, in steady-state (e.g., distortion, noise, power, gain, impedance, etc.). Time-step integration engines, as linear multistep methods, or Runge-Kutta methods, which were tailored for finding the circuit’s transient response, are not adequate for computing the steady-state because they have to pass through the lengthy process of integrating all transients and expecting them to vanish. In circuits presenting extremely different time constants, or high Q resonances, as is typically the case of RF and microwave circuits, time-step integration can be very inefficient. Indeed, in such cases, frequencies in steady-state response are much higher than the rate at which the circuit approaches steady-state or the ratio between the highest and the lowest frequency is very large. Thus, the number of discretization time steps used by the numerical integration scheme will be enormous because the time interval over which the differential equations must be numerically integrated is set by the lowest frequency or by how long the circuit takes to achieve steady-state, while the size of the time steps is constrained by the highest frequency component.It must be noted that there are several different kinds of steady-state behavior that may be of interest. The first one is DC steady-state. Here, the solution does not vary with time. Stable linear circuits driven by sinusoidal sources may exhibit a sinusoidal steady-state regime, which is characterized as being purely sinusoidal except, possibly, for some DC offset. If the steady-state response of a circuit consists of generic waveforms presenting a common period, then the circuit is said to be in a periodic steady-state. Directly computing the periodic steady-state response of an electronic circuit, without having to first integrate its transient response, involves finding the initial condition,y(t0), for the differential system that describe the circuit’s operation, such that the solution at the end of one period matches the initial condition, that is, y(t0)=y(t0+T), where T is the period. Problems of this form, those of finding the solution to a system of ordinary differential equations that satisfies constraints at two or more distinct points in time, are referred to as boundary value problems. In this particular case, we have a periodic boundary value problem that can be formulated as
(2)p(y(t))+dq(y(t))dt=x(t),y(t0)=y(t0+T),t0≤t≤t0+T,y(t)∈ℝn,
where the condition y(t0)=y(t0+T) is known as the periodic boundary condition.In the following, we will focus our attention to the most widely used technique for computing the periodic steady-state solution of RF and microwave electronic circuits: the harmonic balance method [3–5].


### 2.3. Harmonic Balance

Harmonic balance (HB) is a mature computer steady-state simulation tool that operates in the frequency domain [3]. Frequency-domain methods differ from time-domain steady-state techniques in the way that, instead of representing waveforms as a collection of time samples, they represent them using coefficients of sinusoids in trigonometric series. The main advantage of the trigonometric series approach is that the steady-state solution can often be represented accurately with a small number of terms. For example, if the circuit is linear and its inputs are all sinusoidal of the same frequency, only two terms (magnitude and phase) of the trigonometric series will represent the solution exactly, whereas an approximate time-domain solution would require a much larger number of sample points.Another advantage of operating directly in the frequency-domain is that linear dynamic operations, like differentiation or integration, are converted into simple algebraic operations, such as multiplying or dividing by frequency, respectively. For example, when analyzing linear time-invariant circuit devices, the coefficients of the response are easily evaluated by exploiting superposition within phasor analysis [6]. Computing the response of nonlinear devices is obviously more difficult than for linear devices, in part because superposition no longer applies, and also because, in general, the coefficients of the response cannot be computed directly from the coefficients of the stimulus. Nevertheless, in the case of moderate nonlinearities, the steady-state solution is typically achieved much more easily in frequency-domain than in time-domain simulators.HB handles the circuit, its excitation and its state variables in the frequency domain, which is the format normally adopted by RF designers. Because of that, it also benefits from allowing the direct inclusion of distributed devices (like dispersive transmission lines) or other circuit elements described by frequency-domain measurement data, for which we cannot find an exact time-domain representation.In order to provide a brief and illustrative explanation of the conventional HB theory, let us start by considering again the boundary value problem of (2), describing the periodic steady-state regime of an electronic circuit. For simplicity, let us momentarily suppose that we are dealing with a scalar problem, that is, that we have a simple circuit described with a unique state variable y(t), and that this circuit is driven by a single source x(t), verifying the periodic condition x(t)=x(t+T). Since the steady-state response of the circuit will be also periodic with period T, both the excitation and the steady-state solution can be expressed as the Fourier series
(3)x(t)=∑k=-∞+∞Xkejkω0t,y(t)=∑k=-∞+∞Ykejkω0t,
where ω0=2π/T is the fundamental frequency. By substituting (3) into (2), and adopting a convenient harmonic truncation at some order k=K, we will obtain
(4)p(∑k=-K+KYkejkω0t)+ddt[q(∑k=-K+KYkejkω0t)]=∑k=-K+KXkejkω0t.The HB method consists in converting this differential system into the frequency domain, in way to obtain an algebraic system of2K+1 equations, in which the unknowns are the Fourier coefficients Yk. It must be noted that since p and q are, in general, nonlinear functions, it is not possible to directly compute the Fourier coefficients Yk in this system. In fact, we only know a priori the trivial solution y(t)=0 for x(t)=0. So, we can possibly guess an initial estimate to y(t) and then adopt an iterative procedure to compute the steady-state response of the circuit. For that, we use a first-order Taylor-series expansion, in which each initial expansion point corresponds to the previous iterated solution. Indeed, we expand the left hand side of the DAE system in (2) to obtain
(5)p(y[r](t))+dq(y[r](t))dt+dp(y)dy|y=y[r]︸g(y[r])[y[r+1](t)-y[r](t)]+ddtdq(y)dy|y=y[r]︸c(y[r])[y[r+1](t)-y[r](t)]=x(t),
which results in
(6)p(∑k=-KKYk[r]ejkω0t)+ddt[q(∑k=-KKYk[r]ejkω0t)]+g(∑k=-KKYk[r]ejkω0t)[∑k=-KK(Yk[r+1]-Yk[r])ejkω0t]+ddt[c(∑k=-KKYk[r]ejkω0t)[∑k=-KK(Yk[r+1]-Yk[r])ejkω0t]]=∑k=-KKXkejkω0t.The difficulty now arising in solving (6) is that we want to transform this system entirely into the frequency domain, but we do not know how to compute the Fourier coefficients of p(·), q(·), g(·), and c(·) at each iteration r. So, one possible way to do that consists of computing each of these nonlinear functions in the time domain and then calculate their Fourier coefficients. Therefore, according to the properties of the Fourier transform, the time-domain products g(y[r])·[y[r+1](t)-y[r](t)] and c(y[r])·[y[r+1](t)-y[r](t)] will become spectral convolutions, which can be represented as matrix-vector products using the conversion matrix formulation [5, 7]. This way, and because of the orthogonality of the Fourier series, (6) can be expressed in the form
(7)P[r]+jΩQ[r]+G[r][Y[r+1]-Y[r]]+jΩC[r][Y[r+1]-Y[r]]=X,
where
(8)Y=[Y-K⋮Y0⋮YK],X=[X-K⋮X0⋮XK],jΩ=diag(-jKω0,…,0,…,jKω0).In (7), P and Q are vectors containing the Fourier coefficients of p(y(t)) and q(y(t)), respectively, and G and C denote the (2K+1)×(2K+1) conversion matrices (Toeplitz) [5, 7] corresponding to g(y(t)) and c(y(t)). If we rewrite (7) as
(9)P[r]+jΩQ[r]-X︸F(Y[r])+[G[r]+jΩC[r]︸J(Y[r])][Y[r+1]-Y[r]]=0,
we can obtain
(10)F(Y[r])+dF(Y)dY|Y=Y[r][Y[r+1]-Y[r]]=0,
in which
(11)F(Y)=P(Y)+jΩQ(Y)-X=0
is known as the harmonic balance equation, and the (2K+1)×(2K+1) composite conversion matrix
(12)J(Y)=dF(Y)dY=G(Y)+jΩC(Y)
is known as the Jacobian matrix of the error function F(Y).The iterative procedure of (5)–(12) is the so-called harmonic-Newton algorithm. In order to achieve the final solution of the problem, we have to do the following operations at each iteration r: (i) perform inverse Fourier transformation to obtain y[r](t) from Y[r]; (ii) evaluate p(y[r](t)), q(y[r](t)), g(y[r](t)), and c(y[r](t)) in time domain; (iii) calculate their Fourier coefficients to obtain P(Y[r]), Q(Y[r]), G(Y[r]), and C(Y[r]), and thus F(Y[r]) and J(Y[r]); (iv) solve the linear system of (2K+1) algebraic equations of (10) to compute the next estimate Y[r+1]. Consecutive iterations will be conducted until a final solution Y[f] satisfies the HB equation of (11) with a desired accuracy, that is, until
(13)∥F(Y[f])∥=∥P(Y[f])+jΩQ(Y[f])-X∥<tol,
where tol is an allowed error ceiling and ∥F(·)∥ stands for some norm of the error function F(·).Since in a digital computer, both time and frequency domains are represented by discrete quantities, the mathematical tools used to perform Fourier and inverse Fourier transformations are, respectively, the discrete Fourier transform (DFT) and the inverse discrete Fourier transform (IDFT) or their fast algorithms, the fast Fourier transform (FFT) and the inverse fast Fourier transform (IFFT).The system of (10) is typically a sparse linear system in the case of a generic circuit with n state variables. In general, several methods can be used to solve this system, such as direct solvers, sparse solvers, or iterative solvers. However, for very large systems, iterative solvers are usually preferred. Krylov subspace techniques [8] are a class of iterative methods for solving sparse linear systems of equations. An advantage of Krylov techniques is that (10) does not need to be fully solved in each iteration. The iterative process needs only to proceed until Y[r+1]-Y[r] is such that Y[r+1] decreases the error function. This approach to the solution, called inexact Newton, can provide significantly improved efficiency. Today, there is a general consensus that a technique called the generalized minimum residual (GMRES) [9] is the preferred one among the many available Krylov subspace techniques, for harmonic-balance analysis [10–12].The generalization of the above described harmonic-Newton algorithm to the case of a generic electronic circuit withn state variables is obviously straightforward. Indeed, in such case we will simply have
(14)Y=[Y1T,Y2T,…,YnT]T,
where each one of the Yv, v=1,…,n,  is a (2K+1)×1 vector containing the Fourier coefficients of the corresponding state variable yv(t). The jΩ matrix will be defined as
(15)jΩ=diag(-jKω0,…,jKω0︸v=1,-jKω0,…,jKω0︸v=2,…,-jKω0,…,jKω0︸v=n),
and the Jacobian matrix J(Y)=dF(Y)/dY will have a block structure, consisting of an n×n matrix of square submatrices (blocks), each of one with dimension (2K+1). Each block contains information on the sensitivity of changes in a component of the error function F(Y),  resulting from changes in a component of Y. The general block of row m and column l can be expressed as
(16)dFm(Y)dYl=dPm(Y)dYl+jΩdQm(Y)dYl,
where dPm(Y)/dYl and dQm(Y)/dYl denote, respectively, the Toeplitz conversion matrices [7] of the vectors containing the Fourier coefficients of dpm(y(t))/dyl(t) and dqm(y(t))/dyl(t).


## 2.1. Mathematical Model of an Electronic Circuit

The behavior of an electronic circuit can be described with a system of equations involving voltages, currents, charges, and fluxes. This system of equations can be constructed from a circuit description using, for example, nodal analysis, which involves applying the Kirchhoff current law to each node in the circuit, and applying the constitutive or branch equations to each circuit element. Systems generated this way have, in general, the following form(1)p(y(t))+dq(y(t))dt=x(t),
where x(t)∈Rn and y(t)∈Rn stand for the excitation (independent voltage or current sources) and state variable (node voltages and branch currents) vectors, respectively. p:ℝn→ℝn stands for all linear or nonlinear elements, as resistors, nonlinear voltage-controlled current sources, and so forth, while q:ℝn→ℝn models dynamic linear or nonlinear elements, as capacitors (represented as linear or nonlinear voltage-dependent electric charges), or inductors (represented as linear or nonlinear current-dependent magnetic fluxes).The system of (1) is, in general, a differential algebraic equations’ (DAE) system, which represents the general mathematical formulation of lumped problems. However, as reviewed in [1], this DAE circuit model formulation could even be extended to include linear distributed elements. For that, these are substituted, one-by-one, by their lumped-element equivalent circuit  models or are replaced, as whole sub-circuits, by reduced order models derived from their frequency-domain characteristics whenever larger distributed linear networks are dealt with.The substitution of distributed devices by lumped-equivalent models is especially reasonable when the size of the circuit elements is small in comparison to the wavelengths, as is the case of most emerging RF technologies (e.g., new systems on chip (SoCs), or systems in package (SiPs), integrating digital high-speed CMOS baseband processing and RFCMOS hardware).


## 2.2. Steady-State Simulation

The most natural way of simulating an electronic circuit is to numerically time-step integrate, in time domain, the ordinary differential system describing its operation. This straightforward technique was used in the first digital computer programs of circuit analysis and is still widely used nowadays. It is the core of all SPICE (which means simulation program with integrated circuit emphasis) [2] or SPICE-like computer programs. The dilemma is that these tools focus on transient analysis, and sometimes electronics designers, as is the case of RF and microwave designers,  are not interested in the circuits’ transient response, but, instead, in their steady-state regimes. This is because certain aspects of circuits’ performance are better characterized, or simply only defined, in steady-state (e.g., distortion, noise, power, gain, impedance, etc.). Time-step integration engines, as linear multistep methods, or Runge-Kutta methods, which were tailored for finding the circuit’s transient response, are not adequate for computing the steady-state because they have to pass through the lengthy process of integrating all transients and expecting them to vanish. In circuits presenting extremely different time constants, or high Q resonances, as is typically the case of RF and microwave circuits, time-step integration can be very inefficient. Indeed, in such cases, frequencies in steady-state response are much higher than the rate at which the circuit approaches steady-state or the ratio between the highest and the lowest frequency is very large. Thus, the number of discretization time steps used by the numerical integration scheme will be enormous because the time interval over which the differential equations must be numerically integrated is set by the lowest frequency or by how long the circuit takes to achieve steady-state, while the size of the time steps is constrained by the highest frequency component.It must be noted that there are several different kinds of steady-state behavior that may be of interest. The first one is DC steady-state. Here, the solution does not vary with time. Stable linear circuits driven by sinusoidal sources may exhibit a sinusoidal steady-state regime, which is characterized as being purely sinusoidal except, possibly, for some DC offset. If the steady-state response of a circuit consists of generic waveforms presenting a common period, then the circuit is said to be in a periodic steady-state. Directly computing the periodic steady-state response of an electronic circuit, without having to first integrate its transient response, involves finding the initial condition,y(t0), for the differential system that describe the circuit’s operation, such that the solution at the end of one period matches the initial condition, that is, y(t0)=y(t0+T), where T is the period. Problems of this form, those of finding the solution to a system of ordinary differential equations that satisfies constraints at two or more distinct points in time, are referred to as boundary value problems. In this particular case, we have a periodic boundary value problem that can be formulated as
(2)p(y(t))+dq(y(t))dt=x(t),y(t0)=y(t0+T),t0≤t≤t0+T,y(t)∈ℝn,
where the condition y(t0)=y(t0+T) is known as the periodic boundary condition.In the following, we will focus our attention to the most widely used technique for computing the periodic steady-state solution of RF and microwave electronic circuits: the harmonic balance method [3–5].


## 2.3. Harmonic Balance

Harmonic balance (HB) is a mature computer steady-state simulation tool that operates in the frequency domain [3]. Frequency-domain methods differ from time-domain steady-state techniques in the way that, instead of representing waveforms as a collection of time samples, they represent them using coefficients of sinusoids in trigonometric series. The main advantage of the trigonometric series approach is that the steady-state solution can often be represented accurately with a small number of terms. For example, if the circuit is linear and its inputs are all sinusoidal of the same frequency, only two terms (magnitude and phase) of the trigonometric series will represent the solution exactly, whereas an approximate time-domain solution would require a much larger number of sample points.Another advantage of operating directly in the frequency-domain is that linear dynamic operations, like differentiation or integration, are converted into simple algebraic operations, such as multiplying or dividing by frequency, respectively. For example, when analyzing linear time-invariant circuit devices, the coefficients of the response are easily evaluated by exploiting superposition within phasor analysis [6]. Computing the response of nonlinear devices is obviously more difficult than for linear devices, in part because superposition no longer applies, and also because, in general, the coefficients of the response cannot be computed directly from the coefficients of the stimulus. Nevertheless, in the case of moderate nonlinearities, the steady-state solution is typically achieved much more easily in frequency-domain than in time-domain simulators.HB handles the circuit, its excitation and its state variables in the frequency domain, which is the format normally adopted by RF designers. Because of that, it also benefits from allowing the direct inclusion of distributed devices (like dispersive transmission lines) or other circuit elements described by frequency-domain measurement data, for which we cannot find an exact time-domain representation.In order to provide a brief and illustrative explanation of the conventional HB theory, let us start by considering again the boundary value problem of (2), describing the periodic steady-state regime of an electronic circuit. For simplicity, let us momentarily suppose that we are dealing with a scalar problem, that is, that we have a simple circuit described with a unique state variable y(t), and that this circuit is driven by a single source x(t), verifying the periodic condition x(t)=x(t+T). Since the steady-state response of the circuit will be also periodic with period T, both the excitation and the steady-state solution can be expressed as the Fourier series
(3)x(t)=∑k=-∞+∞Xkejkω0t,y(t)=∑k=-∞+∞Ykejkω0t,
where ω0=2π/T is the fundamental frequency. By substituting (3) into (2), and adopting a convenient harmonic truncation at some order k=K, we will obtain
(4)p(∑k=-K+KYkejkω0t)+ddt[q(∑k=-K+KYkejkω0t)]=∑k=-K+KXkejkω0t.The HB method consists in converting this differential system into the frequency domain, in way to obtain an algebraic system of2K+1 equations, in which the unknowns are the Fourier coefficients Yk. It must be noted that since p and q are, in general, nonlinear functions, it is not possible to directly compute the Fourier coefficients Yk in this system. In fact, we only know a priori the trivial solution y(t)=0 for x(t)=0. So, we can possibly guess an initial estimate to y(t) and then adopt an iterative procedure to compute the steady-state response of the circuit. For that, we use a first-order Taylor-series expansion, in which each initial expansion point corresponds to the previous iterated solution. Indeed, we expand the left hand side of the DAE system in (2) to obtain
(5)p(y[r](t))+dq(y[r](t))dt+dp(y)dy|y=y[r]︸g(y[r])[y[r+1](t)-y[r](t)]+ddtdq(y)dy|y=y[r]︸c(y[r])[y[r+1](t)-y[r](t)]=x(t),
which results in
(6)p(∑k=-KKYk[r]ejkω0t)+ddt[q(∑k=-KKYk[r]ejkω0t)]+g(∑k=-KKYk[r]ejkω0t)[∑k=-KK(Yk[r+1]-Yk[r])ejkω0t]+ddt[c(∑k=-KKYk[r]ejkω0t)[∑k=-KK(Yk[r+1]-Yk[r])ejkω0t]]=∑k=-KKXkejkω0t.The difficulty now arising in solving (6) is that we want to transform this system entirely into the frequency domain, but we do not know how to compute the Fourier coefficients of p(·), q(·), g(·), and c(·) at each iteration r. So, one possible way to do that consists of computing each of these nonlinear functions in the time domain and then calculate their Fourier coefficients. Therefore, according to the properties of the Fourier transform, the time-domain products g(y[r])·[y[r+1](t)-y[r](t)] and c(y[r])·[y[r+1](t)-y[r](t)] will become spectral convolutions, which can be represented as matrix-vector products using the conversion matrix formulation [5, 7]. This way, and because of the orthogonality of the Fourier series, (6) can be expressed in the form
(7)P[r]+jΩQ[r]+G[r][Y[r+1]-Y[r]]+jΩC[r][Y[r+1]-Y[r]]=X,
where
(8)Y=[Y-K⋮Y0⋮YK],X=[X-K⋮X0⋮XK],jΩ=diag(-jKω0,…,0,…,jKω0).In (7), P and Q are vectors containing the Fourier coefficients of p(y(t)) and q(y(t)), respectively, and G and C denote the (2K+1)×(2K+1) conversion matrices (Toeplitz) [5, 7] corresponding to g(y(t)) and c(y(t)). If we rewrite (7) as
(9)P[r]+jΩQ[r]-X︸F(Y[r])+[G[r]+jΩC[r]︸J(Y[r])][Y[r+1]-Y[r]]=0,
we can obtain
(10)F(Y[r])+dF(Y)dY|Y=Y[r][Y[r+1]-Y[r]]=0,
in which
(11)F(Y)=P(Y)+jΩQ(Y)-X=0
is known as the harmonic balance equation, and the (2K+1)×(2K+1) composite conversion matrix
(12)J(Y)=dF(Y)dY=G(Y)+jΩC(Y)
is known as the Jacobian matrix of the error function F(Y).The iterative procedure of (5)–(12) is the so-called harmonic-Newton algorithm. In order to achieve the final solution of the problem, we have to do the following operations at each iteration r: (i) perform inverse Fourier transformation to obtain y[r](t) from Y[r]; (ii) evaluate p(y[r](t)), q(y[r](t)), g(y[r](t)), and c(y[r](t)) in time domain; (iii) calculate their Fourier coefficients to obtain P(Y[r]), Q(Y[r]), G(Y[r]), and C(Y[r]), and thus F(Y[r]) and J(Y[r]); (iv) solve the linear system of (2K+1) algebraic equations of (10) to compute the next estimate Y[r+1]. Consecutive iterations will be conducted until a final solution Y[f] satisfies the HB equation of (11) with a desired accuracy, that is, until
(13)∥F(Y[f])∥=∥P(Y[f])+jΩQ(Y[f])-X∥<tol,
where tol is an allowed error ceiling and ∥F(·)∥ stands for some norm of the error function F(·).Since in a digital computer, both time and frequency domains are represented by discrete quantities, the mathematical tools used to perform Fourier and inverse Fourier transformations are, respectively, the discrete Fourier transform (DFT) and the inverse discrete Fourier transform (IDFT) or their fast algorithms, the fast Fourier transform (FFT) and the inverse fast Fourier transform (IFFT).The system of (10) is typically a sparse linear system in the case of a generic circuit with n state variables. In general, several methods can be used to solve this system, such as direct solvers, sparse solvers, or iterative solvers. However, for very large systems, iterative solvers are usually preferred. Krylov subspace techniques [8] are a class of iterative methods for solving sparse linear systems of equations. An advantage of Krylov techniques is that (10) does not need to be fully solved in each iteration. The iterative process needs only to proceed until Y[r+1]-Y[r] is such that Y[r+1] decreases the error function. This approach to the solution, called inexact Newton, can provide significantly improved efficiency. Today, there is a general consensus that a technique called the generalized minimum residual (GMRES) [9] is the preferred one among the many available Krylov subspace techniques, for harmonic-balance analysis [10–12].The generalization of the above described harmonic-Newton algorithm to the case of a generic electronic circuit withn state variables is obviously straightforward. Indeed, in such case we will simply have
(14)Y=[Y1T,Y2T,…,YnT]T,
where each one of the Yv, v=1,…,n,  is a (2K+1)×1 vector containing the Fourier coefficients of the corresponding state variable yv(t). The jΩ matrix will be defined as
(15)jΩ=diag(-jKω0,…,jKω0︸v=1,-jKω0,…,jKω0︸v=2,…,-jKω0,…,jKω0︸v=n),
and the Jacobian matrix J(Y)=dF(Y)/dY will have a block structure, consisting of an n×n matrix of square submatrices (blocks), each of one with dimension (2K+1). Each block contains information on the sensitivity of changes in a component of the error function F(Y),  resulting from changes in a component of Y. The general block of row m and column l can be expressed as
(16)dFm(Y)dYl=dPm(Y)dYl+jΩdQm(Y)dYl,
where dPm(Y)/dYl and dQm(Y)/dYl denote, respectively, the Toeplitz conversion matrices [7] of the vectors containing the Fourier coefficients of dpm(y(t))/dyl(t) and dqm(y(t))/dyl(t).


## 3. Hybrid Time-Frequency Simulation

### 3.1. Modulated Signals

Signals containing components that vary at two or more widely separated rates are usually referred to as multirate signals and have a special incidence in RF and microwave applications, such as mixers (up/down converters), modulators, demodulators, power amplifiers, and so forth. Multirate signals can appear in RF systems due to the existence of excitation regimes of widely separated time scales (e.g., baseband stimuli and high frequency local oscillators) or because the stimuli can be, themselves, multirate signals (e.g., circuits driven by modulated signals). The general form of an amplitude and phase-modulated signal can be defined as(17)x(t)=e(t)cos(ωCt+ϕ(t)),
where e(t) and ϕ(t) are, respectively, the amplitude, or envelope, and phase slowly varying baseband signals, modulating the cos(ωCt) fast-varying carrier. Circuits driven by this kind of signals, or presenting themselves state variables of this type, are common in RF and microwave applications. Since the baseband signals have a spectral content of much lower frequency than the carrier, that is, because they are typically slowly varying signals while the carrier is a fast-varying entity, simulating nonlinear circuits containing this kind of signals is often a very challenging issue. Because the aperiodic nature of the signals obviates the use of any steady-state technique, one might think that conventional time-step integration would be the natural method for simulating such circuits. However, the large time constants of the bias networks determine long transient regimes and, as a result, the obligation of simulating a large number of carrier periods. In addition, computing the RF carrier oscillations long enough to obtain information about its envelope and phase properties is, itself, a colossal task. Time-step integration is thus inadequate for simulating this kind of problems because it is computationally expensive or prohibitive.


### 3.2. Hybrid Time-Frequency ETHB Technique

Theenvelope transient harmonic balance (ETHB) [13–16] is a hybrid time-frequency technique that was conceived to overcome the inefficiency revealed by SPICE-like engines (time-step integration schemes) when simulating circuits driven by modulated signals or presenting state variables of this type. It consists in calculating the response of the circuit to the baseband and the carrier by treating the envelope and phase in the time domain and the carrier in the frequency domain. For that, it assumes that the envelope and phase baseband signals are extremely slow when compared to the carrier, so that they can be considered as practically constant during many carrier periods. Taking this into account, ETHB samples the baseband signals in an appropriately slow time rate and assumes a staircase version of both amplitude and phase, which will conduct to a new modulated version of these signals. The steady-state response of the circuit to this new modulated version is then computed at each time step with the frequency-domain HB engine.In order to provide a very brief theoretical description of the ETHB technique, let us suppose that we have a circuit driven by a single source of the form ofx(t) in (17). If we rewrite x(t) as
(18)x(t)=e(t)∑k=-11Akejk[ωCt+ϕ(t)]=∑k=-11e(t)Akejkϕ(t)ejkωCt=∑k=-11Xk(t)ejkωCt
and assume that the circuit is stable, then all its state variables can be expressed as time-varying Fourier series
(19)y(t)=∑kYk(t)ejkωCt,
where Yk(t) represents the time-varying Fourier coefficients of y(t), which are slowly varying in the baseband time scale. Now, if we take into consideration the disparity between the baseband and the carrier time scales and assume that they are also uncorrelated, which is normally the case, then we can rewrite (17) and (19) as
(20)x(tE,tC)=e(tE)cos(ωCtC+ϕ(tE)),(21)y(tE,tC)=∑kYk(tE)ejkωCtC,
where tE is the slow baseband time scale and tC is the fast carrier time scale. Then, if we discretize the slow baseband time scale using a grid of successive time instants tE,i and adopt a convenient harmonic truncation at some order k=K, we will obtain for each tE,i a periodic boundary value problem that can be solved in the frequency domain with HB. In order to compute the whole response of the circuit, a set of successive HB equations of the form
(22)F(Y(tE,i))=P(Y(tE,i))+jΩQ(Y(tE,i))-X(tE,i)=0
has to be solved, in which X(tE,i) and Y(tE,i) represent the vectors containing the time-varying Fourier coefficients of the excitation and the solution, respectively.Two different ways can be conceived to evidence the system’s dynamics to the time-varying envelope, depending on whether the circuit’s elements’ constitutive relations are described in the frequency domain or they can be formulated in the time domain.In one possibility, we rely on the frequency-domain description of each of the constitutive elements, and so of the entire system represented in (22). Assuming that the envelope time evolution is much slower than that of the carrier, we no longer consider that each harmonic component of the carrier occupies a single frequency (constant amplitude and phase carrier) but spreads through its vicinity (slowly varying amplitude and phase modulation). For example, any dynamic linear component whose frequency-domain representation is
(23)Q(Y(ω))=H(ω)Y(ω)
can be approximated by a Taylor series (or any other polynomial or rational function) in the vicinity of each of the carrier harmonics, kωC, that is, ω=kωC+ϖ, where ϖ is a slight frequency perturbation, as
(24)H(ω)-H(kωC)≃dH(ω)dω|ω=kωCϖ+12!dH2(ω)dω2|ω=kωCϖ2+13!dH3(ω)dω3|ω=kωCϖ3+⋯=H1(kωC)ϖ+H2(kωC)ϖ2+H3(kωC)ϖ3+⋯=H~k(ϖ),
which leads to
(25)Q(Y~k(ϖ=ω-kωC))=H~k(ϖ)Y~k(ϖ)≃H1(kωC)ϖY~k(ϖ)+H2(kωC)ϖ2Y~k(ϖ)+H3(kωC)ϖ3Y~k(ϖ)+⋯≃1jjϖH1(kωC)Y~k(ϖ)+1j2(jϖ)2H2(kωC)Y~k(ϖ)+1j3(jϖ)3H3(kωC)Y~k(ϖ)+⋯,
with H~ and Y~ being the low-pass equivalent of H and Y. Since Hm(kωC)/jm is a constant, and (jϖ)mY~k(ϖ) can be interpreted as the m’th order derivative of the time-domain Yk(tE) with respect to time tE, (25) can be rewriten as
(26)Q(Yk(tE))≃H1(kωC)jdYk(tE)dtE+H2(kωC)j2d2Yk(tE)dtE2+H3(kωC)j3d3Yk(tE)dtE3+⋯,
which, substituted in (22), would evidence the desired system’s dynamics to the amplitude and phase modulations. Therefore, the ETHB technique consists in the transient simulation, in an envelope time-step by time step basis, tE,i,tE,i+1,…, of the harmonic balance equation of (22).This formulation of ETHB is, nowadays, a mature technique in the RF simulation community. However, its basic assumption constitutes also its major drawback. By requiring the envelope and phase to be extremely slowly varying signals when compared to the carrier frequency, this mixed frequency-time technique becomes restricted to circuits whose stimuli occupy only a small fraction of the available bandwidth.In an alternative ETHB formulation, we assume that every element can be described in the time domain. Hence, we can substitute the time-varying Fourier description of (21) into (1) and then treat the carrier time, tC, in the frequency domain—converting the DAE system into an algebraic one—but keeping the envelope time, tE, in the time domain. This way, we obtain another hybrid time-frequency description of the system that no longer suffers from the narrow bandwidth restriction just mentioned and whose formulation and solution will be discussed in more detail in Section 3.4.


### 3.3. Multivariate Formulation

We will now introduce a powerful strategy for analyzing nonlinear circuits handling amplitude and/or phase modulated signals, as with any other kind of multirate signals. This strategy consists in using multiple time variables to describe the multirate behavior, and it is based on the fact that multirate signals can be represented much more efficiently if they are defined as functions of two or more time variables, that is, if they are defined asmultivariate functions [17, 18]. With this multivariate formulation, circuits will be no longer described by ordinary differential algebraic equations in the one-dimensional time t but, instead, by partial differential algebraic systems.Let us consider the amplitude and phase-modulated signal of (17), and let us define its bivariate form as
(27)x^(t1,t2)=e(t1)cos(ωCt2+ϕ(t1)),
where t1 is the slow envelope time scale and t2 is the fast carrier time scale. As can be seen, x^(t1,t2) is a periodic function with respect to t2 but not to t1, that is,
(28)x^(t1,t2)=x^(t1,t2+T2),T2=2πωC,
and, in general, this bivariate form requires far fewer points to represent numerically the original signal, especially when the  t1 and t2 time scales are widely separated [17, 18].Let us now consider the differential algebraic equations’ (DAEs) system of (1), describing the behavior of a generic RF circuit driven by the envelope-modulated signal of (17). Taking the above considerations into account, we will adopt the following procedure: for the slowly varying parts (envelope time scale) of the expressions of vectors x(t) and y(t), t is replaced by t1; for the fast-varying parts (RF carrier time scale), t is replaced by t2. The application of this bivariate strategy to the DAE system of (1) converts it into the following multirate partial differential algebraic equations’ (MPDAEs) system [17, 18]:
(29)p(y^(t1,t2))+∂q(y^(t1,t2))∂t1+∂q(y^(t1,t2))∂t2=x^(t1,t2).
The mathematical relation between (1) and (29) establishes that if x^(t1,t2) and y^(t1,t2) satisfy (29), then the univariate forms x(t)=x^(t,t) and y(t)=y^(t,t) satisfy (1) [18]. Therefore, the univariate solutions of (1) are available on diagonal lines t1=t, t2=t, along the bivariate solutions of (29), that is, y(t) may be retrieved from its bivariate form y^(t1,t2), by simply setting t1=t2=t. Consequently, if one wants to obtain the univariate solution in a generic [0,tFinal] interval due to the periodicity of the problem in the t2 dimension we will have
(30)y(t)=y^(t,tmodT2)
on the rectangular domain [0,tFinal]×[0,T2], where tmodT2 represents the remainder of division of t by T2. The main advantage of this MPDAE approach is that it can result in significant improvements in simulation speed when compared to DAE-based alternatives [17–20].Envelope-modulated responses to excitations of the form of (17) correspond to a combination of initial and periodic boundary conditions for the MPDAE. This means that the bivariate forms of these solutions can be obtained by numerically solving the following initial-boundary value problem [18]
(31)p(y^(t1,t2))+∂q(y^(t1,t2))∂t1+∂q(y^(t1,t2))∂t2=x^(t1,t2),y^(0,t2)=g(t2),y^(t1,0)=y^(t1,T2),
on the rectangle [0,tFinal]×[0,T2]. g(·) is a given initial-condition function defined on [0,T2], satisfying g(0)=g(T2)=y(0), and the periodic boundary condition y^(t1,0)=y^(t1,T2) is due to the periodicity of the problem in the t2 fast carrier time scale. The reason why bivariate envelope-modulated solutions do not need to be evaluated on the entire  [0,tFinal]×[0,tFinal] domain (which would be computationally very expensive and would turn the multivariate strategy useless), and are restricted to the rectangle [0,tFinal]×[0,T2], is because the solutions repeat along the t2 time axis.


### 3.4. Multitime Envelope Transient Harmonic Balance

Multitime envelope transient harmonic balance is an improved version of the previously described ETHB technique, which is based on the multivariate formulation [21, 22]. For achieving an intuitive explanation of the multitime envelope transient harmonic balance let us consider the initial-boundary value problem of (31), and let us also consider the semidiscretization of the rectangular domain [0,tFinal]×[0,T2] in the t1 slow time dimension defined by the grid
(32)0=t1,0<t1,1<⋯<t1,i-1<t1,i<⋯<t1,K1=tFinal,h1,i=t1,i-t1,i-1,
where K1 is the total number of steps in t1. If we replace the derivatives of the MPDAE in t1 with a finite-differences approximation (e.g., the Backward Euler rule), then we obtain for each slow time instant t1,i, from i=1 to i=K1, the periodic boundary value problem defined by
(33)p(y^i(t2))+q(y^i(t2))-q(y^i-1(t2))h1,i+dq(y^i(t2))dt2=x^(t1,i,t2),y^i(0)=y^i(T2),
where y^i(t2)≃y^(t1,i,t2). This means that, once y^i-1(t2) is known, the solution on the next slow time instant, y^i(t2), is obtained by solving (33). Thus, for obtaining the whole solution y^ in the entire domain [0,tFinal]×[0,T2], a total of K1 boundary value problems have to be solved. With multitime ETHB, each one of these periodic boundary value problems is solved using the harmonic balance method. The corresponding HB system for each slow time instant t1,i is the n×(2K+1) algebraic equations set given by
(34)P(Y^(t1,i))+Q(Y^(t1,i))-Q(Y^(t1,i-1))h1,i+jΩQ(Y^(t1,i))=X^(t1,i),
where X^(t1,i) and Y^(t1,i) are the vectors containing the Fourier coefficients of the excitation sources and of the solution (the state variables), respectively, at t1=t1,i. P(·) and Q(·) are unknown functions, jΩ is the diagonal matrix (15), and the Y^(t1,i) vector can be expressed as
(35)Y^(t1,i)=[Y^1(t1,i)T,Y^2(t1,i)T,…,Y^n(t1,i)T]T,
where each one of the state variable frequency components, Y^v(t1,i),  v=1,…,n,  is a (2K+1)×1 vector defined as
(36)Y^v(t1,i)=[Yv,-K(t1,i),…,Yv,0(t1,i),…,Yv,K(t1,i)]T.As seen in Section2.3, since p(·) and q(·) are in general nonlinear functions, one possible way to compute P(·) and Q(·) in (34) consists in evaluating p(·) and q(·) in the time domain and then calculate its Fourier coefficients. The HB system of (34) can be rewriten as
(37)F(Y^(t1,i))=P(Y^(t1,i))+Q(Y^(t1,i))-Q(Y^(t1,i-1))h1,i+jΩQ(Y^(t1,i))-X^(t1,i)=0,
or, in its simplified form, as
(38)F(Y^(t1,i))=0,
in which F(Y^(t1,i)) is the error function at t1=t1,i. In order to solve the nonlinear algebraic system of (38) a Newton-Raphson iterative solver is usually used. In this case, the Newton-Raphson algorithm conducts us to
(39)F(Y^[r](t1,i))+dF(Y^(t1,i))dY^(t1,i)|Y^(t1,i)=Y^[r](t1,i)[Y^[r+1](t1,i)-Y^[r](t1,i)]=0,
which means that at each iteration r, we have to solve a linear system of n×(2K+1) equations to compute the new estimate Y^[r+1](t1,i). Consecutive Newton iterations will be computed until a desired accuracy is achieved, that is, until ∥F(Y^(t1,i))∥<tol, where tol is the allowed error ceiling.The system of (39) involves the derivative of the vector F(Y^(t1,i)), with respect to the vector Y^(t1,i). The result is a matrix, the so-called Jacobian of F(Y^(t1,i)),
(40)J(Y^(t1,i))=dF(Y^(t1,i))dY^(t1,i)=[dF1(Y^(t1,i))dY^1(t1,i)dF1(Y^(t1,i))dY^2(t1,i)⋯dF1(Y^(t1,i))dY^n(t1,i)dF2(Y^(t1,i))dY^1(t1,i)dF2(Y^(t1,i))dY^2(t1,i)⋯dF2(Y^(t1,i))dY^n(t1,i)⋯⋯⋯⋯dFn(Y^(t1,i))dY^1(t1,i)dFn(Y^(t1,i))dY^2(t1,i)⋯dFn(Y^(t1,i))dY^n(t1,i)].
In the same way as in Section 2.3, this matrix has a block structure, consisting of an n×n matrix of square submatrices (blocks), each one with dimension (2K+1). The general block of row m and column l can now be expressed as
(41)dFm(Y^(t1,i))dY^l(t1,i)=dPm(Y^(t1,i))dY^l(t1,i)+1h1.idQm(Y^(t1,i))dY^l(t1,i)+jΩdQm(Y^(t1,i))dY^l(t1,i).In summary, multitime ETHB handles the solution dependence ont2 in frequency domain, while treating the course of the solution to t1 in time domain. So, it is a hybrid time-frequency technique which is similar to the ETHB engine previously reported in Section 3.2. However, an important advantage of multitime ETHB over conventional ETHB is that it does not suffer from bandwidth limitations [21]. For example, in circuits driven by envelope modulated signals, the only restriction that has to be imposed is that the modulating signal and the carrier must not be correlated in time (which is typically the case).


## 3.1. Modulated Signals

Signals containing components that vary at two or more widely separated rates are usually referred to as multirate signals and have a special incidence in RF and microwave applications, such as mixers (up/down converters), modulators, demodulators, power amplifiers, and so forth. Multirate signals can appear in RF systems due to the existence of excitation regimes of widely separated time scales (e.g., baseband stimuli and high frequency local oscillators) or because the stimuli can be, themselves, multirate signals (e.g., circuits driven by modulated signals). The general form of an amplitude and phase-modulated signal can be defined as(17)x(t)=e(t)cos(ωCt+ϕ(t)),
where e(t) and ϕ(t) are, respectively, the amplitude, or envelope, and phase slowly varying baseband signals, modulating the cos(ωCt) fast-varying carrier. Circuits driven by this kind of signals, or presenting themselves state variables of this type, are common in RF and microwave applications. Since the baseband signals have a spectral content of much lower frequency than the carrier, that is, because they are typically slowly varying signals while the carrier is a fast-varying entity, simulating nonlinear circuits containing this kind of signals is often a very challenging issue. Because the aperiodic nature of the signals obviates the use of any steady-state technique, one might think that conventional time-step integration would be the natural method for simulating such circuits. However, the large time constants of the bias networks determine long transient regimes and, as a result, the obligation of simulating a large number of carrier periods. In addition, computing the RF carrier oscillations long enough to obtain information about its envelope and phase properties is, itself, a colossal task. Time-step integration is thus inadequate for simulating this kind of problems because it is computationally expensive or prohibitive.


## 3.2. Hybrid Time-Frequency ETHB Technique

Theenvelope transient harmonic balance (ETHB) [13–16] is a hybrid time-frequency technique that was conceived to overcome the inefficiency revealed by SPICE-like engines (time-step integration schemes) when simulating circuits driven by modulated signals or presenting state variables of this type. It consists in calculating the response of the circuit to the baseband and the carrier by treating the envelope and phase in the time domain and the carrier in the frequency domain. For that, it assumes that the envelope and phase baseband signals are extremely slow when compared to the carrier, so that they can be considered as practically constant during many carrier periods. Taking this into account, ETHB samples the baseband signals in an appropriately slow time rate and assumes a staircase version of both amplitude and phase, which will conduct to a new modulated version of these signals. The steady-state response of the circuit to this new modulated version is then computed at each time step with the frequency-domain HB engine.In order to provide a very brief theoretical description of the ETHB technique, let us suppose that we have a circuit driven by a single source of the form ofx(t) in (17). If we rewrite x(t) as
(18)x(t)=e(t)∑k=-11Akejk[ωCt+ϕ(t)]=∑k=-11e(t)Akejkϕ(t)ejkωCt=∑k=-11Xk(t)ejkωCt
and assume that the circuit is stable, then all its state variables can be expressed as time-varying Fourier series
(19)y(t)=∑kYk(t)ejkωCt,
where Yk(t) represents the time-varying Fourier coefficients of y(t), which are slowly varying in the baseband time scale. Now, if we take into consideration the disparity between the baseband and the carrier time scales and assume that they are also uncorrelated, which is normally the case, then we can rewrite (17) and (19) as
(20)x(tE,tC)=e(tE)cos(ωCtC+ϕ(tE)),(21)y(tE,tC)=∑kYk(tE)ejkωCtC,
where tE is the slow baseband time scale and tC is the fast carrier time scale. Then, if we discretize the slow baseband time scale using a grid of successive time instants tE,i and adopt a convenient harmonic truncation at some order k=K, we will obtain for each tE,i a periodic boundary value problem that can be solved in the frequency domain with HB. In order to compute the whole response of the circuit, a set of successive HB equations of the form
(22)F(Y(tE,i))=P(Y(tE,i))+jΩQ(Y(tE,i))-X(tE,i)=0
has to be solved, in which X(tE,i) and Y(tE,i) represent the vectors containing the time-varying Fourier coefficients of the excitation and the solution, respectively.Two different ways can be conceived to evidence the system’s dynamics to the time-varying envelope, depending on whether the circuit’s elements’ constitutive relations are described in the frequency domain or they can be formulated in the time domain.In one possibility, we rely on the frequency-domain description of each of the constitutive elements, and so of the entire system represented in (22). Assuming that the envelope time evolution is much slower than that of the carrier, we no longer consider that each harmonic component of the carrier occupies a single frequency (constant amplitude and phase carrier) but spreads through its vicinity (slowly varying amplitude and phase modulation). For example, any dynamic linear component whose frequency-domain representation is
(23)Q(Y(ω))=H(ω)Y(ω)
can be approximated by a Taylor series (or any other polynomial or rational function) in the vicinity of each of the carrier harmonics, kωC, that is, ω=kωC+ϖ, where ϖ is a slight frequency perturbation, as
(24)H(ω)-H(kωC)≃dH(ω)dω|ω=kωCϖ+12!dH2(ω)dω2|ω=kωCϖ2+13!dH3(ω)dω3|ω=kωCϖ3+⋯=H1(kωC)ϖ+H2(kωC)ϖ2+H3(kωC)ϖ3+⋯=H~k(ϖ),
which leads to
(25)Q(Y~k(ϖ=ω-kωC))=H~k(ϖ)Y~k(ϖ)≃H1(kωC)ϖY~k(ϖ)+H2(kωC)ϖ2Y~k(ϖ)+H3(kωC)ϖ3Y~k(ϖ)+⋯≃1jjϖH1(kωC)Y~k(ϖ)+1j2(jϖ)2H2(kωC)Y~k(ϖ)+1j3(jϖ)3H3(kωC)Y~k(ϖ)+⋯,
with H~ and Y~ being the low-pass equivalent of H and Y. Since Hm(kωC)/jm is a constant, and (jϖ)mY~k(ϖ) can be interpreted as the m’th order derivative of the time-domain Yk(tE) with respect to time tE, (25) can be rewriten as
(26)Q(Yk(tE))≃H1(kωC)jdYk(tE)dtE+H2(kωC)j2d2Yk(tE)dtE2+H3(kωC)j3d3Yk(tE)dtE3+⋯,
which, substituted in (22), would evidence the desired system’s dynamics to the amplitude and phase modulations. Therefore, the ETHB technique consists in the transient simulation, in an envelope time-step by time step basis, tE,i,tE,i+1,…, of the harmonic balance equation of (22).This formulation of ETHB is, nowadays, a mature technique in the RF simulation community. However, its basic assumption constitutes also its major drawback. By requiring the envelope and phase to be extremely slowly varying signals when compared to the carrier frequency, this mixed frequency-time technique becomes restricted to circuits whose stimuli occupy only a small fraction of the available bandwidth.In an alternative ETHB formulation, we assume that every element can be described in the time domain. Hence, we can substitute the time-varying Fourier description of (21) into (1) and then treat the carrier time, tC, in the frequency domain—converting the DAE system into an algebraic one—but keeping the envelope time, tE, in the time domain. This way, we obtain another hybrid time-frequency description of the system that no longer suffers from the narrow bandwidth restriction just mentioned and whose formulation and solution will be discussed in more detail in Section 3.4.


## 3.3. Multivariate Formulation

We will now introduce a powerful strategy for analyzing nonlinear circuits handling amplitude and/or phase modulated signals, as with any other kind of multirate signals. This strategy consists in using multiple time variables to describe the multirate behavior, and it is based on the fact that multirate signals can be represented much more efficiently if they are defined as functions of two or more time variables, that is, if they are defined asmultivariate functions [17, 18]. With this multivariate formulation, circuits will be no longer described by ordinary differential algebraic equations in the one-dimensional time t but, instead, by partial differential algebraic systems.Let us consider the amplitude and phase-modulated signal of (17), and let us define its bivariate form as
(27)x^(t1,t2)=e(t1)cos(ωCt2+ϕ(t1)),
where t1 is the slow envelope time scale and t2 is the fast carrier time scale. As can be seen, x^(t1,t2) is a periodic function with respect to t2 but not to t1, that is,
(28)x^(t1,t2)=x^(t1,t2+T2),T2=2πωC,
and, in general, this bivariate form requires far fewer points to represent numerically the original signal, especially when the  t1 and t2 time scales are widely separated [17, 18].Let us now consider the differential algebraic equations’ (DAEs) system of (1), describing the behavior of a generic RF circuit driven by the envelope-modulated signal of (17). Taking the above considerations into account, we will adopt the following procedure: for the slowly varying parts (envelope time scale) of the expressions of vectors x(t) and y(t), t is replaced by t1; for the fast-varying parts (RF carrier time scale), t is replaced by t2. The application of this bivariate strategy to the DAE system of (1) converts it into the following multirate partial differential algebraic equations’ (MPDAEs) system [17, 18]:
(29)p(y^(t1,t2))+∂q(y^(t1,t2))∂t1+∂q(y^(t1,t2))∂t2=x^(t1,t2).
The mathematical relation between (1) and (29) establishes that if x^(t1,t2) and y^(t1,t2) satisfy (29), then the univariate forms x(t)=x^(t,t) and y(t)=y^(t,t) satisfy (1) [18]. Therefore, the univariate solutions of (1) are available on diagonal lines t1=t, t2=t, along the bivariate solutions of (29), that is, y(t) may be retrieved from its bivariate form y^(t1,t2), by simply setting t1=t2=t. Consequently, if one wants to obtain the univariate solution in a generic [0,tFinal] interval due to the periodicity of the problem in the t2 dimension we will have
(30)y(t)=y^(t,tmodT2)
on the rectangular domain [0,tFinal]×[0,T2], where tmodT2 represents the remainder of division of t by T2. The main advantage of this MPDAE approach is that it can result in significant improvements in simulation speed when compared to DAE-based alternatives [17–20].Envelope-modulated responses to excitations of the form of (17) correspond to a combination of initial and periodic boundary conditions for the MPDAE. This means that the bivariate forms of these solutions can be obtained by numerically solving the following initial-boundary value problem [18]
(31)p(y^(t1,t2))+∂q(y^(t1,t2))∂t1+∂q(y^(t1,t2))∂t2=x^(t1,t2),y^(0,t2)=g(t2),y^(t1,0)=y^(t1,T2),
on the rectangle [0,tFinal]×[0,T2]. g(·) is a given initial-condition function defined on [0,T2], satisfying g(0)=g(T2)=y(0), and the periodic boundary condition y^(t1,0)=y^(t1,T2) is due to the periodicity of the problem in the t2 fast carrier time scale. The reason why bivariate envelope-modulated solutions do not need to be evaluated on the entire  [0,tFinal]×[0,tFinal] domain (which would be computationally very expensive and would turn the multivariate strategy useless), and are restricted to the rectangle [0,tFinal]×[0,T2], is because the solutions repeat along the t2 time axis.


## 3.4. Multitime Envelope Transient Harmonic Balance

Multitime envelope transient harmonic balance is an improved version of the previously described ETHB technique, which is based on the multivariate formulation [21, 22]. For achieving an intuitive explanation of the multitime envelope transient harmonic balance let us consider the initial-boundary value problem of (31), and let us also consider the semidiscretization of the rectangular domain [0,tFinal]×[0,T2] in the t1 slow time dimension defined by the grid
(32)0=t1,0<t1,1<⋯<t1,i-1<t1,i<⋯<t1,K1=tFinal,h1,i=t1,i-t1,i-1,
where K1 is the total number of steps in t1. If we replace the derivatives of the MPDAE in t1 with a finite-differences approximation (e.g., the Backward Euler rule), then we obtain for each slow time instant t1,i, from i=1 to i=K1, the periodic boundary value problem defined by
(33)p(y^i(t2))+q(y^i(t2))-q(y^i-1(t2))h1,i+dq(y^i(t2))dt2=x^(t1,i,t2),y^i(0)=y^i(T2),
where y^i(t2)≃y^(t1,i,t2). This means that, once y^i-1(t2) is known, the solution on the next slow time instant, y^i(t2), is obtained by solving (33). Thus, for obtaining the whole solution y^ in the entire domain [0,tFinal]×[0,T2], a total of K1 boundary value problems have to be solved. With multitime ETHB, each one of these periodic boundary value problems is solved using the harmonic balance method. The corresponding HB system for each slow time instant t1,i is the n×(2K+1) algebraic equations set given by
(34)P(Y^(t1,i))+Q(Y^(t1,i))-Q(Y^(t1,i-1))h1,i+jΩQ(Y^(t1,i))=X^(t1,i),
where X^(t1,i) and Y^(t1,i) are the vectors containing the Fourier coefficients of the excitation sources and of the solution (the state variables), respectively, at t1=t1,i. P(·) and Q(·) are unknown functions, jΩ is the diagonal matrix (15), and the Y^(t1,i) vector can be expressed as
(35)Y^(t1,i)=[Y^1(t1,i)T,Y^2(t1,i)T,…,Y^n(t1,i)T]T,
where each one of the state variable frequency components, Y^v(t1,i),  v=1,…,n,  is a (2K+1)×1 vector defined as
(36)Y^v(t1,i)=[Yv,-K(t1,i),…,Yv,0(t1,i),…,Yv,K(t1,i)]T.As seen in Section2.3, since p(·) and q(·) are in general nonlinear functions, one possible way to compute P(·) and Q(·) in (34) consists in evaluating p(·) and q(·) in the time domain and then calculate its Fourier coefficients. The HB system of (34) can be rewriten as
(37)F(Y^(t1,i))=P(Y^(t1,i))+Q(Y^(t1,i))-Q(Y^(t1,i-1))h1,i+jΩQ(Y^(t1,i))-X^(t1,i)=0,
or, in its simplified form, as
(38)F(Y^(t1,i))=0,
in which F(Y^(t1,i)) is the error function at t1=t1,i. In order to solve the nonlinear algebraic system of (38) a Newton-Raphson iterative solver is usually used. In this case, the Newton-Raphson algorithm conducts us to
(39)F(Y^[r](t1,i))+dF(Y^(t1,i))dY^(t1,i)|Y^(t1,i)=Y^[r](t1,i)[Y^[r+1](t1,i)-Y^[r](t1,i)]=0,
which means that at each iteration r, we have to solve a linear system of n×(2K+1) equations to compute the new estimate Y^[r+1](t1,i). Consecutive Newton iterations will be computed until a desired accuracy is achieved, that is, until ∥F(Y^(t1,i))∥<tol, where tol is the allowed error ceiling.The system of (39) involves the derivative of the vector F(Y^(t1,i)), with respect to the vector Y^(t1,i). The result is a matrix, the so-called Jacobian of F(Y^(t1,i)),
(40)J(Y^(t1,i))=dF(Y^(t1,i))dY^(t1,i)=[dF1(Y^(t1,i))dY^1(t1,i)dF1(Y^(t1,i))dY^2(t1,i)⋯dF1(Y^(t1,i))dY^n(t1,i)dF2(Y^(t1,i))dY^1(t1,i)dF2(Y^(t1,i))dY^2(t1,i)⋯dF2(Y^(t1,i))dY^n(t1,i)⋯⋯⋯⋯dFn(Y^(t1,i))dY^1(t1,i)dFn(Y^(t1,i))dY^2(t1,i)⋯dFn(Y^(t1,i))dY^n(t1,i)].
In the same way as in Section 2.3, this matrix has a block structure, consisting of an n×n matrix of square submatrices (blocks), each one with dimension (2K+1). The general block of row m and column l can now be expressed as
(41)dFm(Y^(t1,i))dY^l(t1,i)=dPm(Y^(t1,i))dY^l(t1,i)+1h1.idQm(Y^(t1,i))dY^l(t1,i)+jΩdQm(Y^(t1,i))dY^l(t1,i).In summary, multitime ETHB handles the solution dependence ont2 in frequency domain, while treating the course of the solution to t1 in time domain. So, it is a hybrid time-frequency technique which is similar to the ETHB engine previously reported in Section 3.2. However, an important advantage of multitime ETHB over conventional ETHB is that it does not suffer from bandwidth limitations [21]. For example, in circuits driven by envelope modulated signals, the only restriction that has to be imposed is that the modulating signal and the carrier must not be correlated in time (which is typically the case).


## 4. Advanced Hybrid Time-Frequency Simulation

One limitation of the ETHB and multitime ETHB engines is that they do not perform any distinction between nodes or blocks within the circuit, that is to say that they treat all the circuit’s state variables in the same way. Thus, if the circuit evidences some heterogeneity, as is the case of modern wireless architectures combining radio frequency, baseband analog, and digital blocks in the same circuit, these tools cannot benefit from such feature. To overcome this difficulty an innovative mixed mode time-frequency technique was recently proposed by the authors [23, 24]. This technique splits the circuit’s state variables (node voltages and mesh currents) into fast and slowly varying subsets, treating the former with multitime ETHB and the later with a SPICE-like engine (a time-step integration scheme). This way, the strong nonlinearities of the circuit are appropriately evaluated in the time domain, while the moderate ones are computed in the frequency domain [23, 24].


### 4.1. Time-Domain Latency within the Multivariate Formulation

In order to provide an illustrative explanation of the issues under discussion in this section, let us start by considering an RF circuit in which some of its state variables (node voltages and branch currents) are fast carrier envelope-modulated waveforms, while the remaining state variables are slowly varying aperiodic signals. For concreteness, let us suppose that the signals(42)y1(t)=∑k=-KKYk(t)ejkωCt,y2(t)=e(t)
are two distinct state variables in different parts of the circuit. Yk(t) represents the Fourier coefficients of y1(t), which are slowly varying in the baseband time scale, ωC is the carrier frequency, and e(t) is a slowly varying aperiodic baseband function. We will denote signals of the form of y1(t) as active and signals of the form of y2(t) as latent. The latency revealed by y2(t) indicates that this variable belongs to a circuit block where there are no fluctuations dictated by the fast carrier. Consequently, due to its slowness, it can be represented efficiently with much less sample points than y1(t). On the other hand, since it does not evidence any periodicity, it cannot be processed with harmonic balance. On the contrary, if the number of harmonics K is not too large, the fast carrier oscillation components of y1(t) can be efficiently computed in the frequency domain. Therefore, it is straightforward to conclude that if we want to simulate circuits having such signal format disparities in an efficient way, distinct numerical strategies will be required.Let us now consider the bivariate forms ofy1(t) and y2(t) denoted by y^1(t1,t2) and y^2(t1,t2) and defined as
(43)y^1(t1,t2)=∑k=-KKYk(t1)ejkωCt2,y^2(t1,t2)=e(t1),
where t1 and t2 are, respectively, the slow envelope time dimension and the fast carrier time dimension. As we can see, y^2(t1,t2) has no dependence on t2, so it has no fluctuations in the fast time axis. In fact, it is so because y2(t) does not oscillate at the carrier frequency. Consequently, for each slow time instant t1,i defined on the grid of (32), while y^1(t1,i,t2) is a waveform that has to be represented by a certain quantity k=-K,…,K of harmonic components, y^2(t1,i,t2) is merely a constant (DC) signal that can be simply represented by the k=0 component. Therefore, there is no necessity to perform the conversion between time and frequency domains for y^2(t1,i,t2), which means that this state variable can be processed in a purely time-domain scheme.


### 4.2. Mixed Mode Time-Frequency Technique

In the above, we illustrated that bivariate forms of latent state variables have no undulations in thet2 fast time scale. So, while active state variables have to be represented by a set of (2K+1) harmonic components arranged in vectors of the form of (36), latent state variables can be represented as scalar quantities, that is,
(44)Y^v(t1,i)=Yv,0(t1,i)=y^v(t1,i).By considering this, it is straightforward to conclude that the size of theY^(t1,i) vector defined by (35) can be considerably reduced, as can be the total number of equations in the HB system of (37). An additional and crucial detail is that there is no longer obligation to perform the conversion between time and frequency domains for the latent state variables expressed in the form of (44), as well as for the components of F(Y^(t1,i)) corresponding to latent blocks of the circuit. Since the k=0 order Fourier coefficient Yv,0(t1,i) is exactly the same as the constant t2 time value y^v(t1,i), the use of the discrete Fourier transform (DFT) and the inverse discrete Fourier transform (IDFT)—or their fast algorithms, that is, the fast Fourier transform (FFT) and the inverse fast Fourier transform (IFFT)—will be required only for components in the HB system of (37) having dependence on active state variables. Significant Jacobian matrix size reductions will be achieved, too. In effect, by taking into consideration this multirate characteristic (the subset circuit latency), some of the blocks of (40) will be merely 1×1 scalar elements that contain dc information on the sensitivity of changes in components of F(Y^(t1,i)) resulting from changes in latent components of Y^(t1,i).With this strategy of partitioning the circuit into active and latent subcircuits (blocks), significant computation and memory savings can be achieved when finding the solution of (37). Indeed, with the state variable Y^(t1,i) and the error function F(Y^(t1,i)) vector size reductions, as also the resulting Jacobian J(Y^(t1,i)) matrix size reduction, it is possible to avoid dealing with large linear systems in the iterations of (39). Thus, a less computationally expensive Newton-Raphson iterative solver is required.


## 4.1. Time-Domain Latency within the Multivariate Formulation

In order to provide an illustrative explanation of the issues under discussion in this section, let us start by considering an RF circuit in which some of its state variables (node voltages and branch currents) are fast carrier envelope-modulated waveforms, while the remaining state variables are slowly varying aperiodic signals. For concreteness, let us suppose that the signals(42)y1(t)=∑k=-KKYk(t)ejkωCt,y2(t)=e(t)
are two distinct state variables in different parts of the circuit. Yk(t) represents the Fourier coefficients of y1(t), which are slowly varying in the baseband time scale, ωC is the carrier frequency, and e(t) is a slowly varying aperiodic baseband function. We will denote signals of the form of y1(t) as active and signals of the form of y2(t) as latent. The latency revealed by y2(t) indicates that this variable belongs to a circuit block where there are no fluctuations dictated by the fast carrier. Consequently, due to its slowness, it can be represented efficiently with much less sample points than y1(t). On the other hand, since it does not evidence any periodicity, it cannot be processed with harmonic balance. On the contrary, if the number of harmonics K is not too large, the fast carrier oscillation components of y1(t) can be efficiently computed in the frequency domain. Therefore, it is straightforward to conclude that if we want to simulate circuits having such signal format disparities in an efficient way, distinct numerical strategies will be required.Let us now consider the bivariate forms ofy1(t) and y2(t) denoted by y^1(t1,t2) and y^2(t1,t2) and defined as
(43)y^1(t1,t2)=∑k=-KKYk(t1)ejkωCt2,y^2(t1,t2)=e(t1),
where t1 and t2 are, respectively, the slow envelope time dimension and the fast carrier time dimension. As we can see, y^2(t1,t2) has no dependence on t2, so it has no fluctuations in the fast time axis. In fact, it is so because y2(t) does not oscillate at the carrier frequency. Consequently, for each slow time instant t1,i defined on the grid of (32), while y^1(t1,i,t2) is a waveform that has to be represented by a certain quantity k=-K,…,K of harmonic components, y^2(t1,i,t2) is merely a constant (DC) signal that can be simply represented by the k=0 component. Therefore, there is no necessity to perform the conversion between time and frequency domains for y^2(t1,i,t2), which means that this state variable can be processed in a purely time-domain scheme.


## 4.2. Mixed Mode Time-Frequency Technique

In the above, we illustrated that bivariate forms of latent state variables have no undulations in thet2 fast time scale. So, while active state variables have to be represented by a set of (2K+1) harmonic components arranged in vectors of the form of (36), latent state variables can be represented as scalar quantities, that is,
(44)Y^v(t1,i)=Yv,0(t1,i)=y^v(t1,i).By considering this, it is straightforward to conclude that the size of theY^(t1,i) vector defined by (35) can be considerably reduced, as can be the total number of equations in the HB system of (37). An additional and crucial detail is that there is no longer obligation to perform the conversion between time and frequency domains for the latent state variables expressed in the form of (44), as well as for the components of F(Y^(t1,i)) corresponding to latent blocks of the circuit. Since the k=0 order Fourier coefficient Yv,0(t1,i) is exactly the same as the constant t2 time value y^v(t1,i), the use of the discrete Fourier transform (DFT) and the inverse discrete Fourier transform (IDFT)—or their fast algorithms, that is, the fast Fourier transform (FFT) and the inverse fast Fourier transform (IFFT)—will be required only for components in the HB system of (37) having dependence on active state variables. Significant Jacobian matrix size reductions will be achieved, too. In effect, by taking into consideration this multirate characteristic (the subset circuit latency), some of the blocks of (40) will be merely 1×1 scalar elements that contain dc information on the sensitivity of changes in components of F(Y^(t1,i)) resulting from changes in latent components of Y^(t1,i).With this strategy of partitioning the circuit into active and latent subcircuits (blocks), significant computation and memory savings can be achieved when finding the solution of (37). Indeed, with the state variable Y^(t1,i) and the error function F(Y^(t1,i)) vector size reductions, as also the resulting Jacobian J(Y^(t1,i)) matrix size reduction, it is possible to avoid dealing with large linear systems in the iterations of (39). Thus, a less computationally expensive Newton-Raphson iterative solver is required.


## 5. Performance of the Methods

The performance and the efficiency of the ETHB and multitime ETHB techniques were already attested and recognized by the RF and microwave community. In the same way, the performance and the efficiency of the advanced hybrid technique described in the previous section (the mixed mode time-frequency simulation technique) were also already demonstrated through its application to several illustrative examples of practical relevance. Indeed, electronic circuits with distinct configurations and levels of complexity were especially selected to illustrate the significant gains in computational speed that can be achieved when simulating the circuits with this method [23, 24]. Nevertheless, in order to provide the reader with a realistic idea of the potential of this recently proposed technique, we included in this section a brief comparison between this method and the previous state-of-the-art multitime ETHB. For that, we considered two distinct circuits: the resistive FET mixer depicted in Figure 1 and the RF polar transmitter described in [23] and depicted in Figure 2.Figure 1
Simplified resistive FET mixer used in wireless transmitters.Figure 2
RF polar transmitter with a hybrid envelope amplifier [23].The circuits were simulated in MATLAB with the mixed mode time-frequency simulation technique versus the multitime ETHB. In our experiments a dynamic step size control tool was used in thet1 slow time scale, and we considered K=9 as the maximum harmonic order for the HB evaluations. Numerical computation times (in seconds) for simulations in the [0,0.5 μs] and [0,5.0 μs] intervals are presented in Tables 1 and 2.Table 1
Computation times—resistive FET mixer.
Simulation time interval
Mixed mode time-frequency technique
Multitime ETHB
Speedup (approx.)
[0, 0.5 μs]
2.1 s
4.6 s
2.2
[0, 5.0 μs]
19.3 s
42.5 s
2.2Table 2
Computation times—RF polar transmitter.
Simulation time interval
Mixed mode time-frequency technique
Multitime ETHB
Speedup (approx.)
[0, 0.5 μs]
3.2 s
51.5 s
16.1
[0, 5.0 μs]
25.3 s
484.6 s
19.2As we can see, speedups of approximately 2 times were obtained for the simulation of the resistive FET mixer, and speedups of more than one order of magnitude were obtained for the RF polar transmitter. These efficiency gains were achieved without compromising accuracy. Indeed, for both cases, the maximum discrepancy between solutions (for all the circuits’ state variables) was on the order of 10−8.The choice of these two circuits, which have different levels of complexity, was to illustrate how the computational efficiency is more evident as the ratio between the number of active and latent state variables is increased. In the first example, this ratio is 1, whereas in the second one this ratio is 4.5.


## 6. Conclusion

Although significant advancement has been made in RF and microwave circuit simulation along the years, the use of more elaborate functional analysis techniques has kept this subject a hot topic of scientific and practical engineering interest. Indeed, emerging wireless communication technologies continuously bring new challenges to this scientific field, as is now the case of heterogeneous RF circuits containing state variables of distinct formats and running on widely separated time scales. Taking into account the popularity of HB, but mostly ETHB, in the RF and microwave community, in this paper we have briefly reviewed the use of some functional analysis methods to address numerical simulation challenges using hybrid time-frequency techniques. A comparison between two state-of-the-art hybrid techniques in terms of computational speed is also included to evidence the efficiency gains that can be achieved by partitioning heterogeneous circuits into blocks, treating latent blocks in a one-dimensional space, and active ones in a bidimensional space.



---

*Source: 578326-2013-02-17.xml*