Abstract:
Methods for analyzing circuit distortion based on contributions from separate circuit elements are presented. Local approximations that do not require high-order derivatives of device models are developed near an operating point for calculating distortion summaries including compression summaries and second-order intermodulation (IM2) distortion summaries.

Description:
CROSS-REFERENCE TO RELATED APPLICATIONS 
   This application claims the benefit of U.S. Provisional Application No. 60/750,706, filed Dec. 14, 2005, and incorporated herein by reference in its entirety. 

   BACKGROUND OF THE INVENTION 
   The present invention relates to the analysis of circuits generally and more particularly to the distortion analysis of analog and RF (Radio Frequency) circuits by calculating distortion contributions of individual circuit elements. 
   Circuit linearity is a critical specification in analog and RF designs. Intermodulation product and compression are two commonly adopted measures for distortion [6, 11]. As a result, a general methodology for analyzing effects of individual components on circuit distortion and for identifying the major sources of nonlinearity is of practical importance in complex designs. From the designer point of view, it is desirable to separate total distortion into contributions of different elements. While in general such a separation is not possible for most nonlinear effects because of cross terms between devices, insights on distortion dependency on nonlinearity in a single device can still be obtained by conducting some heuristic analysis. At low RF power levels, perturbation methods based on Volterra series can be applied to calculate nonlinear effects [10, 2, 5, 1, 9, 4, 14, 15]. This approach can provide relatively detailed information related to distortion caused by different nonlinear terms in different elements. 
   A drawback of the Volterra method is that it generally requires high order device derivatives, which usually are not provided by device models used in real design simulations. As a result, the application of Volterra series has typically been limited to simple models and simple circuits. To circumvent this difficulty, Li and Pileggi proposed a post-simulation post-processing approach based on linear-centric circuit models, where the nonlinear current of an element is computed by subtracting the current in a linearized device model from the total current calculated by full device model evaluation [8]. According to this approach, the terminal voltage used for device evaluations is taken from the steady state solution which includes nonlinear effects of the whole circuit. Therefore, the nonlinear current obtained contains cross terms implicitly and is not totally independent of other elements. 
   In general, current methods do not adequately model circuit distortion in terms of contributions from separate device elements. Thus, there is a need for improved methods and systems for circuit distortion analysis, especially with applications to analog and RF circuits. 
   SUMMARY OF THE INVENTION 
   In one embodiment of the present invention, a method of analyzing distortion in a circuit includes: determining an operating point for the circuit, where the circuit has a linear offset and a nonlinear offset for analyzing distortion at the operating point; determining first-order output values by using the linear offset to relate an input having an input frequency with the first-order output values; specifying one or more devices for analyzing the circuit; determining first nonlinear-offset values for the one or more devices at twice the input frequency and at a zero frequency by evaluating the nonlinear offset at values based on the first-order output values; determining second-order output values for the one or more devices at twice the input frequency and at the zero frequency by using the linear offset to relate the second-order output values with the first nonlinear-offset values; determining second nonlinear-offset values for the one or more devices at the input frequency by evaluating the nonlinear offset at values based on the first order-order output values and the second-order output values; determining third-order output values for the one or more devices at the input frequency by using the linear offset to relate the third-order output values with the second nonlinear-offset values; and saving at least some values based on the output values. 
   According to one aspect of this embodiment, determining the operating point for the circuit may include approximating the operating point by a time-domain method or a frequency-domain method so that the linear offset defines a linear function for inputs about the operating point and the nonlinear offset defines a nonlinear function for inputs about the operating point. 
   According to another aspect, the method may further include determining one or more compression summaries for the one or more devices by calculating one or more comparisons between third-order output values and corresponding first-order output values. Additionally, calculating the one or more comparisons may include calculating logarithmic ratios between third-order output values and corresponding first-order output values. 
   In another embodiment of the present invention, a method of analyzing distortion in a circuit includes: determining an operating point for the circuit, where the circuit has a linear offset and a nonlinear offset for analyzing distortion at the operating point; determining first-order output values by using the linear offset to relate an input having a first frequency and a second frequency with the first-order output values; specifying one or more devices for analyzing the circuit; determining nonlinear-offset values for the one or more devices at a difference between the first frequency and the second frequency by evaluating the nonlinear offset at values based on the first order-order output values; determining contributions to distortion for the one or more devices by projecting the nonlinear-offset values based on the difference between the first frequency and the second frequency; and saving at least some values based on the contributions to distortion for the one or more devices. 
   According to one aspect of this embodiment, determining the operating point for the circuit may include approximating the operating point by a time-domain method or a frequency-domain method so that the linear offset defines a linear function for inputs about the operating point and the nonlinear offset defines a nonlinear function for inputs about the operating point. 
   According to another aspect, the method may further include determining a second-order intermodulation (IM2) distortion summary for the one or more devices by combining the contributions to distortion for the one or more devices. Additionally, combining the contributions to distortion for the one or more devices may include summing the contributions to distortion for the one or more devices. 
   According to another aspect, the method may further include determining a transfer function that projects the nonlinear-offset values onto one or more output nodes of the system for calculating the contributions to distortion for the one or more devices. Additionally, projecting the nonlinear-offset values may include calculating one or more inner products between the transfer function and the nonlinear offset values. 
   Additional embodiments relate to an apparatus for carrying out any one of the above-described methods, where the apparatus may include a computer for executing instructions related to the method. For example, the computer may include a processor with memory for executing at least some of the instructions. Additionally or alternatively the computer may include a specialized microprocessor or other hardware for executing at least some of the instructions. Additional embodiments also relate to a computer-readable medium that stores (e.g., tangibly embodies) a computer program for carrying out any one of the above-described methods with a computer. 
   In this way the present invention enables improved methods and systems for circuit distortion analysis with applications to analog and RF circuits. 

   
     BRIEF DESCRIPTION OF THE DRAWINGS 
       FIG. 1  shows an exemplary system which is applicable to embodiments of the present invention. 
       FIG. 2  shows an exemplary method for computing compression distortion summaries. 
       FIG. 3  shows computed values related to the embodiment shown in  FIG. 2 . 
       FIG. 4  shows computed values related to the embodiment shown in  FIG. 2 . 
       FIG. 5  shows an exemplary method for computing IM2 distortion summaries. 
       FIG. 6  shows computed values related to the embodiment shown in  FIG. 5 . 
   

   DETAILED DESCRIPTION OF EXEMPLARY EMBODIMENTS 
   1. Introduction 
     FIG. 1  shows a system  102  (e.g., a circuit model) that receives an RF input  104  with an input frequency ω RF  and produces an output  106  where the input frequency has been modulated by the system (e.g., shifted by a multiple of the LO (Local Oscillator) frequency ω LO . Under non-ideal conditions, the output  106  includes distortions due to the nonlinear effects in the system  102 . As described below, the present invention enables quantitative characterizations of how component devices contribute to distortion in the output signal  106 . Specific embodiments include calculations for compression distortion and second-order intermodulation products (IM2). This application is related to CALCULATING INTERMODULATION PRODUCTS AND INTERCEPT POINTS FOR CIRCUIT DISTORTION ANALYSIS (application Ser. No. 11/303,049, filed Dec. 14, 2005), which is incorporated herein by reference in its entirety. 
   2. Compression Distortion Summary 
   The circuit  102  in  FIG. 1  has a periodic time-dependent operating point V 0 (t) at the LO frequency ω LO  More generally, systems with DC (direct current) operating points such as an LNA (Low Noise Amplifier) can be considered as special cases with ω LO =0. When the RF input signal  104  is applied, the output signal  106  at ω RF +nω LO  is a function of the input magnitude V RF  and given as:
 
 V   out (ω RF   +nω   LO )= c   1   V   RF   +c   3   V   RF   3 + . . .   (1)
 
   The first term on the right-hand side of eq. (1) is the linear response to RF input, and the coefficient c 1  can be obtained by small signal calculation. Other terms model distortion due to third and higher order terms in the nonlinear response. The compression distortion can be measured in dB (decibels) as 
   
     
       
         
           
             
               
                 dB 
                 = 
                 
                   10 
                   · 
                   
                     
                       log 
                       ⁡ 
                       
                         ( 
                         
                           
                             V 
                             out 
                           
                           
                             
                               c 
                               1 
                             
                             ⁢ 
                             
                               V 
                               RF 
                             
                           
                         
                         ) 
                       
                     
                     . 
                   
                 
               
             
             
               
                 ( 
                 2 
                 ) 
               
             
           
         
       
     
   
   While it&#39;s desirable to know the contribution of each individual device to total compression, it is mathematically intractable to make such a separation for third or higher order nonlinearities. In order to quantify the effect of a specific device on the nonlinear response at the output signal, we assume the rest of the circuit responds to the RF signal linearly and compute the distortion in the output caused solely by nonlinear response at this device. With this approach, although the circuit is hypothetical, the result provides useful information about how sensitive distortion is to the device. The effects of the device on the nonlinear response and the transfer function from device nodes to output nodes are reflected in the calculation. 
   Notably, since the operating point is time-dependent, all the other devices are still nonlinear in the calculations. Only their response to RF signal is linearized. In particular, for mixer cases, transistors in LO block still function nonlinearly to generate the clock signal. However, their contribution to the output distortion summary is expected to be very small because the RF signal  104  usually does not disturb the LO part of the circuit. Both linear and nonlinear responses in LO block are typically almost zero. 
   We formulate compression distortion summary calculations by writing circuit equation as 
                         ⅆ     ⅆ   t       ⁢     q   ⁡     (     υ   ⁡     (   t   )       )         +     i   ⁡     (     υ   ⁡     (   t   )       )         =       B   ⁡     (   t   )       +     s   ⁡     (   t   )           ,           (   3   )               
where B is the operating point source such as LO at ω LO  and s is RF input signal at ω RF  [13]. For simplicity, we rewrite eq. (3) as
   F (υ)= B+s.   (4) 
Assume V 0 (t) is the steady-state (or DC) solution at zero RF input so that V 0  satisfies equation
   F ( V   0 ( t ))= B ( t ).  (5) 
   Next, expanding eq. (4) around V 0  gives
 
 L·v+F   NL ( v )= s,   (6)
 
where v is the circuit response to RF signal s. The linear operator L in the first term is defined as
 
                   L   ·   v     =         ⅆ     ⅆ   t       ⁢     (       ∂     ∂   υ       ⁢       q   ⁡     (     V   0     )       ·   v       )       +       ∂     ∂   υ       ⁢       i   ⁡     (     V   0     )       ·   v                 (   7   )               
and L·v is just the linearized LHS of eq. (3). The second term F NL  represents the sum of nonlinear terms in all the devices and is defined as
   F   NL ( V   0   ,v )= F ( V   0   +v )− F ( V   0 )− L·v.   (8) 
That is, linearizing about the DC or operating point V 0  determines linear and nonlinear offsets (or offset functions), L and F NL , for analyzing local perturbations of the output.
 
   The linear response v (1)  is given by equation
 
 L·v   (1)   =s   (9)
 
Note that F NL  is the sum of nonlinear contribution from each device. To calculate distortion caused solely by device d, we omit F NL  in all other devices except in d. Eq. (6) then becomes
 
 L·v+F   NL   (d) ( v )= s   (10)
 
   Compression due to d is computed with v and v (1)  at frequency ω RF +nω LO  (n is an integer): 
                 dB   =     10   ·     log   ⁡     [       v   ⁡     (       ω   RF     +     n   ⁢           ⁢     ω   LO         )           v     (   1   )       ⁡     (       ω   RF     +     n   ⁢           ⁢     ω   LO         )         ]                 (   11   )               
In most analog designs, the circuit functions in a nearly linear region. The nonlinear term F NL   (d)  in eq. (10) is small compared to the linear term L·v and v is not far away from v (1) . Therefore, we can treat F NL   (d)  as perturbation and apply Born approximation to solve v iteratively
   u   (n)   =v   (1)   −L   −1   ·F   NL   (d) ( u   (n-1) ),  (12) 
where u (n)  is the approximated v accurate up to nth order of V RF  for n≧2 and u (1)  is the same as v (1) . Eq. (12) is equivalent to AC or periodic AC analysis with F NL   (d) (v (n-1) ) being the small signal. As shown in eq. (8), evaluation of F NL   (d)  takes a function evaluation of F at device d and its first derivative. No modification is needed in the device model. In compression distortion summary, we solve v up to 3rd order perturbation and use u (3)  as v in eq. (11).
 
   To simplify eq. (12) we define RF harmonic of a multi-tone function ƒ(t) as 
                     f   ω     ⁡     (   t   )       =       ∑   m     ⁢           ⁢       f   ⁡     (     ω   +     m   ⁢           ⁢     ω   LO         )       ⁢     ⅇ       j   ⁡     (     ω   +     m   ⁢           ⁢     ω   LO         )       ⁢   t                   (   13   )               
where ω is the frequency a particular RF harmonic and ƒ(ω+mω LO ) is the Fourier transform of ƒ(t). Eq. (12) can be written in RF harmonic form as
   u   ω   (n)   =v   ω   (1)   −L   −1   ·F   NL,ω   (d) ( u   (n−1) )  (14) 
RF harmonic F NL,ω   (d)  of F NL   (d)  can be obtained by DFT (Discrete Fourier Transform).
 
   For compression distortion, we solve (14) at the first RF harmonic ω=ω RF  to 3rd order.
 
 u   ω     RF     (3)   =v   ω     RF     (1)   −L   −1   ·F   NL,ω     RF     (d) ( u   (2) )  (15)
 
   Notice that as defined in (8) F NL (v) is the nonlinear part of function F(v) and its lowest order polynomial is v 2 . Third order terms in u ω     RF     (3)  equation arise from 2nd and 3rd polynomials in F NL  which involve only RF harmonics at ω RF  and −ω RF  in v (1)  and 2ω RF  and 0 in u (2) . 
   Thus, without losing accuracy we can write eq. (15) as
 
 u   ω     RF     (3)   =v   ω     RF     (1)   −L   −1   ·F   NL,ω     RF     (d) ( v   ω     RF     (1)   +v   −ω     RF     (1)   +u   2ω     RF     (2)   +u   0   (2) )  (16)
 
Eq. (16) shows that to solve u ω     RF     (3)  at third order we only need to solve u 2ω     RF     (2)  and u 0   (2)  in second-order Born approximation, for which we have
 
 u   2ω     RF     (2)   =−L   −1   ·F   NL,2ω     RF     (d) ( v   (1) )  (17)
 
and
 
 u   0   (2)   =−L   −1   ·F   NL,0   (d) ( v   (1) ).  (18)
 
Other RF harmonics in u (2)  contribute to higher order terms, and therefore can be ignored in 3rd order perturbation. (Note that in eq. (16) v ω     RF     (1)  and v −ω     RF     (1)  are complex conjugates so that only one solution must be obtained.)
 
     FIG. 2  shows a method for calculating a compression distortion summary  202  consistent with the above discussion. First, a DC or operating point is determined  204  (e.g., approximated in the time domain [13] or the frequency domain [16]). Next eq. (9) is solved for the first-order linear response v (1)  at ω RF    206 . Note that eq. (9) can be solved in the time domain [13] or alternatively in the frequency domain [16] by standard methods. 
   The process then continues for one or more devices (e.g, component elements d of the system  102 ). First the nonlinear terms are calculated at frequencies 2ω RF  and zero  208 . Next, second-order solutions are calculated for these frequencies  210 . Next the nonlinearities in Eq. (16) are evaluated at ω RF    212 , and the equation is solved for the third-order solution u (3)    214 . This process can be continued for all devices in the system  102  or a limited subset of component devices. 
   Simulation results corresponding to an implementation of this method are shown in  FIGS. 3 and 4 .  FIG. 3  shows a compression distortion summary (in dB) for mixer ne600 (Philips Semiconductors) at −15 dbm at node Pif. Compressions due to all devices and due to each individual device are listed. Results based on time-domain shooting and harmonic-balance implementations are given and they are consistent with each other.  FIG. 4  shows similar results (based on time-domain shooting) for a compression distortion summary (in dB) of a PA (Power Amplifier) at 0.1V at node RFOUT. 
   3. IM2 Distortion Summary 
   In direct conversion mixer design, second-order intermodulation (IM2) between two RF input signals is used to measure circuit distortion. As a second order effect, IM2 does not involve crossing nonlinear terms between devices. It is simply the sum of second order terms in each individual device. This can be demonstrated by second order Born approximation for IM2 at frequency ω 1 -ω 2  where ω 1  and ω 2  are frequencies of the two RF input signals. 
   
     
       
         
           
             
               
                 
                   u 
                   
                     
                       ω 
                       1 
                     
                     - 
                     
                       ω 
                       2 
                     
                   
                   
                     ( 
                     2 
                     ) 
                   
                 
                 = 
                 
                   
                     
                       - 
                       
                         L 
                         
                           - 
                           1 
                         
                       
                     
                     · 
                     
                       
                         F 
                         
                           NL 
                           , 
                           
                             
                               ω 
                               1 
                             
                             - 
                             
                               ω 
                               2 
                             
                           
                         
                       
                       ⁡ 
                       
                         ( 
                         
                           v 
                           
                             ( 
                             1 
                             ) 
                           
                         
                         ) 
                       
                     
                   
                   = 
                   
                     
                       - 
                       
                         L 
                         
                           - 
                           1 
                         
                       
                     
                     · 
                     
                       
                         ∑ 
                         d 
                       
                       ⁢ 
                       
                           
                       
                       ⁢ 
                       
                         
                           F 
                           
                             NL 
                             , 
                             
                               
                                 ω 
                                 1 
                               
                               - 
                               
                                 ω 
                                 2 
                               
                             
                           
                           
                             ( 
                             d 
                             ) 
                           
                         
                         ⁡ 
                         
                           ( 
                           
                             v 
                             
                               ( 
                               1 
                               ) 
                             
                           
                           ) 
                         
                       
                     
                   
                 
               
             
             
               
                 ( 
                 19 
                 ) 
               
             
           
         
       
     
   
   As shown in eq. (19), evaluation of F NL   (d) (v (1) ) takes only linear response v (1)  at second order perturbation. Hence, F NL   (d) (v (1) ) is solely due to nonlinearity in device d and no other device is involved in it. The same argument also holds in Volterra series based perturbation. 
   We apply the adjoint analysis to compute contribution of each device to IM2 at the relevant output nodes. Using transfer function X ω     1     -ω     2   , eq. (19) can be rewritten as inner products of X ω     1     -ω     2    and F NL,ω     1     -ω     2     (d) (v (1) ): 
                     V   out     ⁡     (       ω   1     -     ω   2       )       =       ∑   d     ⁢     (       X       ω   1     -     ω   2         ,     -       F     NL   ,       ω   1     -     ω   2           (   d   )       ⁡     (     v     (   1   )       )           )               (   20   )               
where V out (ω 1 -ω 2 ) is IM2 measured at the output nodes at frequency ω 1 -ω 2 . Note that in this context the transfer function X ω     1     -ω     2    is a projection of the solution operator (e.g., Green&#39;s function) onto a combination of one or more output nodes where V out (ω 1 -ω 2 ) is being measured (e.g. a difference of voltages at two output nodes). This approach reduces the computational cost since X ω     1     -ω     2    is independent of F NL . Once X ω     1     -ω     2    is solved by the linear adjoint analysis, then instead of solving eq. (19) repeatedly, the contribution of element d can be computed by simply calculating the inner product of X ω     1     -ω     2    and F NL,ω     1     ω     2     (d) (v (1) ).
 
   Alternatively, the complete solution of eq. (19) can be calculated before making a projection onto the relevant output nodes. However, the solution of eq. (19) is generally only required at a scalar combination of output nodes where the distortion is being measured. Therefore, calculating the IM2 distortion summary according to eq. (20) is generally preferred. 
     FIG. 5  shows a method for calculating an IM2 distortion summary  502  consistent with the above discussion. First, a DC or operating point is determined  504  (e.g., approximated in the time domain [13] or the frequency domain [16]). Next eq. (9) is solved for the first-order linear response v (1)  at ω 1  and ω 2    506 . As noted above, eq. (9) can be solved in the time domain [13] or alternatively in the frequency domain [16] by standard methods. Next the transfer function X ω     1     -ω     2    is characterized as a projection of the solution operator onto the combination of outputs where the distortion is being measured  508 . Note that this transfer function X ω     1     -ω     2    can be separately calculated or calculated implicitly when determining the inner products as discussed below  512 . 
   The process then continues for one or more devices (e.g, component elements d of the system  102 ). The nonlinearities in eq. (20) are evaluated at ω 1 -ω 2    510 , and the inner products are evaluated  512  for characterizing the distortion V out (ω 1 -ω 2 ). This process can be continued for all devices in the system  102  or a limited subset of component devices. 
   Simulation results corresponding to an implementation of this method are shown in  FIG. 6 .  FIG. 6  shows IM2 distortion summary results in V for a direct conversion mixer. Total IM2 at output nodes and contribution of each individual device are listed. Results from time-domain shooting-based and harmonic balance-based implementations are given and they are consistent to each other. 
   4. Conclusion 
   Additional embodiments relate to an apparatus for carrying out any one of the above-described methods, where the apparatus may include a computer for executing computer instructions related to the method. In this context the computer may be a general-purpose computer including, for example, a processor, memory, storage, and input/output devices (e.g., monitor, keyboard, disk drive, Internet connection, etc.). However, the computer may include circuitry or other specialized hardware for carrying out some or all aspects of the method. In some operational settings, the apparatus may be configured as a system that includes one or more units, each of which is configured to carry out some aspects of the method either in software, in hardware or in some combination thereof. 
   At least some values based on the computed distortion summaries can be saved, either in memory (e.g, RAM (Random Access Memory)) or permanent storage (e.g., a hard-disk system) for later use. For example, values for the left-hand sides of eqs. (11) and (20) can be used directly as distortion summaries. Alternatively, values for components of the right-hand sides of these equations can also be saved for later use. And similarly related values (e.g., averages) can also be saved depending on the requirements of the operational setting. 
   Additional embodiments also relate to a computer-readable medium that stores (e.g., tangibly embodies) a computer program for carrying out any one of the above-described methods by means of a computer. The computer program may be written, for example, in a general-purpose programming language (e.g., C, C++) or some specialized application-specific language. The computer program may be stored as an encoded file in some useful format (e.g., binary, ASCII). 
   Although only certain exemplary embodiments of this invention have been described in detail above, those skilled in the art will readily appreciate that many modifications are possible in the exemplary embodiments without materially departing from the novel teachings and advantages of this invention. For example, aspects of embodiments disclosed above can be combined in other combinations to form additional embodiments. Accordingly, all such modifications are intended to be included within the scope of this invention. 
   5. References 
   The following references are related to the disclosed subject matter:
     [1] J. Bussgang, L. Ehrman, and J. Graham. Analysis of nonlinear systems with multiple inputs.  Proc. IEEE,  62:1088-1119, August 1974.   [2] S. Chisholm and L. Nagel. Efficient computer simulation of distortion in electronic circuits.  IEEE Trans. on Circuits and Systems,  20(6):742-745, November 1973.   [3] P. Feldmann, B. Melville, and D. Long. Efficient frequency domain analysis of large nonlinear analog circuits.  Proc. Custom Integrated Circuits Conf ., pages 461-464, 1996.   [4] J. Haywood and Y. L. Chow. Intermodulation distortion analysis using a frequency-domain harmonic balance technique.  IEEE Trans. on Microwave Theory and Techniques,  36(8):1251-1257, August 1988.   [5] Y. Kuo. Distortion analysis of bipolar transistor circuits.  IEEE Trans. on Circuits and Systems,  20(6):709-716, November 1973.   [6] T. Lee.  The Design of CMOS Radio - Frequency Integrated Circuits . Cambridge University Press, Cambridge, UK, 1998.   [7] B. Leon and D. Schaefer. Volterra series and picard iteration for nonlinear circuits and systems.  IEEE Trans. on Circuits and Systems,  25(9):789-793, September 1978.   [8] P. Li and L. Pileggi. Efficient per-nonlinearity distortion analysis for analog and rf circuits.  IEEE Trans. on Computer - Aided Design of Integrated Circuits and Systems,  22(10):1297-1309, October 2003.   [9] S. Maas.  Nonlinear Microwave Circuits . Artech. House, Norwood, Mass., 1988.   [10] S. Narayanan. Application of volterra series to intermodulation distortion analysis of transistor feedback amplifiers.  IEEE Trans. on Circuits and Systems,  17(4):518-527, November 1970.   [11] B. Razavi.  RF Microelectronics . Prentice Hall, Upper Saddle River, N.J., 1998.   [12] J. Sakurai.  Modern Quantum Mechanics . The Benjamin/Cummings Publishing Company, Inc., Menlo Park, Calif., 1985.   [13] R. Telichevesky, K. Kundert, and J. White. Efficient AC and noise analysis of two-tone RF circuits.  Proc. of  33 th DAC  pages 292-297, June 96.   [14] P. Wambacq, G. Gielen, P. Kinget, and W. Sansen. High-frequency distortion analysis of analog integrated circuits.  IEEE Trans. on Circuits and Systems II,  46(3):335-345, March 1999.   [15] F. Yuan and A. Opal. Distortion analysis of periodically switched nonlinear circuits using time-varying Volterra series.  IEEE Trans. on Circuits and Systems I,  48(6):726-738, June 2001.   D. Long, R. Melville, K. Ashby, and B. Horton. Full-chip Harmonic Balance.  Custom Integrated Circuits Conference,  1997 , Proceedings of the IEEE  1997, 5-8 May 1997, pp. 379-382.