Hybrid Time-Frequency Numerical Simulation of Electronic Radio Frequency Systems

Abstract

This chapter is devoted to the discussion of a hybrid frequency-time CAD tool especially designed for the efficient numerical simulation of nonlinear electronic radio frequency circuits operating in an aperiodic slow time scale and a periodic fast time scale. Circuits driven by envelope-modulated signals, in which the baseband signal (the information) is aperiodic and has a spectral content of much lower frequency than the periodic carrier, are typical examples of practical interest involving such time evolution rates. The discussed method is tailored to take advantage of the circuits and signals heterogeneity and so will benefit from the time-domain latency of some state variables in the circuits. Because the aperiodic slowly varying state variables are treated only in time domain, the proposed method can be seen as a hybrid scheme combining multitime envelope transient harmonic balance based on a multivariate formulation, with a purely time-step integration scheme.

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

$\frac{dq\left(x\left(t\right)\right)}{dt}+f\left(x\left(t\right)\right)=b\left(t\right),$E1

in which b(t)∈ℝⁿ stands for the excitation vector (independent voltage or current sources) and x(t)∈ℝⁿ represents the state variable vector (node voltages and branch currents). f(·) models the memoryless elements of the circuit, as is the case of linear or nonlinear resistors, linear or nonlinear controlled sources, etc. q(·) models the dynamic elements, as capacitors or inductors, represented as voltage-dependent electric charges or current-dependent magnetic fluxes, respectively.

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 $\left[t_0,t_{Final}\right]$ with a specific initial condition x(t₀)=x₀ 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₀), for the DAE system modeling the circuit that leads to a solution obeying the final condition (right boundary) x(t₀)=x(t₀+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,

$\frac{dq\left(x\left(t\right)\right)}{dt}+f\left(x\left(t\right)\right)=b\left(t\right),x\left(t_0\right)=x\left(t_0+T\right),t_0\le t\le t_0+T,$E2

where the condition x(t₀)=x(t₀+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

$b_s\left(t\right)={\displaystyle \sum _{k=-K}^{+K}{B}_k{e}^{jk{\omega }_{0}t}},x_v\left(t\right)={\displaystyle \sum _{k=-K}^{+K}{X}_k{e}^{jk{\omega }_{0}t}},$E3

where ${\omega }_0=2\pi /T$ is the fundamental frequency and 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(2K+1) algebraic equations system

$j\Omega Q^{\left[r\right]}+F^{\left[r\right]}+j\Omega C^{\left[r\right]}\left[X^{\left[r+1\right]}-X^{\left[r\right]}\right]+G^{\left[r\right]}\left[X^{\left[r+1\right]}-X^{\left[r\right]}\right]=B,$E4

in which the unknowns are the Fourier coefficients of the state variables of the circuit

$X={\left[X_1{}^T,X_2{}^T,\dots ,X_n{}^T\right]}^T,$E5

where each one of the $X_v,v=1,\dots ,n$, is a $\left(2K+1\right)\times 1$ vector containing the Fourier coefficients of the corresponding state variable x_v(t). F and Q are vectors containing the Fourier coefficients of f(x(t)) and q(x(t)), respectively, and G and C denote the $n\left(2K+1\right)\times n\left(2K+1\right)$ conversion matrices (Toeplitz) [7,9] corresponding to g(x)=df(x)/dx and c(x)=dq(x)/dx. j is the imaginary unit and Ω is a diagonal matrix defined as

$\Omega =\text{diag}\left(\underset{v=1}{\underbrace{-K{\omega }_0,\dots ,K{\omega }_0}},\underset{v=2}{\underbrace{-K{\omega }_0,\dots ,K{\omega }_0}},\dots ,\underset{v=n}{\underbrace{-K{\omega }_0,\dots ,K{\omega }_0}}\right).$E6

The system of (4) can be rewritten as

$\underset{H\left(X^{\left[r\right]}\right)}{\underbrace{j\Omega Q^{\left[r\right]}+F^{\left[r\right]}-B}}+\left[\underset{J\left(X^{\left[r\right]}\right)}{\underbrace{j\Omega C^{\left[r\right]}+G^{\left[r\right]}}}\right]\left[X^{\left[r+1\right]}-X^{\left[r\right]}\right]=0,$E7

or, in its simplified form, as

$H\left(X^{\left[r\right]}\right)+{\frac{dH\left(X\right)}{dX}|}_{X=X^{\left[r\right]}}\left[X^{\left[r+1\right]}-X^{\left[r\right]}\right]=0,$E8

in which

$H\left(X\right)=j\Omega Q\left(X\right)+F\left(X\right)-B=0$E9

is known as the harmonic balance system, and the $n\left(2K+1\right)\times n\left(2K+1\right)$ composite conversion matrix

$J\left(X\right)=\frac{dH\left(X\right)}{dX}=j\Omega C\left(X\right)+G\left(X\right)$E10

is known as 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

$\Vert H\left(X^{\left[f\right]}\right)\Vert =\Vert j\Omega Q\left(X^{\left[f\right]}\right)+F\left(X^{\left[f\right]}\right)-B\Vert <\delta ,$E11

where δ is the allowed residual size, and $\Vert H(•)\Vert$ stands for some norm of the error function 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

$b\left(t\right)=e\left(t\right)\cos \left(2\pi f_{C}t\right),$E12

where e(t) is the envelope (a slowly varying signal), while $\cos \left(2\pi f_{C}t\right)$ is the fast-varying RF carrier. Simulating a circuit with this kind of stimulus, using conventional time-step integration schemes (e.g., Runge-Kutta schemes, or linear multistep methods), is computationally very expensive. The main reason is that the solution has to be computed during a long time interval imposed by the slowly varying envelope, whereas the step length is severely constrained by the high-frequency RF carrier.

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

$\widehat{b}\left(t_E,t_C\right)=e\left(t_E\right)\cos \left(2\pi f_C{t}_C\right),$E13

where t_E is the slow envelope time scale and t_C is the fast carrier time scale. In this particular case, $\widehat{b}\left(t_E,t_C\right)$ is a periodic function with respect to t_C but not to t_E, i.e.,

$\widehat{b}\left(t_E,t_C\right)=\widehat{b}\left(t_E,t_C+T_C\right),T_C=1/f_C.$E14

The univariate form, b(t), with $e\left(t\right)=2e^{-80\times {10}^{6}t}$ and $f_C=2\text{ GHz}$, is plotted in Figure 1 for the [0, 25 ns] time interval. The corresponding bivariate form, $\widehat{b}\left(t_E,t_C\right)$, is depicted in Figure 2 for the rectangular region [0, 25 ns] × [0, 0.5 ns]. There, it can be appreciated that $\widehat{b}\left(t_E,t_C\right)$ does not have as many undulations as b(t), thus allowing a more compact representation with fewer samples. Furthermore, due to the periodicity of $\widehat{b}\left(t_E,t_C\right)$ in 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.

Figure 1.

Envelope-modulated signal in the univariate time.

Figure 2.

Bivariate representation of the envelope-modulated signal.

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]:

$\frac{\partial q\left(\widehat{x}\left(t_E,t_C\right)\right)}{\partial t_E}+\frac{\partial q\left(\widehat{x}\left(t_E,t_C\right)\right)}{\partial t_C}+f\left(\widehat{x}\left(t_E,t_C\right)\right)=\widehat{b}\left(t_E,t_C\right).$E15

The mathematical relation between (1) and (15) establishes that if $\widehat{b}\left(t_E,t_C\right)$ and $\widehat{\mathit{x}}\left(t_E,t_C\right)$ satisfy (15), then the univariate forms $b\left(t\right)=\widehat{b}\left(t,t\right)$ and $x\left(t\right)=\widehat{x}\left(t,t\right)$ satisfy (1) [11]. Therefore, univariate solutions of (1) are available on diagonal lines $t_E=t,t_C=t,$ along the bivariate solutions of (15), i.e., x(t) may be retrieved from its bivariate form $\widehat{x}\left(t_E,t_C\right)$, by simply setting 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

$x\left(t\right)=\widehat{x}\left(t,t\mathrm{mod}{T}_C\right)$E16

on the rectangular domain $\left[0,t_{Final}\right]\times \left[0,T_C\right]$, where 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

$\begin{array}{c}\frac{\partial q\left(\widehat{x}\left(t_E,t_C\right)\right)}{\partial t_E}+\frac{\partial q\left(\widehat{x}\left(t_E,t_C\right)\right)}{\partial t_C}+f\left(\widehat{x}\left(t_E,t_C\right)\right)=\widehat{b}\left(t_E,t_C\right),\\ \widehat{x}\left(0,t_C\right)=i\left(t_C\right),\widehat{x}\left(t_E,0\right)=\widehat{x}\left(t_E,T_C\right),\end{array}$E17

on the rectangle $\left[0,t_{Final}\right]\times \left[0,T_C\right]$. i(·) is a given initial-condition function defined on [0,T_C], satisfying i(0)=i(_TC)=x(0), and $\widehat{\mathbf{x}}$(t_E0)=$\widehat{\mathbf{x}}$(t_E,T_C) is the periodic boundary condition due to the periodicity of the problem in the t_c fast carrier time scale.

The reason why bivariate envelope-modulated solutions do not need to be evaluated on the entire $\left[0,t_{Final}\right]\times \left[0,t_{Final}\right]$ domain (which would be computationally very expensive and would turn the multivariate strategy useless), and are restricted to the rectangle $\left[0,t_{Final}\right]\times \left[0,T_C\right]$, is because the solutions repeat along the 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 $\left[0,t_{Final}\right]\times \left[0,T_C\right]$ in the t_E slow time dimension described by the following general non uniform grid

$0=t_{E,0}<t_{E,1}<\cdots <t_{E,i-1}<t_{E,i}<\cdots <t_{E,K_E}=t_{Final},h_{E,i}=t_{E,i}-t_{E,i-1},$E18

in which 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

$\begin{array}{c}\frac{q\left({\widehat{x}}_i\left(t_C\right)\right)-q\left({\widehat{x}}_{i-1}\left(t_C\right)\right)}{h_{E,i}}+\frac{dq\left({\widehat{x}}_i\left(t_C\right)\right)}{d{t}_C}+f\left({\widehat{x}}_i\left(t_C\right)\right)=\widehat{b}\left(t_{E,i},t_C\right),\\ {\widehat{x}}_i\left(0\right)={\widehat{x}}_i\left(T_C\right),\end{array}$E19

where ${\widehat{\mathit{x}}}_i\left(t_C\right)$ is an approximation to the exact solution $\widehat{\mathbf{x}}$(t_E,i,t_C). Thus, once ${\widehat{\mathit{x}}}_{i-1}\left(t_C\right)$ is computed, the solution on the next slow time instant, ${\widehat{x}}_i\left(t_C\right)$, is evaluated by solving (19). Consequently, it is straightforward to conclude that we have to solve a set of K_E periodic boundary value problems if we want to obtain the solution $\widehat{x}\left(t_E,t_C\right)$ in the entire $\left[0,t_{Final}\right]\times \left[0,T_C\right]$ domain. With multitime ETHB, each one of these periodic boundary value problems is solved using the harmonic balance method. For each slow time instant t_E,i, the resultant HB system is the n(2K + 1) algebraic equation system defined by

$\frac{Q\left(\widehat{X}\left(t_{E,i}\right)\right)-Q\left(\widehat{X}\left(t_{E,i-1}\right)\right)}{h_{E,i}}+j\Omega Q\left(\widehat{X}\left(t_{E,i}\right)\right)+F\left(\widehat{X}\left(t_{E,i}\right)\right)=\widehat{B}\left(t_{E,i}\right),$E20

where $\widehat{B}\left(t_{E,i}\right)$ and $\widehat{X}\left(t_{E,i}\right)$ are the vectors containing the Fourier coefficients of the excitation sources and of the solution (the state variables), respectively, at t_E = t_E,i. $F(•)$ and $Q(•)$ are unknown functions that can be computed by evaluating f(·) and q(·) in the time domain and then calculating their Fourier coefficients. Ω is the diagonal matrix (6), and the $\widehat{X}\left(t_{E,i}\right)$ vector can be described as

$\widehat{X}\left(t_{E,i}\right)={\left[{\widehat{X}}_1{\left(t_{E,i}\right)}^T,{\widehat{X}}_2{\left(t_{E,i}\right)}^T,\dots ,{\widehat{X}}_n{\left(t_{E,i}\right)}^T\right]}^T,$E21

where each one of the state variable frequency components, ${\widehat{X}}_v\left(t_{E,i}\right)$, v = 1,…,n, is a (2K+1)×1 vector defined as

${\widehat{X}}_v\left(t_{E,i}\right)={\left[X_{v,-K}\left(t_{E,i}\right),\dots ,X_{v,0}\left(t_{E,i}\right),\dots ,X_{v,K}\left(t_{E,i}\right)\right]}^T.$E22

The HB system of (20) can be rewritten as

$H\left(\widehat{X}\left(t_{E,i}\right)\right)=\frac{Q\left(\widehat{X}\left(t_{E,i}\right)\right)-Q\left(\widehat{X}\left(t_{E,i-1}\right)\right)}{h_{E,i}}+j\Omega Q\left(\widehat{X}\left(t_{E,i}\right)\right)+F\left(\widehat{X}\left(t_{E,i}\right)\right)-\widehat{B}\left(t_{E,i}\right)=0,$E23

or simply as

$H\left(\widehat{X}\left(t_{E,i}\right)\right)=0,$E24

in which $H\left(\widehat{X}\left(t_{E,i}\right)\right)$ is the error function at t_E = t_E,i. In general, the nonlinear algebraic system of (24) is iteratively solved using Newton’s method

$H\left({\widehat{X}}^{\left[r\right]}\left(t_{E,i}\right)\right)+{\frac{dH\left(\widehat{X}\left(t_{E,i}\right)\right)}{d\widehat{X}\left(t_{E,i}\right)}|}_{\widehat{X}\left(t_{E,i}\right)={\widehat{X}}^{\left[r\right]}\left(t_{E,i}\right)}\left[{\widehat{X}}^{\left[r+1\right]}\left(t_{E,i}\right)-{\widehat{X}}^{\left[r\right]}\left(t_{E,i}\right)\right]=0,$E25

which requires that we have to solve a linear system of n(2K + 1) equations at each iteration r to compute the new estimate ${\widehat{X}}^{\left[r+1\right]}\left(t_{E,i}\right)$. Consecutive Newton iterations will be computed until a desired accuracy is achieved, i.e., until $\Vert H\left(\widehat{X}\left(t_{E,i}\right)\right)\Vert <\delta$, where δ is the allowed residual size.

The system of (25) requires the computation of the Jacobian matrix $J\left(\widehat{X}\left(t_{E,i}\right)\right)$, i.e., the derivative of the vector $H\left(\widehat{X}\left(t_{E,i}\right)\right)$, with respect to the vector $\widehat{X}\left(t_{E,i}\right)$,

$J\left(\widehat{X}\left(t_{E,i}\right)\right)=\frac{dH\left(\widehat{X}\left(t_{E,i}\right)\right)}{d\widehat{X}\left(t_{E,i}\right)}=\left[\begin{array}{cccc}\frac{\partial H_1\left(\widehat{X}\left(t_{E,i}\right)\right)}{\partial {\widehat{X}}_1\left(t_{E,i}\right)}& \frac{\partial H_1\left(\widehat{X}\left(t_{E,i}\right)\right)}{\partial {\widehat{X}}_2\left(t_{E,i}\right)}& \cdots & \frac{\partial H_1\left(\widehat{X}\left(t_{E,i}\right)\right)}{\partial {\widehat{X}}_n\left(t_{E,i}\right)}\\ \frac{\partial H_2\left(\widehat{X}\left(t_{E,i}\right)\right)}{\partial {\widehat{X}}_1\left(t_{Ei}\right)}& \frac{\partial H_2\left(\widehat{X}\left(t_{E,i}\right)\right)}{\partial {\widehat{X}}_2\left(t_{E,i}\right)}& \cdots & \frac{\partial H_2\left(\widehat{X}\left(t_{E,i}\right)\right)}{\partial {\widehat{X}}_n\left(t_{E,i}\right)}\\ \cdots & \cdots & \cdots & \cdots \\ \frac{\partial H_n\left(\widehat{X}\left(t_{E,i}\right)\right)}{\partial {\widehat{X}}_1\left(t_{E,i}\right)}& \frac{\partial H_n\left(\widehat{X}\left(t_{E,i}\right)\right)}{\partial {\widehat{X}}_2\left(t_{E,i}\right)}& \cdots & \frac{\partial H_n\left(\widehat{X}\left(t_{E,i}\right)\right)}{\partial {\widehat{X}}_n\left(t_{E,i}\right)}\end{array}\right]$E26

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

$\frac{d{H}_m\left(\widehat{X}\left(t_{E,i}\right)\right)}{d{\widehat{X}}_l\left(t_{E,i}\right)}=\frac{1}{h_{E,i}}\frac{d{Q}_m\left(\widehat{X}\left(t_{E,i}\right)\right)}{d{\widehat{X}}_l\left(t_{E,i}\right)}+j\Omega \frac{d{Q}_m\left(\widehat{X}\left(t_{E,i}\right)\right)}{d{\widehat{X}}_l\left(t_{E,i}\right)}+\frac{d{F}_m\left(\widehat{X}\left(t_{E,i}\right)\right)}{d{\widehat{X}}_l\left(t_{E,i}\right)}.$E27

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

$\begin{array}{l}\frac{\partial q_A\left({\widehat{x}}_A\left(t_E,t_C\right),{\widehat{x}}_L\left(t_E,t_C\right)\right)}{\partial t_E}+\frac{\partial q_A\left({\widehat{x}}_A\left(t_E,t_C\right),{\widehat{x}}_L\left(t_E,t_C\right)\right)}{\partial t_C}+f_A\left({\widehat{x}}_A\left(t_E,t_C\right),{\widehat{x}}_L\left(t_E,t_C\right)\right)=\widehat{b}\left(t_E,t_C\right)\\ \frac{\partial q_L\left({\widehat{x}}_A\left(t_E,t_C\right),{\widehat{x}}_L\left(t_E,t_C\right)\right)}{\partial t_E}+\frac{\partial q_L\left({\widehat{x}}_A\left(t_E,t_C\right),{\widehat{x}}_L\left(t_E,t_C\right)\right)}{\partial t_C}+f_L\left({\widehat{x}}_A\left(t_E,t_C\right),{\widehat{x}}_L\left(t_E,t_C\right)\right)=\widehat{b}\left(t_E,t_C\right)\end{array}$E28

with

$x\left(t\right)=\widehat{x}\left(t,t\right)=\left[\begin{array}{c}{x}_A\left(t\right)={\widehat{x}}_A\left(t,t\right)\\ x_L\left(t\right)={\widehat{x}}_L\left(t,t\right)\end{array}\right],x_A\left(t\right)\in {ℝ}^{n_A},x_L\left(t\right)\in {ℝ}^{n_L},n_A+n_L=n,$E29

where 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

$x_A\left(t\right)=x_{A,v}\left(t\right)={\displaystyle \sum _{k=-K}^K{X}_{k,v}\left(t\right)e^{jk2\pi f_{C}t}},v=1,\dots ,n_A$E30

while state variables in the baseband block can be seen as slowly varying aperiodic functions of the form

$x_L\left(t\right)=x_{L,v}\left(t\right)={\psi }_v\left(t\right),v=1,\dots ,n_L.$E31

In (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 ${\widehat{x}}_{A,v}\left(t_E,t_C\right)$ and ${\widehat{x}}_{L,v}\left(t_E,t_C\right)$, and defined as

${\widehat{x}}_{A,v}\left(t_E,t_C\right)={\displaystyle \sum _{k=-K}^K{X}_{k,v}\left(t_E\right)}{e}^{jk2\pi f_C{t}_C}$E32

and

${\widehat{x}}_{L,v}\left(t_E,t_C\right)={\psi }_v\left(t_E\right),$E33

where t_E and t_C are, respectively, the slow envelope time dimension and the fast carrier time dimension. As can be seen, since the ${\widehat{x}}_{L,v}\left(t_E,t_C\right)$ state variables have no dependence on 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 ${\widehat{x}}_{L,v}\left(t_{E,i},t_C\right)$ is merely a constant 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 ${\widehat{x}}_{L,v}\left(t_{E,i},t_C\right)$, which means that these state variables can be processed in a purely time-domain scheme. In contrast, for each slow time instant t_Ei, each of the ${\widehat{x}}_{A,v}\left(t_{E,i},t_C\right)$ is a waveform that has to be represented as a Fourier series adopting a convenient harmonic truncation at some order k = −K,…,K, i.e., each of the ${\widehat{x}}_{A,v}\left(t_{E,i},t_C\right)$ is a waveform that requires a total of 2K + 1 harmonic components for a convenient frequency domain representation. In summary, while active state variables have to be represented by a set of 2K + 1 Fourier coefficients arranged in (2K + 1)×1 vectors of the form

${\widehat{X}}_{A,v}\left(t_{E,i}\right)={\left[X_{A,v,-K}\left(t_{E,i}\right),\dots ,X_{A,v,0}\left(t_{E,i}\right),\dots ,X_{A,v,K}\left(t_{E,i}\right)\right]}^T,v=1,\dots ,n_A,$E34

latent state variables can be represented as 1×1 scalar quantities, i.e., they can be simply represented as

${\widehat{X}}_{L,v}\left(t_{E,i}\right)=X_{L,v,0}\left(t_{E,i}\right),v=1,\dots ,n_L.$E35

By considering this, we can easily deduce that the size of the $\widehat{X}\left(t_{E,i}\right)$ vector defined by (21) will be significantly decreased, as well as the total number of equations in the HB system of (23). Furthermore, another crucial aspect is that we are no longer forced to carry out the conversion between time and frequency domains for the latent state variables expressed in the form of (35), as well as for the components of $H\left(\widehat{X}\left(t_{E,i}\right)\right)$ corresponding to latent blocks of the circuit. Given that the k = 0 order Fourier coefficient X_v,0(t_E,i) is exactly the same as the constant t_C time value ${\widehat{x}}_v\left(t_{E,i}\right)$, components of the HB system of (23) that have no dependence on active state variables will not be required for any direct or inverse Fourier transformation operations.

Considerable Jacobian $J\left(\widehat{X}\left(t_{E,i}\right)\right)$ matrix size reductions will also be achieved with this technique. Indeed, by considering the latency of state variables in some parts of the circuit, some blocks of the Jacobian matrix (26) are simply reduced to 1×1 scalar elements. These scalar elements contain the dc sensitivity of $H\left(\widehat{X}\left(t_{E,i}\right)\right)$ to the latent components of $\widehat{X}\left(t_{E,i}\right)$.

With the state variable $\widehat{X}\left(t_{E,i}\right)$ and the error function $H\left(\widehat{X}\left(t_{E,i}\right)\right)$ vector size reductions, as also the resulting Jacobian $J\left(\widehat{X}\left(t_{E,i}\right)\right)$ matrix size reduction, it is possible to avoid dealing with large linear systems in the iterations of (25). Thus, a less computationally expensive Newton-Raphson iterative solver is required to solve (23). In conclusion, partitioning the circuit into active and latent sub-circuits (blocks) can lead us to significant computation and memory savings when computing the solution.

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⁻⁷ for all the state variables of the circuit.

Figure 3.

Simplified schematic diagram of a mixer (frequency translation device).

Simulationtimeinterval	Partitionedtime-frequencytechnique	MultitimeETHB	Univariate time-stepintegration(Gear-2 method)
[0, 1.0 μs]	00:00:04.9	00:00:11.3	00:19:21
[0, 10.0 μs]	00:00:39.6	00:01:35.1	02:47:33

Table 1.

CPU time (h:min:sec)—simulation of the circuit depicted in Figure 3.

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.

Acknowledgments

This work is funded by National Funds through FCT - Fundação para a Ciência e a Tecnologia, under the project UID/EEA/50008/2013.