Method and apparatus for simulating quasi-periodic circuit operating conditions using a mixed frequency/time algorithm

Described is a process for performing an improved mixed frequency-time algorithm to simulate responses of a circuit that receives a periodic sample signal and at least one information signal. The process selects a set of evenly spaced distinct time points and a set of reference time points. Each of the reference points is associated with a distinct time point, and a reference time point is a signal period away from its respective distinct time point. The process finds a first set of relationships between the values at the distinct time points and the values at the reference time points. The process also finds a second set of relationships between the values at the distinct time points and the values at the reference time points. The process then combines the first and second sets of relationships to establish a system of nonlinear equations in terms of the values at the distinct time points only. By solving the system of nonlinear equations, the process finds simulated responses of the circuit in time domain. The process then converts the simulated circuit responses from time domain to frequency domain.

FIELD OF THE INVENTION

The present invention relates generally to analog circuit design simulations, and more specifically to analog circuit design simulations using a mixed frequency/time approach.

BACKGROUND OF THE INVENTION

Using a description language such as a netlist and device models an analog circuit can be first designed in terms of its predetermined inputs and expected outputs. The analog circuit design is then simulated before it is physically fabricated on a silicon chip.

One of the most difficult challenges in analog circuit stimulation is the analysis of the circuits that operate on multiple time scales. Typical examples of this type of circuits are switched-capacitor filters and circuits used in RF (radio frequency) communications systems. Applying standard transient analysis to a circuit of this type requires simulation of the detailed responses of the circuit over hundreds of thousands of clock cycles (millions of time points).

Many circuits of engineering interest are designed to operate near a time-varying, but quasi-periodic, operating point. Some of these circuits can be analyzed under the assumption that one of the circuit inputs produces a periodic response that can be directly calculated by steady-state algorithms, thus avoiding long transient simulation times. Under this assumption, all other time-varying circuit inputs are treated as small-signal by linearizing the circuit around the periodic operating point.

Existing algorithms are able to find periodic operating points and to perform periodic time-varying small-signal analysis. However, many circuits cannot be analyzed with the periodic-operating-point-plus-small-signal approach, because the above-described assumption may not apply. For example, predicting intermodulation distortion of a narrowband circuit, such as a mixer-plus-filter circuit, involves calculating the nonlinear response of the mixer circuit, driven by an LO (local oscillator), to two high-frequency inputs that are closely spaced in frequency. The steady-state response of such a circuit is quasi-periodic.

The analog circuit simulation is further complicated by the fact that many multi-timescale circuits have a response (again mixers and switched-capacitor filters are typical examples) that is highly nonlinear with respect to at least one of the exciting inputs, and so steady-state approaches, such as the multi-frequency harmonic balance approach, do not perform well. To circumvent these difficulties, mixed frequency-time (MFT) algorithms have been proposed. Specifically, the MFT algorithms exploit the fact that many circuits of engineering interest have a strongly nonlinear response to only one input, such as the clock in the case of a switched-capacitor circuit, or local oscillator in the case of a mixer, but respond only in a weakly nonlinear manner to other inputs.

Unfortunately, existing MFT algorithms suffer from several drawbacks that prevent their application to practical circuits, particularly large circuits. In existing MFT algorithms, poor sample point selection leads to ill-conditioned simulation environment, in which simulation values may be unsolveable with acceptable accuracy. In addition, existing MFT algorithms are based on a matrix-explicit linear solver (via Gaussian elimination) whose computational cost (or time) is proportional to an order of N3for each Newton iteration, where N is the number of nodes of the circuit in simulation.

A new class of algorithms have been developed for simulating multi-timescale circuits by converting the circuit DAE (differential-algebraic equation) into an equivalent multi-variable partial differential equations (M-PDE). However, the effectiveness of the M-PDE method to simulate large circuits has yet to be proven. In addition, there is evidence that, for some circuits, the M-PDE method generates inaccurate simulation results.

There is, therefore, a need in the art for a method and apparatus that utilizes the MFT method to accurately simulate large circuits.

There is another need in the art for a method and apparatus utilizing the MFT method to simulate large circuits with reduced computational cost and increased speed.

There is still another need in the art for a method and apparatus for generating an efficient linear problem solver structured such that the MFT method can accurately simulate large circuits with improved convergence.

The present invention provides a method and apparatus to meet these and other needs.

SUMMARY OF THE INVENTION

To overcome the shortcomings of the available art, the present invention discloses a novel method and apparatus for simulating analog circuits by using a mixed frequency/time approach.

In broad terms, the present invention provides a method for simulating responses of a circuit, the circuit receiving a periodic sample signal and at least one information signal. The method comprises the steps of: selecting a set of distinct time points; defining a set of reference time points, wherein each of the reference time points is associated with one of the distinct time pints; establishing a first set of relationships between the values at the distinct time points and the values at the reference time points; establishing a second set of relationships between the values at the distinct time points and the values at the reference points; combining the first and second relationships to establishing a system of equations in terms of the values at the distinct time points; and finding responses of the circuit at the distinct time points by solving the established system of equations.

The present invention also provides a corresponding apparatus for performing steps in the method described above.

DETAILED DESCRIPTION OF PREFERRED EMBODIMENTS

FIG. 1shows circuit block102including N circuit nodes, which will be used to illustrate a simulation process in accordance with the present invention. Circuit102receives S information signals I1, I2, . . . , Is from inputs104.1,104.2, . . . ,104.s, respectively. Circuit102also receives a clock (or reference) signal C from input106. In response to receipt of the input signals, circuit102generates a signal at output108.

1. The Improved MFT Algorithm

The behavior of a circuit (such as circuit102inFIG. 1) can be described by a set of nonlinear differential-algebraic equations (DAEs) that can be written as

ⅆⅆt⁢Q⁡(v⁡(t))+I⁡(v⁡(t))+u⁡(t)=0,(1)
Where Q(v(t))∈Nis typically the vector of sums of capacitor charges at each node, I(v(t))∈Nis the vector of sums of resistive currents at each node, u(t)∈Nis the vector of inputs, v(t)∈Nis the vector of node voltages, and N is the number of circuit nodes.

The present invention is particularly advantageous in situations where the input signal u(t) is quasiperiodic. A signal is L-quasiperiodic if it can be written as a Fourier series with L fundamental frequencies. RF circuits are generally influenced by one periodic timing signal, often referred to as the LO (local oscillator) or the clock, and one or more information signals. If fcdenotes the clock signal frequency (such as the signal at input106shown inFIG. 1), and f1. . . fsare the S information signals (shown inFIG. 1), then the (S+1)-quasiperiodic input can be written as:

In a preferred embodiment, the present invention utilizes two conditions to improve the simulation of quasi-periodic circuit operating conditions. The first condition is that the circuit of interest possesses a quasiperiodic steady-state response. That is, v(t) is an S+1 quasiperiodic signal with fundamentals f1, . . . fs, fc. The second condition is that all physical circuits have a finite bandwidth. Using these two conditions, the present invention selects only a finite number of Fourier series terms to approximate v(t) while maintaining the necessary accuracy. Thus:

Assume that v(t) is sampled at a discrete set of points t′n=t0+nTc, where Tc=1/fcis the clock period, t0∈[0,Tc) and n runs over the integers, to obtain a discrete signalv(t′). Since

FIG. 2shows an exemplary envelop in response to the signals at input104.1(assuming that the inputs at104.2, . . . ,104.S are zero inputs) and input106. InFIG. 2, the hollow, circular dots represent samples (or distinct points), while the continuous waveform (dashed line) is the waveform havingVas its Fourier coefficients, or equivalently, obtained by Fourier interpolation of the sample points.

In principle, because there are only

K=∏s=1s⁢(2⁢K3+1)
Fourier coefficients to representv, once the value of K distinct points t1, . . . ,tkalong the sample envelope are known, then the full envelope can be recovered. The envelope corresponding to the quasiperiodic operating point is obtained by obtaining K sample points that lie on the solution to the DAE given by equation (1).

FIG. 3shows a particular scheme of selecting K (K=2) sample points along the sample envelope shown inFIG. 3to obtain five distinct points (or distinct time points) [5=(K×2)+1], in accordance with the present invention. InFIG. 3, the distinct points are denoted by circles, and the reference points (reference time points) are denoted by squares.

Lets us define the state transition function φ (v0, tk, tf)=v(tf)=v(tf):v(t) that satisfies equation (1) for t∈[tk,tf] and v(tk)=v0. In particular, define the vector
v0=[vT(t1), . . . ,vT(tK)]T=[vT(t1), . . . ,vT(tK)]T,  (6)
where superscript T denotes matrix transpose, to containvat the K sample points tk=t0+nkTc, k=1 . . . K, nk∈Z. The value of the K points that follow by one cycle can be obtained from the transition function,

v_TcT=[vT⁡(t1+Tc),…⁢,vT⁡(tK+Tc)]T=[ϕ⁡(v⁡(t1),t1,t1+Tc)T,…⁢,ϕ⁡(v⁡(tK),tK+Tc)T]T(7)
which may be written more compactly by introducing the multi-cycle transition function that is the collection of the K transition functions from tkto tk+Tc, as
vTc=ΦTTc(v0).  (8)

Note that for each mode n, the vector of signals on that node, at the sample time plus one clock cycle,

v_Tn,
is a delayed version of the signals at the sample points (this will be discussed in greater detail below). there exists, therefore, a linear operator DTcthat maps

Note that DTcis a real matrix and independent of node n. Hence equation (9) holds for each n=1, . . . , N, and represents a boundary condition on a solution to (1).

Combining equations (8) and (9) gives
(DTc{circumflex over (x)}IN)v0−ΦTc(v0)=0,  (10)
where {circle around (x)} is the Kronecker product1and INis the N by N identity matrix. Equation (10) is a system of KN nonlinear equations and KN unknowns that can be solved for the envelope sample points. From these sample points and the transition functions, the circuit's quasiperiodic operating point (in particular, the spectrum of v) can be recovered.

2. Sample Point Selection

To contruct the matrix DTc, referred to as the delay matrix, consider the Fourier series ofv0andvTc. Referring to equation (4), equation (11) holds

Thus if Γ is the matrix mapping sample points on the envelope to Fourier coefficients, then the delay matrix may be constructed as
DTc=Γ−1ΩTcΓ.  (13)
In particular Γ may be constructed as the Kronecker product of one-dimensional (2Ks+1)-point Fourier-transform matrices
Γmn(s)=ej2πmfstn/(2Ks+1)(14)
as
Γ=Γ(1){circumflex over (x)} . . . {circumflex over (x)}Γ(S)(15)
From the properties of Kronecker products, Γ−1is likewise a Kronecker product of the inverses of the Γ(s). In the existing MFT algorithms, no particular consideration was given to the choice of the sample points tk, so that the Γ(s)'s there are ill-conditioned matrices corresponding to an “almost-periodic” Fourier transform. By contrast, the improved MFT algorithm of the present invention performs a process of choosing well-conditioned sample points.

Assume the K sample points can be arranged into an S-dimensional array τ(k1, . . . ,kS), −Ks≦ks≦Ks, 1≦s≦S, such that for a given dimension s, there exists an integer p, and

τ⁡(…⁢,ks+1,…)-τ⁡(…⁢,ks,…)=Ts2⁢Ks+1+pTs(16)
holds. In this case, the entries of the Γ(s)matrices are:

Γmn(s)=ⅇj⁢⁢2⁢⁢π⁢⁢mn/(2⁢Ks+1)(17)
That is, they are the DFT matrices, and the matrix Γ:C2K1+1x . . . x C2Ks+1C2K1+1x . . . x C2Ks−1represents an S-dimensional DFT. Thus Γ has a condition number of one; it is perfectly well-conditioned.

The Newton's method can now be employed to solve
(DTc{circumflex over (x)}IN)v0−ΦTc(v0)=0,  (18)
At iteration i, the Jacobian matrix is given by

v_0-∂Φ∂v_0v_0=v_0i,(19)
Recall from (13) DTc=Γ−1ΩTc, which is fixed through all Newton iterations. Let

J=∂ϕ∂v_0|v_0=v_0i
be obtained from the multicycle transition function by

∂Φ∂v_0v_0=v_0i=[∂ϕ1∂v_0(t1⁢)v_0⁡(t1)=v_0⁡(t1)i⋰∂ϕK∂v_0(tK⁢)v_0⁡(tK)=v_0⁡(tK)i](20)
Note that J is block-diagonal. Defining

b=-(DTc⊗IN)⁢v_0i-ϕ⁡(v_0i),
the Newton iteration is performed by solving the equation

((DTc⊗IN)-J)⁢Δ⁡(v_0i)=b,(21)
using an iterative Generalized Minimal Residual (GMRES) solver, and setting

v_0i+1=v0i+Δ⁢v_0i.(22)
Each iteration of GMRES requires a matrix-vector multiplication. For a vector q∈KN, the term (DTc{circle around (x)} IN)q is calculated by first applying a K dimensional DFT N times, then scaling each row with ΩTc, and finally applying a K dimensional inverse DFT N times.

Let q be partitioned into q=[q1T, . . . ,qKT]T,qk∈N, for 1≦k≦K. Then

∂Φ∂v_⁢q=[∂ϕ1∂v_0⁡(t1)⁢q⁢⁢1⋮∂ϕK∂v_0⁡(tK)⁢qK].(23)
The calculation of each

∂ϕk∂v_0⁡(tk)⁢qk
can be carried through matrix-vector multiplication and backsolving without explicitly forming the

For many problems, the GMRES algorithm is not efficient for solving equation (21) without an effective preconditioner. To analyze the reason, consider the case where the state transition function of the circuit, over one clock cycle, is approximately linear, that is φ(χ,t,t+Tc)≃Hχ(t). Linear circuits are an obvious example of a case where this is true, and while nonlinear circuits will have nonlinear state-transition functions, if the method performs poorly for linear circuits it surely will not work well for nonlinear circuits either. However, many nonlinear circuits have a state-transition function that is nearly linear, a fact which is exploited below to construct an effective preconditioner. The convergence of the GMRES method will depend on the location of the eigenvalues of the Jacobian matrix, DTc−J. If λHis an eigenvalue of the matrix H, then eiωtc−λH, where ω=2π(k1f1+k2f2+ . . . ksfs) will be an eigenvalue of DTc−J, for every k1,k2, . . . in the MFT analysis. Thus unless all the secondary input frequencies are nearly commensurate with the clock frequency, the eigenvalues of DTc−J will be “fanned out” by delay matrix. This will cause severe convergence problems for the GMRES solver. Roughly speaking, the GMRES algorithm in the MFT algorithm with K total harmonics will take K times as many iterations to coverage than the GMRES iteration for the steady-state problem with only the clock excitation applied. This follows because the eigenvalues of H are typically inside the unit circle of the complex plane. The delay matrix replicates the eigenvalue structure K times, each shift being a complex number of order unity, and generally causing the convex hull of the eigenvalues of DTc−J to enclose the origin.

The following lemmas about the properties of Kronecker products are needed to perform the formal analysis.

If λHis an eigenvalue of ∂Φ/∂v0, then ej2πωkTcis an eigenvalue of the MFT Jacobian matrix.

The proof is as follows. For linear circuits, the diagonal blocks of

∂ϕK∂v_0⁡(t1).
Denote a diagonal block as H, then the Jacobian matrix is equal to

⁢(Γ-1⁢ΩTc⁢Γ)⊗IN-(IK⊗H)(24)=⁢(Γ-1⊗IN)⁢(ΩTc⊗IN)⁢(Γ⊗IN)-(IK⊗H)(25)=⁢(Γ-1⊗IN)⁢{(ΩTc⊗IN)-(Γ-1⊗IN)-1⁢(IK⊗H)⁢(Γ⊗IN)-1}⁢(Γ⊗IN)(26)=⁢(Γ-1⊗IN)⁢{(ΩTc⊗IN)-(Γ⊗IN)⁢(IK⊗H)⁢(Γ⊗IN)-1}⁢(Γ⊗IN)(27)=⁢(Γ-1⊗IN)⁢{(ΩTc⊗IN)-(IK⊗H)⁢(Γ⊗IN)⁢(Γ⊗IN)-1}⁢(Γ⊗IN)(28)=⁢(Γ-1⊗IN)⁢{(ΩTc⊗IN)-(IK⊗H)}⁢(Γ⊗IN).(29)
Equation (24) to equation (25) holds because of IN=INININand Lemma 1. Equation (26) to equation (27) holds due to Lemma 2(b), and equation (27) to equation (28) holds due to Lemma 5.2(a). Since (Γ−1{circle around (x)} IN) is unitary and its inverse is (Γ {circle around (x)} IN)−1, the right hand side of equation (29) has the same spectrum as (ΩTc{circle around (x)} IN)−(IK{circle around (x)} H). It is easy to verify that (ΩTc{circle around (x)} IN)−(IK{circle around (x)} H) is block diagonal, hence its eigenvalues are the union of eigenvalues of all the blocks, ej2πωkTcIN−H, for k=1, . . . ,K.

The preceding analysis suggests a good way of preconditioning for solving the Newton equation (21). Solving equation (21) is equivalent to solving
{ΩTc−((Γ{circle around (x)} IN)JΦ(Γ−1{circle around (x)} IN))}γ=Γb,(30)
where γ=ΓΔvi. A good choice of preconditioner is P=(ΩTc{circle around (x)} IN)−(IK{circle around (x)} H) where H can be chosen as the Jacobian matrix from the steady-state analysis in the initial guess stage discussed in Section 5, or any of the diagonal blocks

∂ϕ∂v_0.
In particular, if the signal-cycle state-transition function is linear and time invariant, then the Newton equation can be solved in a single GMRES iteration. Note that the preconditioner presented here is effective if the Jacobian of the state-transition function is nearly constant over multiple cycles. The circuit behavior inside each clock cycle is hidden from the preconditioner. This is not the case in, for example, the time- or frequency-averaged preconditioners typically used in modern harmonic balance codes. For this reason the preconditioner presented here may perform well under much weaker assumptions about the circuit behavior, in particular at higher power levels.

For each GMRES iteration, a system Pu=v has to be solved. Since P is block diagonal, it needs to solve a sequence of K systems (ej2πωkTcIN−H)uk=vk, for k=1, . . . , K, where uT=[u1T, . . . ,uKT]Tand vT=[v1T, . . . ,vKT]T. The preconditioner can be applied very efficiently by incorporating a Krylov subspace reuse algorithm, as the linear system to be solved are the same as arise in the small-signal analysis for periodic time-varying systems. The basic idea of the algorithm is that the Krylov subspaces associated with the matrices ejωTc−H are very similar for different ωk. Essentially, the Krylov-subspace re-use algorithm allows the preconditioner for the matrix (DTc{circle around (x)} IN)−J to be applied with only slightly more cost than an iterative solve with the matrix H.

FIG. 4shows the effectiveness of the preconditioning process in compressing the eigenvalues for an RF receiver circuit. Specifically,FIG. 4shows eigenvalue distribution before and after the preconditioning process. InFIG. 4, the eigenvalues are very tightly clustered around unity, indicating excellent performance of the preconditioner and very rapid GMRES convergence.

FIG. 5shows the effectiveness of the preconditioning process in reducing the number of GMRES iterations needed to solve each MT Newton update equation, for the same RF circuit mentioned above. Only three iterations are needed to reduce the residual by a factor of 10−2, whereas without the preconditioning process, over 400 iterations are needed to achieve any reduction in the residual at all. Since the MFT circuit equations are not solved exactly, on average there is a performance advantage in the MFT method to using approximate solves of the Newton update equation, and therefore GMRES is converged to a relatively loose tolerance.

5. Improving Newton Convergence

Rapid convergence of Newton's method can only be assured with a good initial estimation. To achieve a good initial estimation, the present invention first calculates the periodic steady state response of the circuit with the clock signal applied, while suppressing other non-DC signals. Using the steady state solution as an operating point, a small-signal analysis is performed by treating non-clock fundamentals as small signals. As a result of the small signal analysis, amplitudes at fs+ksfc, for −Ks≦Ks≦Ks, 1≦s≦S, are generated. These amplitudes are transformed into time domain initial conditions via inverse multidimensional discrete Fourier transform (DFT). At higher input power levels, using a Newton continuation method, with the amplitude of the non-clock signals as the continuation parameter, is generally effective in securing convergence.

After the solution is converged, the valuesv=[v(t1)T, v(t2)T, . . . ,v(tK)T]Tand the integration solution inv=[v(t1)T, v(t2)T, . . . ,v(tK)T]Tare available. From these pieces of information, the spectrum v(t) can be obtained. Let

v⁡(.,kc)=1Tc⁢∫0Tc⁢{(Ω⁡(τ)-1⁢Γ)⊗IN}⁢v_⁡(τ)⁢ⅇj⁢⁢2⁢⁢π⁢⁢kc⁢fc⁢τ⁢ⅆτ(33)
Forming {(Ω(T)−1Γ){circle around (x)}IN}v(τ) requires the values for v(t1+τ), . . . ,v(tk+τ), or synchronized time steps between cycles. The total computational costs is one KN-vector integration and M Fourier transforms, where M is the number of synchronized time points.

The synchronized time step requirement may not be easily met in practice. One alternative is to use interpolation schemes. However, these schemes potentially lose accuracy. Another alternative is to utilize integration instead of multidimensional discrete Fourier transforms. Specifically, it is easy to verify that

7. Simulation Utilizing A Preferred Embodiment Of MFT Method

The first example is a low-pass switched-capacitor filter of 4 kHz bandwidth and having 238 nodes, resulting in 337 equations. To analyze this circuit, the improved MFT of the present invention analysis was performed with an 8-phase 100 kHz clock and a 1V sinusoidal input at 100 Hz.

The 1000 to 1 clock to signal ratio makes this circuit difficult for traditional circuit simulators to analyze. In the improved MFT method, three harmonics were used to model the input signal. The eight-phase clock resulted in the need to use about 1250 timepoints in each transient integration. This brings the total number of variables solved by the analysis to slightly less than three million (337×(2×3+1)×1250=2,948,750). The simulation took a little less than 20 minutes CPU time to finish, on a Sun UltraSparc1 workstation with 128 Megabyte memory and a 167 MHz CPU clock.FIG. 6shows the output spectrum of the filter.

The second example is a high-performance image rejection receiver. It consists of a low-noise amplifier, a splitting network, two double-balanced mixers, and two broadband Hilbert transform output filters combined with a summing network that is used to suppress the undesired side-band. A limiter in the LO path is used for controlling the amplitude of the LO. It is a rather large RF circuit that contains 167 bipolar transistors and uses 378 nodes. This circuit generates 987 equations in the simulator.

To determine the intermodulation distortion characteristics, the circuit was driven by a 780 MHz LO and two 50 mV closely placed RF inputs, at 840 MHz and 840 MHz+10 KHz, respectively. Three harmonics were used to model each of the RF signals. 200 time points were used in each transient clock-cycle integration, considered to be conservative in terms of accuracy for this circuit. As a result, nearly ten million unknowns (987×(2×3+1)2×200=9,672,600) were generated. It took 55 CPU minutes to finish on a Sun UltraSpare10 workstation with 128 Megabytes of physical memory and a 300 MHz CPU clock.FIG. 7shows 3rd and 5th order distortion products.

To understand the efficiency of the improved MFT method of the present invention, consider that traditional transient analysis would need at least 80,000 cycles of the LO to compute the distortion, a simulation time of over two days. In contrast, the MFT method of the present invention is able to resolve very small signal levels, such as the 5th order distortion products show inFIG. 7.

Solving the MFT equations by direct factorization methods is also impractical, as the storage needed for the factored rank—50,000(987×(2×3+1)2=48,363) MFT Jacobian of Equation 19 is several gigabytes. Forming the Jacobian matrix by direct methods would also require computation time proportional to the cost of 50,000 transient integration cycles, again a number on the order of days.

FIG. 8shows a flowchart illustrating a process of simulating the responses of a circuit such as the circuit shown inFIG. 1, in accordance with the a preferred embodiment of the present invention.

Step804selects a set of evenly spaced distinct time points shown as the circle dots inFIG. 3. The details of step804can be found in equation (16) and related descriptions.

Step806defines a set of reference time points shown as the square dots inFIG. 3. As shown inFIG. 3, each of the reference time points is associated with a respective distinct time point, and each of the reference time points is one signal period (or clock cycle) away from its respective distinct time point.

Step808establishes a first set of relationships between the values at the distinct time points and the values at the reference time points. The details of step808can be found in equation (8) and related descriptions.

Step810establishes a second set of relationships between the values at the distinct time points and the values at the reference time points. The details of step810can be found in equation (9) and related descriptions.

Step812combines the first and second sets of relationships to establish a system and equations that contain the values at the distinct time points only. The details of step808can be found in equations (10) and (18) and related descriptions.

Step814finds (or generates) the simulated responses of the circuit at the distinct time points by solving the established system of equations. If a circuit includes N internal circuit nodes and M outputs, step814can find (or generate) the simulated responses for all of the N internal circuit nodes and M outputs. The details of step814can be found in equations (18)-(22) and (30) and related descriptions.

FIG. 9shows a flowchart illustrating an exemplary process of solving the established system of equations in step814ofFIG. 8, in accordance with the present invention.

Step904selects a set of estimated values to reflect estimated circuit responses at the distinct time points. The details of step904can be found in Section 5.

Step906establishes a system of linear equations at the estimated values. The details of step906can be found in equation (21) and related descriptions.

Step908preconditions the system of linear equations to improve the convergence of solution to the system of linear equations. The details of step908can be found Section 4.

Step910solves the system of linear equations to generate the correction values to adjust the estimated circuit responses at the distinct time points. The details of step910can be found in equations (21)-(22) and related descriptions.

Step912adjusts the estimated values as newly estimated values to reflect the estimated circuit responses at the distinct time points. The details of step912can be found in equations (21)-(22) and related descriptions.

Step914determines whether the adjusted estimated values have an acceptable accuracy to represent the circuit responses. If the determination is negative, the process is led to step906. If the determination is positive, the process is led to step916. The estimated values and adjusted estimated values are in time domain.

Step916converts the estimated values from time domain to frequency domain. The details of step916can be found in Section 6, equations (31)-(34).

8. Hardware Platform

FIG. 10shows a block diagram illustrating an exemplary computer system1000, which can be used as a hardware platform for executing the program that performs the processes shown inFIGS. 8 and 9.

The hard disk1008is coupled to disk drive interface1006; monitor display1012is coupled to display interface1010; and mouse1016and keyboard1018are coupled to bus interface1014. Coupled to system bus1001are processing unit1002, memory device1004, disk drive interface1006, display interface1010, and network communication interface1020.

The memory device1004stores data and programs. Operating together with disk drive interface1006, hard disk1008also stores data and programs. However, memory device1004has faster access speed than hard disk1008, while hard disk1008has higher capacity than memory device1004.

Operating together with the display interface1010, display monitor1012provides visual interfaces between the programs being executed and users, and displays the outputs generated by the programs.

Operating together with bus interface1014, mouse1016and keyboard1018provide inputs to computer system1000.

The network communication interface1020provides an interface between computer system1000and network104in accordance with predetermined networking protocols.

The processing unit1002, which may include more than one processor, controls the operations of computer system1000by executing the programs stored in memory device1004and hard disk1008. The processing unit also controls the transmissions of data and programs between memory device1004and hard disk1008.

In the present invention, the program for performing the steps shown inFIGS. 8 and 9can be stored in memory device1014or hard disk1018. The program can be executed by processing unit1002.

The present invention improves the existing MFT method. The MFT method of the present invention is an efficient approach to analyzing multi-frequency nonlinear effects such as intermodulation distortion. Making the MFT method computationally efficient on problems of engineering interest required careful construction of the delay matrix, matrix-implicitly Krylov subspace iterative linear solvers, and a preconditioner tailored to the MFT method and the circuits it typically analyzes. As a result, nonlinear systems comprising tens of millions of unknowns can be solved in less than an hour with computational resources commonly available to engineering designers.

One salient advantage of the MFT method in the present invention is in computing the functions Φ and the product of the Jacobian of Φ with some vectors. Both computations are essentially the solution of an initial value problem. Each application of the operator DTc−J or calculation of the Newton residual, involves solving K such initial value problems, that is, integrating K sets of DAEs forward in time over one clock period. Each of the K problems, however, is essentially decoupled. Parallel implementations of the MFT will therefore enjoy very efficient processor utilization. This decoupling also assists the implementation of out-of-core solvers. In fact, it has been observed it is possible to implement the MFT algorithm as an out-of-core algorithm with over 80% average CPU utilization.

While the invention has been illustrated and described in detail in the drawing and foregoing description, it should be understood that the invention may be implemented through alternative embodiments within the spirit of the present invention. Thus, the scope of the present invention is not intended to be limited to the illustration in this specification, but is to be defined by the appended claims.