WO1999036872A1

WO1999036872A1 - Method for simulating circuits using harmonic balance, and apparatus therefore

Info

Publication number: WO1999036872A1
Application number: PCT/RU1998/000323
Authority: WO
Inventors: Sergej Rusakov; Mark Gourary; Sergej Ulyanov; Michael Zarov; Kiran Gullapalli; Brian Mulvaney
Original assignee: Motorola Inc.
Priority date: 1998-10-13
Filing date: 1998-10-13
Publication date: 1999-07-22
Also published as: WO1999036872A8

Abstract

A simulator (100) uses a modified harmonic balance method (400) for analyzing circuit behavior. After receiving a node list (412), providing an initial vector (414) in a frequency domain, and formulating an equation system (416) with a Jacobian matrix, the simulator (100) determines most significant harmonics (500) by adding least significant harmonics to an intermediate sum S and comparing the sum S to a predetermined first threshold. The equation system is then solved for the most significant harmonics. If the calculated results are in a predefined relation to a second threshold, the calculation is repeated. The simulator (100) converts intermediate results to the time domain and provides steady state responses of the circuit.

Description

METHOD FOR SIMULATING CIRCUITS USING HARMONIC BALANCE, AND APPARATUS THEREFORE

Field of the Invention

The present invention generally relates to data processing methods, and, more particularly, to a circuit simulating method and to an apparatus therefore.

Background of the Invention

Computer simulation is now an almost essential step in the design process of electronic circuits. A computer should estimate the circuit behavior within appropriate time and accuracy. However, the complexity of the circuit, the number of its elements, "linearity" or "nonlinearity" and other properties affect the computation costs (e.g., processing time, memory space) and to the working time of the design engineer. Conveniently, the computer simulates the circuit in a time domain in which physical quantities (e.g., voltages, currents) are functions of time or in a frequency domain. A steady-state response of the circuit is a plurality of time domain responses when the circuits receives one or more excitation signals with one or more frequencies. Simulation usually involves solving equation systems with a large number of equations and variables. Equation systems are conveniently solved in iteration steps. Such iterations usually begin with an assumed solution and end upon compliance of a criterion. An example for iteration is the Newton rule.

Using a harmonic balance (HB) method, the computer takes advantage of the well- known Fourier transformation. The computer performs HB nodal analysis according the Kirchhoff's current law in a suitably wide frequency range for substantially every node. For the application of harmonic balance, the following references are useful: [1 ] Kenneth S. Kundert and Alberto Sangiovanni-Vincentelli: "Simulation of Nonlinear Circuits in the Frequency Domain", IEEE Transactions on Computer-Aided Design, Vol. CAD-5, No.-4, October 1986;

[2] Jack T. Yao and Andrew T. Yang: "A Consistent Nonlinear Simulation Environment Based on Improved Harmonic Balance Techniques", IEEE, European Design Automation Conference 1993 and Euro-VHDL 1993, January 1993, pages 90-95; and [3] Kenneth S. Kundert, Jacob K. White and Alberto Sangiovanni-Vincentelli: "Steady- State Methods for Simulating Analog and Microwave Circuits", Kluwer Academic Publishers, Boston/Dordrecht/London, 1990, ISBN 0-7923--9069-5, especially chapter 5 "Harmonic Balance Theory" on pages 81-1 16 and chapter 6 "Implementing Harmonic

Balance" on pages 1 17-156.

The HB method is especially useful for circuits which have a relatively high number of non-linear elements. However, the number of equations to be solved and therefore processing time rises exponentially with the frequency range and the number of circuit elements. Accordingly, the use of the HB method is often limited in the prior art to circuits with a relatively low number of elements. In other words, large equation dimensions are the key problem for HB analysis.

In a prior art approach to this problem, some portions of the circuit are simulated in the time domain and other portions in the frequency domain using, for example, HB.

However, the human operator still has to decide which portions are suitable for time domain simulation (e.g., "linear" portions) and which portions are suitable for frequency domain simulation (e.g., "nonlinear" portions). This is usually not desirable because of the lost time and possibility for error. The present invention seeks to provide a circuit simulator and a method which mitigate or avoid these and other disadvantages and limitations of the prior art.

Brief Description of the Drawings

FIG. 1 illustrates a simplified magnitude-to-time diagram and a simplified amplitude- to-frequency diagram of a continuous signal; FIG. 2 illustrates a simplified block diagram of a memory portion for storing amplitude sets; FIG. 3 illustrates a simplified block diagram of a simulator of the present invention: FIG. 4 illustrates a simplified flow chart diagram of a method of the present invention: FIG. 5 illustrates a simplified flow chart diagram for the determining step of the method of FIG. 4; FIG. 6 illustrates a simplified frequency combination diagram: and FIG. 7 illustrates a simplified flow chart diagram for a modified determining step which is optionally used in the method of FIG. 4.

Detailed Description of a Preferred Embodiment

A prior art simulator generally performs the harmonic balance method with an equal number of frequencies (harmonics) for every node. Thereby, harmonics at some node which do not significantly contribute to the circuit behavior are in any case considered. In the steady-state response for a given plurality of circuit nodes, some nodes exhibit responses having many harmonics while other nodes exhibit responses with less harmonics.

According to apparatus and method of the present invention, the number of harmonics which need to be considered is automatically determined for every node before solving an equation system with a reduced number of equations. This feature substantially saves calculation time and can accelerate the circuit design process. Processing resources saved in this way can be used for additional iterations so that the overall simulation accuracy can be improved. Further, the human operator does not have to specify how many harmonics should be considered for each node. Also, the simulator of the present invention can simulate the behavior of circuits which receive two or more input signals with different frequencies, such as, for example, mixers.

Before explaining these and other features of the present inventions, some writing conventions and notation used herein are explained. Further, a glossary of terms used here and their definitions is provided prior to the claims. The notation ": = ", well known in the processing art, will be used to indicate that a variable on the left side obtains the value which is given on the right side. For example, K' : = K' - 1 means, that K' is decremented by 1. For convenience, the following description distinguishes between a single tone approach for a single input signal and a multi-tone approach for multiple input signals. The term "circuit" refers to any arrangement of active components (e.g., transistors) and passive component (e.g., components) which are related by connections (e.g., nodes, meshes). Circuits can be, for example, operational amplifiers or mixers used in radio applications. For convenience of explanation, the circuit is assumed to have a total number of I nodes.

FIG. 1 illustrates a simplified magnitude versus time diagram (time domain, a(i) = function(t) ) and a simplified amplitude versus frequency diagram (frequency domain, A(i,k) = function(f(k) ) of a continuous signal. As used herein, a "signal" is intended to include physical quantities, such as voltages, currents, charges, etc. and relations between such quantities, which are, for example transfer functions of transistors. In further explanations, the letters "a" and "A" in FIG. 1 can be replaced by other letters or acronyms (e.g., X or Res). The index "i" which distinguishes the signal from other signals is used in further explanations and is therefore introduced for convenience. For the discussion of FIG. 1 , index "i" is, however, not important.

The magnitude versus time diagram shown in the upper portion of FIG. 1 illustrates the magnitude "a(i)" of the signal on a vertical axis in reference to the time t on the horizontal axis. It is assumed that the signal is substantially periodic with a period Timeiϊ). The signal is illustrated as a time-continuous signal. This is convenient, but not essential for the present invention. The signal could also be a time- discrete signal. By applying any of the well-known transformations (e.g., Fourier- transformation, Laplace-transformation, Z-transformation), the signal a(i) can also be represented by the diagram shown in the lower portion of FIG. 1 :

The amplitude versus frequency diagram in the lower portion of FIG. 1 illustrates amplitudes A(i,k) on the vertical axis in reference to discrete frequencies f(k) on the horizontal axis. Such a diagram is also know as "spectrum". Frequencies f(k) (also "harmonics") are multiples of a base frequency f(l), that is: f(k) = f( l) * k (k = 0 to K), ( 1 .1 ) where K is the highest number of harmonics for which the transformation is performed. The frequency f(0) indicates the DC-part. For simplicity of explanation, the DC-part is not discussed here. FIG. 1 uses a single tone approach in which frequencies f(k) are given one-dimensionally. Combinations of frequencies which occur, for example in mixer circuits, are not considered here for simplicity. A two-dimensional spectrum diagram (multi-tone approach) is illustrated in connection with FIG. 6. A plurality of signal amplitudes A(i,k) belonging to all frequencies f ( 1 ) to f(K) is referred to as complete amplitude set { A(i) } , that is:

{ A(i) } = { A(i, l ), A(i,2). . . . , A(i.k) . . . A(i.K) } , (1 .2) with the { } symbols for sets.

According to the present invention, the number of harmonics which is taken into account can be reduced. A plurality of amplitudes A(i,k) belonging to the frequencies f ( 1 ) to f(K'(i) ) is referred to as partial amplitude set {A(i) } ', that is:

{ A(i) } ' = { A(i, l ), A(i,2), . . . , A(i,k) . . . A(i.K'(i) ) } , ( 1 .3) where K'(i) is the reduced number of harmonics for which calculations are actually performed. For convenience, K'(i) is further referred to as "partial set size". Preferably, K'(i) depends on i, so that for different signals (i), the numbers K'(i) can be different. The remaining amplitudes form a truncation amplitude set { A(i) }^τ, that is:

{ A(i) }^τ = { A(i,K'(i)+l ), . . . , A(i,K) } ( 1 .4)

The numbers and the sets are related as: K'(i) < K ( 1 .5)

{ A(i) } ' c { A(i) } ( 1 .6)

{ A(i) } = { A(i) } ' u { A(i) }^τ (1 .7)

In other words, the partial set is a subset of the complete set or both sets are identical

( c ); and partial set and truncation set in combination ( ) form the complete set. The signal representation in the frequency domain (e.g., sets { A(i) } and

{ A(i) } ' ) can be changed to the signal representation in the time domain (a(i) = function(t) ) by inverse transformation without substantially losing information. A feature of the present invention is that the number K'(i) can be automatically determined such, that errors caused by the difference between partial set { A(i) } ' and complete set { A(i) } can be substantially neglected. Taking into account the amplitudes A(i.k) belonging to an amplitude set, a norm of κ,-K- harmonics Norm(κ_v κ₂, i) can be defined as follows: k = κ₂

Nσr (κ, , κ₂,i) = ( Σ | A(i,k) | ² ) (1.8) k = K , The superscript 1/2 stands for "square root", and the | I symbols stand for absolute values. Norm( \ , K, i) is the norm of a complete amplitude set { A(i) } considering all K amplitudes A(i,k) from k = 1 to k = K. Norm(\ , K', i) is the norm of a partial amplitude set { A(i) } for K' amplitudes from k = 1 to k = K'. Residual norm(K' + l), K, i) is the norm of the truncation amplitude set { A(i) }^τ. A common residual norm can be defined for all i (from i = 1 to I), for example, by averaging.

For a prediction of amplitude values, it can be assumed that the amplitudes of the truncation amplitude set { A(i) }^τ approximately follow a geometric series, that is: A(i,k+1 ) = A(i,k) * for k > K' ( 1.9)

A(i,k) = A(i.K') * d^k'κ' for k > K' ( 1 .10) with d being a geometric progression factor. Factor d is, conveniently, smaller than 1 , that is: d < 1 (1.1 1 )

Persons of skill in the art will understand based on the description herein and well known general principles, how to approximate factor d from amplitudes A(i,k) belonging to the partial amplitude set (i.e., k < K'(i) ). There are many ways suitable to do this, such as, for example, regression, or least-square analysis.

The prediction of amplitudes is useful to estimate an allowable norm which is optionally used as a criterion in a comparing step (see FIG. 4, step 450). For example, equation ( 1.8) can be applied to predict amplitudes, so that an allowable norm Anormjpredicled (e.g., valid for all i) can be derived.

FIG. 2 illustrates a simplified block diagram of memory portion 1 10 for storing amplitude sets (e.g., { A(i) } ). FIG. 2 is intended to be a simplified example for illustration. Persons of skill in the art, are able, based on the description herein, to organize memory 1 10 is a different way without departing from the scope of the invention. Memory portion 1 10 comprises I memory fields 1 12-i (for each circuit node, i = 1 to I ). Fields 1 12- 1 to 1 12-1 can store a vector A of complete amplitude sets, that is:

A = [ { A( l ) } , . . ., { A(i) } , . . .. { A(I) } ] (2.1 )

The underscoring of "A" indicates "vector". As used herein, a vector represents the complete circuit. The vector elements are sets { } which represent the circuit nodes. Preferably, a single set { A(i) } is stored in field 1 12-i. Persons of skill in the art know how to divide memory field 1 12-i so that all amplitudes A(i,k) of set { A(i) } can be stored and accessed independently. Memory portion 1 10 also comprises I memory fields 1 14-i (i = 1 to I ) for storing partial set sizes K'(i) (cf. FIG. 1 ). The information in K'(i) indicates which amplitudes A(i,k) are used (e.g., from k = 1 to K') for calculation and which are not used (e.g., for k > K'). This is different from prior art approaches. Optionally, memory portion 1 10 can store only the partial amplitude sets ( A(i) } ' or the truncation amplitude sets { A(i) }^τ. Similarly, these sets form vectors:

A' = [ { A( l ) } ' , . . ., { A(i) } ' , . . ., { A(I) } ' ] (2.2) Δ^τ = [ { A( l ) }^τ , . . ., { A(i) }^τ , . . ., { A(I) }^τ ] (2.3)

FIG. 3 illustrates a simplified block diagram of simulator 100 of the present invention. Simulator 100 comprises memory 1 15 having a plurality of memory portions 1 10 and 1 1 1, processor 120 and interface 130. Processor 120 is coupled to memory portions 1 10 and 1 1 1 and to interface 130. While memory 1 15, processor 120, and interface 130 are shown as being separately connected, persons of skill in the art will understand that they can also be coupled by a common bus or separated busses.

Processor 120 performs calculations, exchanges data with memory portions 1 10 and 1 1 1 and with interface 130. Interface 130 is intended to comprise the elements which are needed for a communication between simulator 100 and the user. For example, interface 130 can comprise input devices (e.g., a keyboard), data loading and storage devices (e.g., a disk drive), and output devices (e.g., a display, or a printer) which are well known in the art.

Memory portion 1 10 has been explained in connection with FIG. 2. Memory portion 1 1 1 stores further variables conveniently used to implement the present invention, such as, for example, an auxiliary sum S; a power variable η, a predetermined threshold value v, iteration counters r and R, the number K or index sets

(e.g., { K,,K₂ } ), and other variables explained later. Also, memory portion 1 1 1 can optionally store control instructions for processor 120.

In other words, simulator 100 of the present invention can be described as an apparatus for simulating the behavior of a circuit which comprises: (a) a memory (e.g., memory 1 10) which stores for each circuit node a signal representation (e.g., set { A(i) } having components (e.g., amplitudes) at different frequencies (e.g., from k = 1 to K); (b) a logic (e.g., processor 120) coupled to the memory which, in repetitions (e.g., see FIGS. 5 and 7), assigns (e.g., adds) components at certain frequencies (e.g., for k > K'(i) ) where the components have minimum values, to an intermediate sum (e.g., sum S) until this intermediate sum reaches a predetermined threshold; and (c) a calculator (e.g., also implemented by processor 120) coupled to the memory which performs nodal analysis (e.g., by the harmonic balance method) only with unassigned components (i.e. for frequencies k < K' (i) ).

FIG. 4 illustrates a simplified flow chart diagram of method 400 according to a preferred embodiment of the present invention. Preferably, method 400 is applied using simulator 100 (FIG. 3). But, this is not essential, and references to simulator 100 serve only as an convenient example for purposes of explanation. A single tone approach is assumed here for convenience of explanation. A modification of method 400 for the multi-tone approach is explained in connection with FIGS. 6-7. In FIG. 4., blocks illustrate method steps, and lines with arrows illustrate a preferred method flow. When two blocks are illustrated in combination (e.g., 412 and 414), then either step can be performed first.

Method 400 comprises storing step 410 and solving step 420 (dashed frames). Storing step 410 comprises receiving node list step 412, providing initial vector X step 414, and formulating equation system step 416. Preferably, steps 412, 414 and 416 are performed consecutively. Bold lines "INPUT" to steps 412 and 414 indicate that simulator 100 receives data. Solving step 420 comprises the following substeps: determining step 500, solving subset of equation system step 430, calculating residual norm step 440, comparing step 450, providing step 46 0, transforming step 470, and checking step 490. Preferably, steps 500, 430, 440, 450. 460, 470, and 490 are performed consecutively. Depending on comparing step 450, steps 500, 430, 440 are repeated (cf. line 451); and depending on checking step 490, steps 500, 430, 440, 450, 460 and 470 are repeated (cf. line 492). An example for determining step 500 is explained with more detail in connection with FIG. 5. FIG. 7 illustrates step 501 as a modification of step 500. A bold line "OUTPUT" from step 470 indicate that simulator 100 provides output data, e.g., by magnitude-time diagrams for each circuit node signal (e.g., as in FIG. 1 ).

Method 400 can be described as a method for analyzing the steady-state response of a circuit in which a plurality of interconnections (e.g., nodes or meshes) are represented by a variable vector X and a residual vector R having multiple elements (e.g., { X(i) } ) in the frequency domain (e.g., with K harmonics for each element), the method comprising the steps of: (a) storing a circuit description (step 410) by a node list in a computer memory (e.g., memory 1 15 of simulator 100) (cf. step 412), providing initial values for X^(r) (r=l) for X (cf. step 414), and formulating equation system (step 416) with a Jacobian matrix; (b) solving the equation system by consecutive iteration steps (identified by superscript (r) ) using, for example the Newton rule, (step 420, see repetition line 491), wherein during each iteration, the following substeps are performed, for example, by processor 120:

(i) determining a sufficient number K' of harmonics (step 500, see also FIG. 5), preferably, for each vector element;

(ii) solving a subset of the equation system for each K' and providing intermediate results (step 430); (iii) calculating a residual norm related to K' (step 440);

(iv) comparing the residual norm to a predetermined threshold value (step 450) and conditionally repeating (cf. line 451) determining, solving, and calculating steps

(cf. line 452 to steps 500, 430, 440) or transforming the intermediate results to the time domain as part of the steady-state response (cf. line 452 to step 470); (v) providing X(r+1) for the next iteration (step 460), checking a completion criterion (step 490) and conditionally repeating the determining, solving and calculating steps with X(r+1) (cf. line 491 to steps 500, 430, 440).

In receiving node list step 412, memory 1 15 of simulator 100 receives a representation of the circuit topology via lists or tables through interface 130. There are many well known ways how to store circuit representations. For example, the representations can be obtained by parsing the representations into a netlist. Persons of skill in the art will know how to carry out step 412.

In providing initial vector X^(1> step 414, processor 120 determines start values for X^(r) (r = 1 ) for further iterative calculations. X^{" is a variable vector representing any physical values (e.g., voltages, currents) due to any given interconnection of circuit components. Superscript (r) is an iteration counter r (e.g., r = 1 to R). Persons of skill in the art know how to provide initial values ( r = 1), for example, when further calculations should be performed using the Newton rule. Step 414 can optionally comprise reading predetermined threshold value v which is used in determining step

500.

In formulating equation system step 416, processor 120 reads the circuit representations (e.g., from memory 1 15) and establishes an equation system. Persons of skill in the art will know how to perform step 416. For example, taking a convenient illustration from chapter 5 of reference [3], the equation system can be expressed by a residual:

Res(V) = CurrentQD + ^Ω * Charge(V) + Y * V (4.1 ) with symbol * for multiplication. The equation system with the vector still represents all K frequencies for all I circuit nodes. Persons of skill in the art can use other equation systems which are suitable to describe the circuit for all frequencies. The system can be modified by adding terms or deleting terms. In the example, symbol V being an example for X stands for an I-element vector of the I voltages V,, V.,, V,, ... , V, of the I nodes of a circuit, all in the frequency domain. Each X can be considered as a complete amplitude set { X(i) } (cf. equation ( 1.2) for all harmonics K ). In the example of equation system (4.1), the variable "Current" is the Fourier transform the current as function of time t; and the variable "Charge" is the

Fourier transform of the charge as function of time t. Symbols Ω for a frequency matrix and Y for an admittance matrix of the linear circuit elements are also explained in chapter 5 of reference [3]. To obtain a steady-state response, system (4.1) should be solved with (writing X for V)

Res(X) = 0 (4.2) with respect to Fourier coefficients (harmonics) of voltage vector N In the following, the equation system is generally be expressed as Res(X) = 0. Preferably, processor 120 writes equation system (4.1 ) in the Jacobian form, that is: J(X^(r)) [AX ⁾] = - Res(X^(r) ) (4.3) wherein "J" stands for the Jacobian Matrix, and [ΔX^(r)] means that the left side of equation (4.3) is applied to X^(r) and X^(r+n of two consecutive iterations. Initially,

ΔX^(r) = X( l) (4.4)

In determining step 500, the numbers of harmonics K'(i) which are to be considered in solving the subset of the equation system in step 430 are determined. Step 500 is an important feature of the present invention. Step 500 is, preferably, applied to all i = 1 to I elements of vector Res(X^(r) ) on the right hand side ("rhs") of equation (4.3). But this is not essential, so that the general descriptor { A(i) } is used for detailed explanation of system (4.3). In step 500, the elements of the rhs with small contribution are neglected. More details are explained in connection with FIG. 5. Reducing the number of harmonics being considered is a big advantage of the present invention over the prior art. The computation time of processor 120 can be substantially saved. Also, determining step 500 does not require any interaction between simulator 100 and a human operator. In further repetitions of step 500, for example, during further iteration repetitions (r), the numbers K'(i) can be increased to obtain higher accuracy, if that is desired.

As it can be seen from FIG. 4, determining step 500 is executed within an inner loop (cf. line 451 ) and an outer loop (cf. line 491). The number of harmonics K'(i) may be varied for different iterations and for different parts of the circuit (i.e. K' depending on i).

In solving subset of equation system step 430, the equation system (4.3) is solved wherein for each node i, only K'(i) < K harmonics are considered. In other words, processor 120 solves a subset of (4.3): J'(X^(r)) [ΔX^(r)] = Res'(X^(r) ) (4.5)

The prime marker (') attached to the J-symbol and to the Res-symbol follows the writing convention of equation (1.3) in which partial amplitude sets have been introduced. The vector elements (e.g., i = 1 to I) of X^w are partial amplitude sets { X(i) } ' (with the prime marker). In other words, equation system (4.5) is a reduced equation system which does not consider all K amplitudes, but instead considers only selected amplitudes (i.e., different K'(i) for each i). Equation system (4.5) comprises fewer equations than the original system (4.3). This is a big advantage of the present invention. For each X(i), all K harmonics should not be considered.

Preferably, processor 120 solves equation system (4.5) using a Krylov subspace iterative technique. For example, processor 120 temporarily determines Krylov subspaces and obtains solutions by iterating. This is convenient, but not essential for the present invention. A useful reference is: [4] Yousef Saad: "Iterative Methods for Sparse Linear Systems", PWS Publishing Company, 1996, ISBN 0-534-94776-X. In reference [4], Krylov approaches have been described in chapter 6 on pages 144-204 and in chapter 7 on pages 205-229. From the methods explained in [4], the

Generalized Minimum Residual Method (GMRES) described in section 6.5 of chapter 6 (pages 158-173) is especially useful applied in step 430. Such a technique minimizes the residual norm. Those of skill in the art can implement step 430 without the need for further explanations. Conveniently, preconditioning can be applied in step 430. Preconditioning is well known in the art and described, for example in section 4.1.2. of chapter 4 ("Iteration Matrices and Preconditioning", pages 102- 103).

In calculating residual norm step 440, processor 120 calculates residual norms Nor ((K'(i)+ l ), K, i) according to equation ( 1.8) for each partial amplitude set of X^(r). Optionally, processor 120 calculates the residual norm as the common residual norm (e.g., average). In comparing step 450, processor 120 compares the residual norm to an allowable norm ANorm. Comparing step 450 is a query step which controls method 400 according to the question ("Norm allowable ?"). If the residual norm (of previous step 440) stays below the allowable norm ("yes"-line 452), then processor 120 continues with step 460. Otherwise, if the residual norm exceeds the allowable norm ("no' ine 451 ), then processor 120 repeats steps 500, 430 and 440. Comparing step 450 is an optional step in method 400 which prevents an accumulation of errors.

For example, processor 120 compares the residual norm to a predetermined allowable norm ANorm_predetermined. Or, the residual norm is compared to an allowable norm ANorm_predicted which is predicted, for example, from a geometric series of amplitudes (cf. equations (1.9) to (1.1 1) ). Other criteria for comparing step 450 can also applied.

In providing step 460, processor 120 provides X^(r+1) for the next iteration. Persons of skill in the art know how to obtain this values. The total number of iterations is conveniently limited by a number r = R.

In transforming step 470, processors 120 performs a transformation of some or all elements of vector X(r) from the frequency domain into the time domain by inverse transformation. For example, processor 120 applies inverse Fourier transformation and provide voltage values for some or all nodes of the circuit. Conveniently, processor 120 sends time domain results as output data to interface 130 or to memory 1 15 for temporary storage.

In checking step 490, processor 120 applies a convergence criterion for the r^lh iterations. Such criteria are well known in the art. Depending on the results of step 490, method 400 is finished (end line 492) or the next repetition of steps 500, 430, 440, 460, 470 and 490 is initiated (cf. line 491 ).

FIG. 5 illustrates a simplified flow chart diagram for determining step 500 of method 400. Line 41 1 at the beginning of step 500 and line 560 at the end of step 500 connect FIG. 5 to FIG. 4. Unless otherwise specified, the explanation of step 500 refers to the current iteration r. Preferably, determining step 500 has the following substeps: initializing step 505, scanning step 510, accumulating step 520, setting step 530, comparing step 540 and conditionally repeating (line 550). Determining step 500 is explained here for a single tone approach.. Optional modifications of determining step 500 for a multi-tone approach are explained in connection with FIG. 6. In initializing step 505, processor 120 sets all K'(i) to a predetermined value of, preferably:

K'(i) : = K (i = 1 to I) (5.1 ) Also, processor 120 initializes an auxiliary sum S to S₀ (preferably S₀ = 0), that is : S : = S₀ (5.2)

In scanning step 510, processor 120 consecutively reads amplitudes A(i,K'(i) ) of { A(i) } ' from memory portion 1 10 (initially { A(i) } ' = { A(i) } ) and determines an index L ( 1 < L < I ), for which A_mιn = A(L,K'(L) ) has a minimum value, that is:

A(L,K' (L) ) < A(i,K'(i) ) V l < i < I (5.3) with V standing as "for all"). In other words, processor 120 identifies a single set

{ A(l) } of vector A whose last amplitude A(K'(L) ) does at least contribute to vector A. In scanning step 510, processor 120, preferably, counts up i (start with i =1 , ends with i = I). This is convenient, but not essential for the present invention, so that processor 120 can use any other order of i.

In accumulating step 520, processor 120 adds A(L, K'(L) ) to the power of η

(superscript) to the sum S, that is: S : = S + [ A(L, K'(L) ]^η (5.4)

Conveniently, accumulating step 520 is performed by adding the squares of A(K'(L) ). that is: η = 2 (5.5)

This is convenient, but not essential for the present inventions, so that other values for η can also be used. In other words, in step 520, processor 120 shifts the least important amplitudes A(i,K'(L) ) from vector A to the intermediate sum S.

In setting step 530, processor 120 decrements K'(L) by 1 and stores the new value in memory portion 1 10, that is:

K'(L) : = K'(L) - 1 (5.6) In comparing step 540, processor 120 compares sum S to the threshold value v, that is, for example:

Is S ≥ v ? (5.7)

Optionally, the threshold value v can depend as v(r) on the current iteration r performed in method 400. This is convenient, but not essential. Depending on a comparison result in step 540, processor 120 continues determining step 500 conditionally. For a first result (e.g., S > v, "yes"-line 550 to step 5 10), processor

120 repeats steps 510, 520, and 530. For a second result (e.g., S < v, "no"-line 560), processor 120 ends determining step 500 and continues method 400 with step 430 (see FIG. 4). Persons of skill in the art are able, based on the description herein to have processor 120 performed comparison step 540 earlier, for example, prior to step 510. In other words, processor 120 repeats steps 510, 520, 530, for example, as long as S stays below the predetermined value v. Or, with further words, processor 120 conditionally repeats steps 510, 520, 530 and 540 until the amplitudes which have been shifted from vector A to sum S have reached or exceeded a threshold value (e.g., v).

At the end of determining step 500, the vector A of complete amplitude sets has been virtually reduced to vector A' of partial sets { A(i) } ' and numbers K'(i) have been modified to a pointer vector K', that is: K' = { K'( 1 ), K'(2), . . . , K'(I) } (5.8)

Preferably, all amplitudes A(i,k) of vector A (all i, all k) remain completely in memory portion 1 10. This is an important feature of the present invention. In further iterations r+1 , the vector A' can be derived from vector A with a different pointer vector K' . FIG. 6 illustrates simplified frequency combination diagram 600 for an amplitude set { A(i,k,,k₂) } . Simulator 100 of the present invention can also be used to simulate circuit 100 which receives two or more frequencies (e.g., mixer). In this multi-tone approach, the signal amplitudes A(i,k) depend on two or more discrete frequencies. Hence, the argument f(k) of A(i.k) = function (f(k) ) (see FIG. 1 ) is generalized to p > 1 frequencies f_p(k_p), that is: A(i ,k , , . . ., k_p) = function ( f, (k ,), . . ., f_p(k_p) ) (6.1 ) by conveniently using additional subscripts p for partial frequencies f_p(k_p) and for indices k . For convenience of explanation, the following description assumes p = 2 for first frequencies f,(k,), second frequencies f (k₂) and amplitudes A(i,k,,k₂). Diagram 600 has a horizontal axis for frequencies f,(k,) and a vertical axis for frequencies f₂(k₂). In a grid defined by the k, and k₂ indices, FIG. 6 illustrates amplitudes A(i,k,.k₂) by • and O symbols. In this two-dimensional representation of

FIG. 6, the amplitude magnitudes are not illustrated. Also, and intended only for the convenience of explanation, the partial frequencies are illustrated as being non-negative. For example, frequencies f,(k,) have magnitudes from, e.g., f,( l ) = 10 kHz to, e.g.. f,(5) = 10*5 = 50 kHz; and frequencies f₂(k₂) have magnitudes from e.g., f₂( l ) = 15 kHz to, e.g., f₂(6) = 6* 15 kHz = 90 kHz. Frequencies f,(0) and f₂(0) stand for the DC-amplitudes. Similar to the definition of K as the highest number of harmonics (see FIG. 1), a complete amplitude set { A(i) } of combined harmonics can be defined by a set { K,, K₂ } (i) of index pairs (k,,k₂) that is, for example,

{ K, ,K₂ } (i) = { (0,0), (0,1 ), . . . , (0,K₂), (6.2) ( 1 ,0), ( 1 , 1 ), . . . , ( 1 ,K₂),

(K, ,0), (K„0), . . . , (K„K₂) }

Amplitudes A(i,k,,k₂) (• and O symbols) of the complete set { A(i) } (for

{ K,.K₂ } ) are enclosed by dashed frame 610. The rectangular shape of frame 610 is convenient for explanation. The complete set { A(i) } can also be defined by other index combinations (k,,k₂), which would result in, for example, diamonds, circles, or other otherwise shaped enclosures.

Similarly as explained in connection with FIG. 1 , a partial amplitude set { A(i) } ' of some amplitudes A(i,k,,k₂) (• symbols within dashed frame 620) can be defined by a subset { K,,K₂ } '(i) c { K,,K₂ }(i), that is, for example: (6.3)

{ K„K₂ } '(i) = { (2,2), (2,3), (2,4), (3,2), (3,3), (3,4), (3,5), (4,3), (4,4) } ln{ A(i) } ', the index combinations are, preferably, neighboring frequencies. For the following explanation, it is convenient to define the terms "border amplitude" and "corner amplitude". A border amplitude is an amplitude having such an index combination (k,,k₂) of the set { K,, K₂ } which would leave a given amplitude set { A(i) } when one of the indices (k, or k₂) is changed and the other index (k₂ or k,) stays constant. For example, border amplitude 614 (indices (5,4) ) would go out of set { A(i) } (frame 610) when k, would be changed by 1 (from 5 to 6) and k₂ would stay constant. Similarly, border amplitude 624 (indices (2,3) ) would leave set { A(i) } ' (frame 620) when k, becomes k,=l and k₂ would stay constant.

A corner amplitude is defined for such index combinations where changes of either index (k, or k₂) can exit the amplitude set (e.g., { A(i) } ). For example, corner amplitude 612 would leave { A(i) } (frame 610); and corner amplitude 622 would leave { A(i) } ' (frame 620). The corner amplitudes are similarly defined as the amplitudes A(i,k) for k = K and k = K' in the one-dimensional representation of FIG. 1. Hence, they can be treated accordingly as explained in the following.

FIG. 7 illustrates a simplified flow chart diagram for modified determining step 501 which is optionally used in method 400. For the multi-tone approach, determining step 500 is modified to determining step 501. Steps 505, 510, 520, and 530 are modified to steps 506, 51 1 , 521 and 531 , respectively. Step 540 with repetition 550 is substantially unchanged. For convenience of explanation, only the modifications in step 501 compared to step 500 are explained here.

In initializing step 506, processor 120 defines the complete amplitude set { A(i) } of combined harmonics, for example sets { K K₂ }(i) (frame 610). Preferably, the set are equal for all i. In scanning step 51 1 , processor 120 first identifies the corner amplitude A(i,k,.k₂)_mιn with has the minimum magnitude of the corner amplitudes for current index i. For example, processor 120 compares corner amplitudes A(i,0,0), A(i,0,6), A(i,5,0), and A(i,5,6) and identifies corner amplitude A(i,5,6) with reference 612 as minimum amplitude A(i,k,,k₂)_m]n. Second, processor 120 compares this minimum amplitude A(i,k,,k₂)_min to the other minimum amplitudes A(i,k,,k₂)_mm for all i = 1 to I and determines index 1 accordingly. In accumulating step 521 , the minimum amplitude (e.g., A(l,5,6) ) is added to sum S. In setting step 531 , a new index set { K,,K₂ } '(1) is defined (for example, frame 610 without corner 612). Thereby, new corner amplitudes are obtained (e.g., (4,6) and (5,5) ). After comparing step 540, steps 51 1 , 521 and 531 are conditionally repeated ("yes"-line 550). Thereby, processor 120 identifies the amplitudes with the O symbols and accumulates them to S. When processor 120 goes to the "no"-line 560 after comparing step 540, the index set has reached its final shape (e.g., frame 620) and the partial amplitude set { A(i) } ' (• symbols) is established. Method 400 continues with step 430. Having described details, the substeps of determining step 500 are now explained as a method. This method is conveniently applied consecutively or substantially simultaneously applied to each complete set { A(i) } , for example, in simulator 100 which solves, e.g., equation system (4.3) with amplitudes in a domain of

(k,, . . . k_p ) discrete frequencies and combinations thereof (single-tone and multi- tone approach). The method determines a final partial set { A(i) } ' from an initially complete set { A(i) } of amplitudes and has the following steps:

(a) scanning all amplitudes A(i,k,, . . . k ) of set { A(i) } to find that frequency combination for which a single amplitude e.g., A(1,K'(1) ) or A(i,k,,k₂)_mιn has a minimum magnitude (e.g., steps 510, 51 1 ); (b) accumulating the η^,h power of the single amplitude ("min") to sum S (initially a predetermined value, of e.g., S=0, cf. step 505), for example in steps 520 or 521 ;

(c) modifying the set { A(i) } to a first partial set { A(i) } " by taking the single minimum amplitude out of set { A(i) } (e.g., by setting step 530, 531 ); and

(d) comparing sum S to predetermined threshold value v and repeating scanning step (a) for first and for further partial sets { A(i) } " until S and v are in a predetermined relation (e.g., larger or equal than relation, step 540) so that the final partial set { A(i) } ' is established.

This method can further comprise: (i) solving a subset (e.g., 4.5) ) of equation system (e.g., (4.3) ) using partial set { A(i) } ' and providing results (e.g., as in step 430) and (ii) conditionally repeating scanning, accumulating and modifying steps (e.g., steps 510/51 1 , 520/521 , 530/531) if the residual norm is too high.

The "partial harmonic balance" method of the present invention can be used to obtain equal or higher accuracy within the same or less total processing time compared to a "complete harmonic balance" of the prior art. More iterations (r) can be performed, or the calculation time can be shortened significantly or a combination thereof.

While the invention has been described in terms of particular structures, devices and methods, those of skill in the art will understand based on the description herein that it is not limited merely to such examples and that the full scope of the invention is properly determined by the claims that follow.

Glossary of Terms

{ } indicates sets underscoring indicates vectors ' stands for "partial"

* multiplication

V "for all" c subset of combination of sets a(i) signal a having index i in the time domain

A(i,k) amplitude of signal a(i) in frequency domain for k^Λ frequency

ANorm allowable norm d geometric progression factor

DC abbreviation for "direct current" E predetermined threshold f(k), f(k,,k₂) discrete frequency i, I index

J Jacobian matrix k index for harmonics K(i) highest number of harmonics

K'(i) reduced number of harmonics, "partial set size" L index min (subscript) minimum value

Norm norm r, R iteration counter

Res residual vector

S auxiliary sum

Time(ϊ) time period

T (superscript) stands for truncation v, V voltage

X(i) vector element in the frequency domain

X vector

Y cf. [3] κ₂ used as k-indices

η indicates power

V threshold value

P index for frequencies

Ω cf. [3]

Claims

1. A method for simulating a circuit having a plurality of nodes described by a vector X with a plurality of N variables X(i) (i=l to N) and wherein each variable X(i) has a predetermined number of harmonics (k=l to K(i) ), said method comprising the steps of:

(a) storing a description of the circuit as an equation system wherein X is a variable;

(b) providing initial values for X(i);

(c) selecting K'(i) < K(i) harmonics for each variable X(i); (d) solving a subset of said equation system for which (k=l to K') to produce an intermediate solution; and (e) repeating said steps (c) and (d) for different values until the intermediate solutions meet a predetermined convergence criterion.

2. The method of claim 1 wherein in said step (a), said equation system is: J ( )[╬öX] = - R (X), with J the Jacobian matrix and R a residual vector.

3. The method of claim 1 wherein in said solving step comprises temporarily determining Krylov subspaces and iterating.

4. The method of claim 1 wherein a single tone input signal determines said variables X(i).

5. The method of claim 1 wherein multi-tone input signal determine said variables X(i).

6. The method of claim 1 wherein said variables X(i) are node voltages.

7. The method of claim 1 wherein said K'(i) is equal for each i.

8. The method of claim 1 further comprising the following steps which are executed after step (d):

(dl) deriving a residual norm from the intermediate solution of said subset; and (d2) repeating said selecting and solving steps when said residual norm exceeds a predetermined norm.

. The method of claim 8 wherein in said predetermined norm is derived by geometric series prediction from amplitudes for harmonics k < K' relative to harmonics k > K ' .

10. A method for analyzing a steady-state response of a circuit in which a plurality of interconnections between circuit components is represented by a variable vector X and a residual vector Res, said vectors X and Res having multiple elements in the frequency domain each with K harmonics, said method comprising the steps of: storing a circuit description by a netlist in a memory of a computer, and providing initial values X^(r) (r=0) for X and representing an equation system:

J(X^(r)) [╬öX^(r)] = - Res(X^(r) ) in said memory, with J(X^(r)) [╬öX^(r)] being the

Jacobian applied initially to ╬öX^(r) =X^<r) (r=0) and ╬öX being a difference between iterations X^(r> and X^(r+I) of vector X; and solving said equation system by consecutive iteration steps identified by (r) wherein each iteration step has the following substeps:

ΓÇó determining a sufficient number K' of harmonics;

ΓÇó solving a subset of said equation system for said K' harmonics and providing intermediate results; ΓÇó calculating a residual norm related to K';

ΓÇó comparing said residual norm to a predetermined threshold value and, depending on comparison result, conditionally repeating said determining, solving and calculating steps

ΓÇó transforming said intermediate results to a time domain as part of said steady- state response; and

ΓÇó providing X(r+1 ), checking a completion criterion and conditionally repeating said determining, solving and calculating steps with X(r+1 ).

1 1. In a simulator which solves an equation system with amplitudes in a domain of

(k . . . k ) discrete frequencies and combinations thereof, a method for determining a final partial set { A(i) } ' from an initially complete set { A(i) } of amplitudes, said method having the following steps: scanning all amplitudes A(i,k . . . k ) of said set { A(i) } to find that frequency combination for which a single amplitude has a minimum magnitude A--;, finding a sum S = S₀ + A_mn ^╬╖ where S₀ is an initial value and ╬╖ is a predetermined exponent; modifying said set { A(i) } by deleting said single amplitude to obtain a first partial set { A(i) } " = { A(i) } - A_m╬╣n; and repeating said scanning step for { A(i) } " until said sum and a predetermined threshold value are in a predetermined relation so that said final partial set

{ A(i) } ' is established.

12. The method of claim 1 1 further comprising:

(i) solving a subset of said equation system using said partial set { A(i) } ' and providing results and (ii) conditionally repeating said steps scanning, accumulating and modifying if a residual norm exceeds an allowable value.

13. The method of claim 1 1 wherein in said scanning step amplitudes are scanned for frequencies f(k) from f(l) to f(K).

14. A method for simulating the steady-state response of non-linear circuits having a plurality of devices and a plurality of interconnecting nodes, said method comprising the steps of: modeling said circuit in a frequency domain; and applying harmonic balance calculation using a variable number of harmonics for said interconnecting nodes.

15. A method for simulating the steady-state response of a nonlinear circuit being represented by a plurality of devices and a plurality of interconnecting nodes, said method comprising the steps of: a) receiving an input netlist with a circuit description and storing said input netlist in a computer memory; b) providing initial values of circuit variables for further iteration; and c) solving an equation system J(X^(r)) ╬öX^(r) = - Res(X^(r) ) with Jacobian matrix J and residual Res in iterative steps (r), wherein in each of said steps: c l) a number of harmonics K' for each node is determined, c2) the solutions for said equation system are provided for K' < K harmonics; c3) new values X(r+1) for the next iteration cycle are obtained;

16. The method of claim 15 wherein in step c), Fourier transformation and inverse Fourier transformation is applied.

17. The method of claim 15 wherein in step c2) the solutions are provided by applying iterative Krylov subspace techniques.

18. An apparatus for simulating the behavior of a circuit, comprising: a memory which stores for each circuit node a representation having signal components at different frequencies; a processor coupled to said memory which assigns components at certain frequencies where the signal components which have minimum values to an intermediate sum until this intermediate sum reaches a predetermined threshold, thereby leaving unassigned components; and then performs nodal analysis on the circuit using the unassigned signal components.