## 1. Introduction

In the last two decades, radio frequency (RF) and microwave system design has been found as a significant part of the electronic semiconductor industry’s portfolio. Over the years, the necessity of continuously providing new wireless systems’ functionalities and higher transmission rates, as also the need to improve transmitters’ efficiency, has been gradually reshaping wireless architectures. Heterogeneous circuits combining baseband blocks, digital blocks, and RF blocks, in the same substrate, are commonly found today. Hence, RF and microwave circuit simulation has been conducted to an increasingly challenging scenario of heterogeneous broadband and strongly nonlinear wireless communication circuits, presenting a wide variety of slowly varying and fast changing state variables (node voltages and branch currents). Thus, RF and microwave design has been an important booster for numerical simulation and device modeling development.

In general, waveforms processed by wireless communication systems can be expressed by a high-frequency RF carrier modulated by some kind of slowly varying baseband aperiodic signal (the information signal). Therefore, the evaluation of any relevant information time window requires the computation of thousands or millions of time instants of the composite modulated signal, turning any conventional numerical time-step integration of the circuits’ systems of differential algebraic equations highly inefficient. However, if the waveforms do not require too many harmonic components for a convenient frequency-domain representation, this category of circuits 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), hybrid time-frequency techniques are playing an important role in RF and microwave circuit simulation.

Beyond overcoming the signals’ time-scale disparity, the partitioned time-frequency technique discussed in Section 3.2 is also able to efficiently simulate highly heterogeneous RF networks, by splitting the circuits into different subsets (blocks) and computing their state variables with distinct numerical schemes.

## 2. Theoretical background material

### 2.1. Mathematical model of an electronic circuit

Dynamic behavior of an electronic circuit can be modeled by a system of differential algebraic equations (DAE) involving electric voltages, currents and charges, and magnetic fluxes. The DAE system can, in general, be formulated as

in which*∈*

**b**(t)*ℝ*stands for the excitation vector (independent voltage or current sources) and

^{n}*∈*

**x**(t)*ℝ*represents the state variable vector (node voltages and branch currents).

^{n}**(·) models the memoryless elements of the circuit, as is the case of linear or nonlinear resistors, linear or nonlinear controlled sources, etc.**

*f***(·) models the dynamic elements, as capacitors or inductors, represented as voltage-dependent electric charges or current-dependent magnetic fluxes, respectively.**

*q*This system of (1) can be constructed from a circuit description using, for example, nodal analysis, which involves applying Kirchhoff currents’ law to each node in the circuit, and applying the constitutive or branch equations to each circuit element. Hence, it represents the general mathematical formulation of lumped problems. However, as reviewed in [1], this DAE circuit model formulation can also include linear distributed elements. For that, the distributed devices are substituted by their lumped-element equivalent circuit models, or are replaced, as whole sub-circuits, by reduced order models derived from their frequency-domain characteristics. It must be noted that 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 integrating digital high-speed CMOS baseband processing and RFCMOS hardware in the same substrate.

### 2.2. Transient simulation

Obtaining the solution of (1) over a specified time interval **x**(t_{0}*)*=**x**_{0} is what is usually known as an *initial value problem*, and evaluating such solution is frequently referred to as *transient analysis*. The most natural way to compute * x*(

*t*) is to numerically time-step integrate (1) directly in time domain. So, it should be of no surprise that this straightforward technique was used in the first digital computer programs of circuit analysis and is still nowadays widely used. It is present in all SPICE (which means simulation program with integrated circuit emphasis) or SPICE-like computer programs [2]. In order to numerically time-step integrate the DAE system of (1) commercial tools use initial value solvers, such as linear multistep methods (LMM) [3–5], or one-step methods, i.e., Runge-Kutta (RK) methods [3–5]. Either LMM or RK families can offer a great diversity of explicit and implicit (iterative) numerical schemes, suitable to compute the numerical solution of different types of initial value problems with a desired accuracy.

### 2.3. Steady-state simulation

Although SPICE-like computer programs (which were initially conceived to compute the transient response of electronic circuits) are still widely used nowadays, RF and microwave designers’ interest normally resides on the steady-state response. The reason for that is some properties of the circuits are better described, or simply only defined, in steady-state (e.g., harmonic or intermodulation distortion, noise, power, gain, impedance, etc.). Initial value solvers, 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.

Computing the periodic steady-state response of an electronic circuit can be formulated as finding out a starting condition (left boundary), * x*(

*t*

_{0}), for the DAE system modeling the circuit that leads to a solution obeying the final condition (right boundary)

**x**(t_{0}

*)*=

**x**(t_{0}+

*T)*, with

*T*being the period. In mathematics, these problems are usually known as

*periodic boundary value problems*. Taking into account the formulation of (1), these problems will have here the following form,

**x**(t_{0}

*)*=

**x**(t_{0}+

*T)*is known as the

*periodic boundary condition*.

In order to numerically solve (2), a solution that simultaneously satisfies the differential system and the two-point periodic boundary condition has to be computed. A particular technique has been found especially useful for RF circuit simulation: the *harmonic balance* (HB) [6–8] method.

### 2.4. Harmonic balance

Harmonic balance (HB) [6–8] handles the circuit, its excitation and its state variables in the frequency-domain. Because of that, it benefits from allowing the direct inclusion of distributed devices (like dispersive transmission lines), or other circuit elements described by frequency-domain measurement data. 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 those using coefficients of sinusoids in trigonometric series. As a consequence, under moderate nonlinearities, the steady-state solution is typically achieved much more easily in the frequency domain than in the time domain.

In periodic steady-state, any stimulus *b*_{s}(*t*) or state variable *x*_{v}(*t*) can be expressed as a Fourier series

*K*is the order adopted for a convenient harmonic truncation. The HB method consists in converting the DAE system of (2) into the frequency domain, to obtain the following

*n*(2

*K*+1) algebraic equations system in which the unknowns are the Fourier coefficients of the state variables of the circuit where each one of the

*x*

_{v}(

*t*).

**and Q are vectors containing the Fourier coefficients of**

*F***(**

*f**(*

**x***t*)) and

**(**

*q**(*

**x***t*)), respectively, and

**G**and

**C**denote the

*=*

**g(x)***and*

**df(x)/dx***=*

**c(x)***d*.

**q**(**x**)/d**x***j*is the imaginary unit and Ω is a diagonal matrix defined as

The system of (4) can be rewritten as

or, in its simplified form, as in which is known as the*harmonic balance system*, and the

*Jacobian matrix*of the error function

**H**(

**X**). This system of (9) is iteratively solved according to (8), until a sufficiently accurate solution

**X**

^{[f]}is achieved, i.e., until where

*δ*is the allowed residual size, and

**H**(•).

### 2.5. Multivariate formulation

The multivariate formulation is a powerful strategy that emerged in the late 1990s, playing an important role in RF circuit simulation today. It was first introduced by Brachtendorf et al. [10] as a sophisticated derivation of quasi-periodic harmonic balance, followed by Roychowdhury [11], who demonstrated that the multivariate formulation can be an efficient strategy to analyze circuits running on distinct time scales. The multivariate formulation uses multiple time variables (artificial time scales) to describe the multirate behavior of the circuits. Thus, it is suitable to describe typical multirate regimes of operation present in RF and microwave systems, as is the case of circuits handling amplitude and/or phase-modulated signals, quasiperiodic signals, or any other kind of multirate signals containing a periodic component.

The main achievements of the multivariate formulation are due to the fact that multirate signals can be represented much more efficiently if they are defined as functions of two or more time variables (artificial time scales), i.e., if they are defined as multivariate functions [11–16]. Therefore, as we see in Section 2.5.2, circuits will be no longer described by ordinary differential algebraic equations in the one-dimensional time *t*, but, instead, by partial differential algebraic systems.

#### 2.5.1. Multivariate representations

The multivariate (multidimensional) strategy is easily illustrated by applying it to a bi-dimensional problem (two distinct time scales). So, let us consider, for example, an amplitude-modulated RF carrier of the form

where*e*(

*t*) is the envelope (a slowly varying signal), while

Consider now the following bidimensional definition for *b*(*t*),

*t*

_{E}is the slow envelope time scale and

*t*

_{C}is the fast carrier time scale. In this particular case,

*t*

_{C}but not to

*t*

_{E}, i.e.,

The univariate form, *b*(*t*), with **Figure 1** for the [0, 25 ns] time interval. The corresponding bivariate form, **Figure 2** for the rectangular region [0, 25 ns] × [0, 0.5 ns]. There, it can be appreciated that *b*(*t*), thus allowing a more compact representation with fewer samples. Furthermore, due to the periodicity of *t*_{C}, we know that its plot repeats over the rest of this time axis. Thus, the bivariate form plotted in **Figure 2** contains all the information necessary to recover the original univariate form depicted in **Figure 1**.

#### 2.5.2. Multirate partial differential algebraic equations’ systems

Let us consider a general nonlinear RF circuit described by the differential algebraic equations’ system of (1), and let us suppose that this circuit is driven by the envelope-modulated signal of (12). Considering the above stated, we are able to reformulate the excitation ** b**(

*t*) and the state variables

**(**

*x**t*) vectors as bidimensional entities, in which

*t*is replaced by

*t*

_{E}for the slowly varying parts (the envelope time scale) and by

*t*

_{C}for the fast-varying parts (the RF carrier time scale). This bidimensional formulation converts the DAE system of (1) into the following

*multirate partial differential algebraic equations’*(MPDAE) system [11]:

The mathematical relation between (1) and (15) establishes that if ** x**(

*t*) may be retrieved from its bivariate form

*t*

_{E}=

*t*

_{C}=

*t*. Consequently, if one wants to obtain the univariate solution in the generic [0,

*t*

_{Final}] interval, due to the periodicity of the problem in the

*t*

_{C}dimension we will have

*t*mod

*T*

_{C}represents the remainder of division of

*t*by

*T*

_{C}. The main advantage of this MPDAE approach is that it can result in significant improvements in simulation speed when compared to DAE-based alternatives.

#### 2.5.3. Initial and boundary conditions for envelope-modulated regimes

Dynamical behavior of RF circuits driven by stimuli of the form of (12) can be described by the MPDAE system of (15) together with a set of initial and periodic boundary conditions. In fact, bivariate forms of the circuits’ state variables can be achieved by computing the solution of the following initial periodic-boundary value problem

on the rectangle*i*(·) is a given initial-condition function defined on [0,

*T*

_{C}], satisfying

*=*

**i**(0)*, and*

**i**(_{TC})=x(0)*(t*=

_{E}0)*(t*is the periodic boundary condition due to the periodicity of the problem in the

_{E},T_{C})*t*

_{c}fast carrier time scale.

The reason why bivariate envelope-modulated solutions do not need to be evaluated on the entire *t*_{C} time axis. The way how univariate solutions are recovered from their multivariate forms was already defined above by (16).

## 3. Hybrid time-frequency techniques for computing the solution of MPDAEs

In this section, we will finally discuss the hybrid time-frequency numerical techniques that can be used to evaluate the solution of MPDAEs describing the operation of nonlinear electronic radio frequency circuits running in an aperiodic slow time scale and a periodic fast time scale. Section 3.1 addresses an efficient technique often referred to as *multitime envelope transient harmonic balance* (multitime ETHB). Then, Section 3.2 presents an advanced partitioned time-frequency technique, which is an improved version of multitime ETHB and has demonstrated to be even more efficient than this technique.

### 3.1. Multitime envelope transient harmonic balance

Let us consider the initial-boundary value problem of (17) and let us define a semi-discretization of the rectangular domain *t*_{E} slow time dimension described by the following general non uniform grid

*K*

_{E}represents the total number of steps in

*t*

_{E}and

*h*

_{E,i}denotes the grid size at each time step

*i*. If we replace the derivatives of the MPDAE in

*t*

_{E}with a finite-differences approximation (e.g., a backward differentiation formula, the modified trapezoidal rule, etc.), then we obtain for each slow time instant

*t*

_{E,i}, from

*i*= 1 to

*i*=

*K*

_{E}, a periodic boundary value problem in

*t*

_{C}. For simplicity, and clarity, let us suppose that the Backward Euler rule is used. In such a case, we obtain where

*(t*. Thus, once

_{E,i},t_{C})*K*

_{E}periodic boundary value problems if we want to obtain the solution

*t*

_{E,i}, the resultant HB system is the

*n*(2

*K*+ 1) algebraic equation system defined by where

*t*

_{E}=

*t*

_{E,i}.

*f*(·) and

*q*(·) in the time domain and then calculating their Fourier coefficients. Ω is the diagonal matrix (6), and the

*v*= 1,…,

*n*, is a (2

*K*+1)×1 vector defined as

The HB system of (20) can be rewritten as

or simply as in which*t*

_{E}=

*t*

_{E,i}. In general, the nonlinear algebraic system of (24) is iteratively solved using Newton’s method which requires that we have to solve a linear system of

*n*(2

*K*+ 1) equations at each iteration

*r*to compute the new estimate

*δ*is the allowed residual size.

The system of (25) requires the computation of the Jacobian matrix

(26) |

This matrix has a block structure, consisting of *n*×*n* square submatrices (blocks), each one with dimension (2*K* + 1). The general block of row *m* and column *l* can be expressed as

### 3.2. Partitioned time-frequency technique

Although multitime ETHB can take advantage of the signals’ time rate disparity, it does not take into account the circuit’s heterogeneities, i.e., it uses the same numerical algorithm to compute all the circuit’s state variables. Thus, if the circuit evidences some heterogeneity (e.g., modern wireless architectures combining RF, baseband analog circuitry, and digital components in the same circuit), this tool cannot benefit from such a feature. This lack of ability to perform some distinction between nodes or blocks within the circuit had already been identified by Rizzoli et al. [17] and is the main limitation of multitime ETHB. To cope with this deficiency, the partitioned time-frequency technique separates the circuit’s state variables (node voltages and branch currents) into fast (*active*) and slowly varying (*latent*) subsets. That implies the MPDAE system of (15) to be first considered as coupled active-latent MPDAE subsystems

(28) |

*x*

_{A}(

*t*) and

*x*

_{L}(

*t*) are the vectors containing, respectively, the fast-varying and the slowly varying state variables. As we will see, with this partition stratagem, fast-varying state variables can be computed with multitime ETHB, while slowly varying ones are being evaluated with a unidimensional time-step integration scheme. This tactic also allows the moderate nonlinearities to be treated in the frequency domain, while severe nonlinearities are appropriately evaluated in the time domain [16].

With the purpose of providing an elucidatory explanation of the partitioned time-frequency technique, let us consider a typical wireless system, composed of RF and baseband blocks. In such a case, the state variables in the RF block can be described as fast carrier envelope modulated waveforms defined as

while state variables in the baseband block can be seen as slowly varying aperiodic functions of the formIn (30), *X*_{k,v}(*t*) represents the Fourier coefficients of *x*_{A,v}(*t*), which are slowly varying in the baseband time scale, and *f*_{C} is the high-frequency carrier. As stated above, signals of the form *x*_{A,v}(*t*) will be denoted as active, whereas signals of the form *x*_{L,v}(*t*) will be designated as latent. The latency (slowness) of *x*_{L,v}(*t*) indicates that these variables belong to a circuit block where there are no fluctuations dictated by the fast carrier. Thus, it is straightforward to conclude that all of the *x*_{L,v}(*t*) can be efficiently represented with much less sample points than any of the *x*_{A,v}(*t*). Moreover, since the *x*_{L,v}(*t*) state variables do not evidence any periodicity, they cannot be evaluated in the frequency domain. In contrast, if the number of harmonics *K* is not too large, the fast carrier oscillation components of *x*_{A,v}(*t*) can be efficiently computed with harmonic balance. Taking the above into account we can easily conclude that distinct numerical strategies will be required if we want to simulate, in an efficient way, circuits having such signal format disparities.

In the following we provide a brief theoretical description of the partitioned time-frequency technique fundamentals. For that, let us now consider the bivariate forms of *x*_{A,v}(*t*) and *x*_{L,v}(*t*), denoted by

*t*

_{E}and

*t*

_{C}are, respectively, the slow envelope time dimension and the fast carrier time dimension. As can be seen, since the

*t*

_{c}, they have no fluctuations in the fast time axis. The reason is that they belong to a circuit block where there are no carrier frequency oscillations. As a result, for each slow time instant

*t*

_{E,i}defined on the grid of (18), each of the

*k*= 0 component. Therefore, there is no necessity to perform the conversion between time and frequency domains for

*t*

_{Ei}, each of the

*k*= −

*K*,…,

*K*, i.e., each of the

*K*+ 1 harmonic components for a convenient frequency domain representation. In summary, while active state variables have to be represented by a set of 2

*K*+ 1 Fourier coefficients arranged in (2

*K*+ 1)×1 vectors of the form latent state variables can be represented as 1×1 scalar quantities, i.e., they can be simply represented as

By considering this, we can easily deduce that the size of the *k* = 0 order Fourier coefficient *X*_{v,0}(*t*_{E,i}) is exactly the same as the constant *t*_{C} time value

Considerable Jacobian

With the state variable

## 4. Efficiency of the partitioned time-frequency technique

The effectiveness of the multitime ETHB technique is nowadays widely recognized by the RF and microwave community. The efficiency of the partitioned time-frequency simulation technique described in the previous section was also already established, as a consequence of the considerable reductions in the computational effort required to obtain the numerical solution of several RF circuits with distinct topologies and levels of complexity [16]. Even so, a brief comparison between this method, the previous state-of-the-art multitime ETHB and a conventional univariate time-step integration scheme (SPICE-like simulation), is included in this section. This will help the reader to get a perception of the potential of the partitioned hybrid technique. For that, we considered the RF mixer (frequency translation device) depicted in **Figure 3** as the illustrative application example. The circuit was simulated in MATLAB with three different techniques: (i) the partitioned time-frequency simulation technique, (ii) the multitime ETHB, and (iii) the Gear-2 multistep method [5] (a time-step integrator commonly used by SPICE-like commercial simulators).

Numerical computation times for simulations in the [0, 1.0 *μ*s] and [0, 10.0 *μ*s] intervals are presented in **Table 1**. For simplicity, in the hybrid time-frequency techniques we assumed a uniform grid in the *t*_{E} slow time scale (we have chosen *h*_{E} = 10 ns as the step size in that dimension) and we considered *K* = 11 as the maximum harmonic order for the HB evaluations in the *t*_{C} fast carrier time scale.

By comparing the CPU times obtained with the methods, we can attest the superiority of the partitioned time-frequency method. Indeed, speedups of more than two times were obtained with this method in comparison to multitime ETHB. We can also attest the inefficiency of univariate time-step integration when dealing with RF problems. Finally, it must be noted that the efficiency gain of the partitioned time-frequency technique was achieved without compromising accuracy. Indeed, the maximum discrepancy between solutions computed with this technique and multitime ETHB was in the order of 10^{−7} for all the state variables of the circuit.

## 5. Conclusions

In this chapter, we have presented a partitioned time-frequency numerical technique especially designed for the efficient simulation of RF circuits operating in a periodic fast time scale and an aperiodic slow time scale. This technique can be viewed as a wise combination of multitime ETHB based on a multivariate formulation, with a conventional univariate time-step integration scheme. With this technique fast changing (active) state variables are computed in a bivariate mixed time-frequency domain, whereas slowly varying (latent) state variables are evaluated in the natural one-dimensional time domain. By partitioning the circuits into active and latent parts and exploiting the fact there is no obligation to perform conversion between time and frequency for latent blocks, considerable reductions in the computational effort can be achieved without compromising the accuracy of the results. Although the speedups obtained with the simulation of the illustrative application example presented in Section 4 are already notable, it must be noted that higher efficiency gains should be expected when simulating RF networks containing a number of latent blocks larger than the active ones.