Simulation Of Microwave And Millimeter-wave Oscillators

Preparing to load PDF file. please wait...

0 of 0
100%
Simulation Of Microwave And Millimeter-wave Oscillators

Transcript Of Simulation Of Microwave And Millimeter-wave Oscillators

SIMULATION OF MICROWAVE AND MILLIMETER-WAVE OSCILLATORS, PRESENT CAPABILITY AND FUTURE DIRECTIONS 1
MICHAEL B. STEER
High Frequency Electronics Laboratory Department of Electrical and Computer Engineering North Carolina State University, Raleigh NC 27695-7911, U.S.A. +919-515-5191 FAX: +919-515-5523 email: [email protected]
ABSTRACT
The current state-of-the art of oscillator simulation techniques is presented. Candidate approaches for the next genertion of oscillator simulation techniques are reviewed. The method is presented which uses an efficient and robust convolution-based procedure to integrate frequency-domain modeling of a distributed linear network in transient simulation. The impulse response of the entire linear distributed network is obtained and the algorithm presented herein ensures that aliasing effects are minimized by introducing a procedure that ensures that the interconnect network response is both time-limited and band-limited. In particular, artificial filtering to bandlimit the response is not required.
I. Introduction
Large signal simulation of microwave oscillators is necessary to provide steady-state characterization of oscillator performance. Such quantities as power and harmonic content information are then readily available. This is particularly important in achieving first pass successful design of monolithicly integrated oscillators. Circuit simulation of microwave oscillators by the method of harmonic balance is reasonably mature with several commercial products available and used on a regular basis and been adapted to some rather unusual applications, e.g. [44]. However large signal oscillator analysis in the time domain using programs such as SPICE [20] enables the build-up of oscillations to be observed. In spite of being time-consuming and the difficulty of determining the time at which steady state is obtained, time-domain simulation techniques have the ability to predict the start-up of oscillation in addition to the frequency of oscillation and non-steady-state behavior. There a many difficulties in applying transient analysis techniques to distributed circuits but these are gradually being addressed. The near future will see a rapid development of these techniques and will be used regularly in microwave oscillator simulation. In this paper we review the current state of microwave oscillator simulation using the harmonic balance approach and describing function methods. We then consider the current state of transient microwave oscillator simulation and focus on transient simulation techniques that have potential for microwave oscillator simulation. We present a new method of transient simulation which integrates the distributed nature of microwave circuits into a transient simulator by using convolution of the impulse response of the linear circuit. Particular attention is given to reducing aliasing effects in deriving the impulse response and in handling of high Q linear circuits.
II. Harmonic Balance Techniques
The obvious approach to free-running oscillator analysis is to use the harmonic balance equations developed for the circuit and to include the oscillation frequency as an additional optimization variable. This method has been used by a number of workers [5], [6]. Generally,
1Keynote talk — Workshop on Nonlinear Microwave CAD, Germany, October 1992.

one of the variables that would be used as an optimization variable in examining a nonautonomous circuit2 is eliminated, for example, by setting the phase of a voltage or current to zero. Usually, with this approach, the simulated results tend to converge to a degenerate solution [7] (e.g. all currents equal to zero is also a solution of Kirchhoff’s current law which is the basis of the harmonic balance equations.), or else the initial setting of the oscillating frequency must be very close to the final result [6].
The degenerate solution can be avoided by incorporating additional criteria in the system objective function. This was done by Sterzer [8], in the early 1960’s, in calculating the output power of a GaAs tunnel diode oscillator by incorporating the Kurokawa oscillation condition [9]. Also, in the early 1980’s, Solbach [10] working with a Gunn diode oscillator and Bates [11] examining an IMPATT diode oscillator, predicted the frequency and output power by solving multifrequency forms of Kurokawa’s oscillation condition [12] using frequency-domain powerseries analysis techniques. However, the work on single diode oscillator circuits cannot be directly extended to general nonlinear oscillator circuits as unique characteristics of the diode or specialized analysis strategies are employed. In recent work from our group we presented method of autonomous circuit simulation which incorporates the Kurakawa conditions in a general purpose harmonic balance framework [40]. Kurokawa’s oscillation condition is coupled with the modified nodal admittance form of the circuit equations to avoid degenerate solutions. The algorithm has been implemented using both harmonic balance and frequencydomain spectral balance techniques. A range of possible oscillation frequencies is specified and a ‘random walk’ frequency search conducted until a convergent starting point was found. This is an inefficient procedure and has been modified at Compact Software, Inc. [1] to use an initial linear circuit bifurcation search to obtain a valid initial frequency guess. For low level oscillations and in simulating high Q oscillator circuits this initial frequency guess is extremely accurate.
Rizzoli et al. [2], [3] proposed a method (implemented in [4]) based on the harmonic balance technique for oscillator synthesis. The oscillation frequency is fixed while one circuit parameter is optimized to ensure that the harmonic balance equations are satisfied at that frequency. This method is numerically efficient and yields well-defined and accurate results. However, this method is not directly amenable to free-running oscillator analysis as the frequency is fixed and a degree of freedom, such as a tuning element or the load impedance, varied until harmonic balance is achieved. By repeating this process for a number of frequencies a curve of frequency versus the degree of freedom is obtained. In this manner oscillators can be analyzed without involving autonomous circuit simulation. This is particularly advantageous as the time required to simulate the linear subcircuit is typically much greater than the time required to simulate the nonlinear subcircuit and the frequency searching approach requires very many linear circuit evaluations. In contrast the synthesis approach requires a single linear subcircuit evaluation for each steady-state point found.
In a hybrid of the frequency searching and synthesis techniques, Elad, Madjar and BarLev [42] begin with a nonautonomous harmonic balance analysis of the oscillator circuit operating under small signal conditions using, as an intial guess, an oscillation frequency determined from small signal analysis. As the amplitude of oscillations is allowed to increase the linear circuit is adjusted to prevent the oscillation frequency from deviating from its small signal value due to small signal effects.
Frequency domain techniques for the simulation of microwave oscillators have most recently been reviewed by Rizzoli, Neri, Ghigi and Mastri [48] where frequency domain oscillator design approaches and stability analysis techniques are considered and are not considered here. aspects.
III. Describing Function Techniques
Many workers have implemented noniterative nonlinear analyses of free-running oscillators using describing function techniques [13], [14] (commonly used in nonlinear control system analysis) and functional expansions (generally based on Volterra series techniques
2A non-autonomous circuit is one in which the frequencies of signals are determined by signal sources.
For example, an amplifier is a non-autonomous circuit.
2

[15], [16], [17], [47] but also using specialized functional expansions [18]). In these techniques, the system equations and oscillation criteria are combined to yield a set of algebraic equations which can be solved recursively. However, these methods are restricted to weakly nonlinear oscillators. Cheng and Everard [19] used a spectral balance approach to analyze a microwave FET oscillator. They used a Volterra series technique to generate the harmonic components but, since they used power series expansions of the nonlinear element characteristics, their technique is not restricted to weakly nonlinear oscillators. In their method they converted an oscillator into a one-port network by making a break somewhere in the circuit. This leads to the oscillation criteria that the impedance looking into this port is zero at the fundamental and at all harmonics. A relaxation algorithm is used to solve for this condition. The relaxation algorithm, however, has poor convergence properties and the use of Volterra-series-based nonlinear analysis is restrictive (limiting the analysis to one dimensional nonlinearities and, unless power-series-like descriptions are available, to weakly nonlinear oscillators).
IV. Transient Techniques
All of the above techniques assume that a periodic steady-state solution of the system equations exist and then proceed to derive it. In contrast, large signal oscillator analysis in the time domain using programs such as SPICE [20] enables the build-up of oscillations to be observed [21] as well as chaotic behavior [35]. In spite of being time-consuming and the difficulty of determining the time at which steady state is obtained, time-domain simulation techniques have the ability to predict the start-up of oscillation in addition to the frequency of oscillation and non-steady-state behavior (e.g. chaotic behavior). It is also easier to incorporate physical device models (e.g. those described by coupled partial differential equations or electron statistics) in time-domain simulations [22], [23], [24].
The Kurokawa oscillation criterion was incorporated in the system error function in an harmonic balance program in which the oscillating frequency is used as an independent variable¿ The problem was formulated so that the resulting set of equations can be solved using the Newton method to speed convergence. The method also has good convergence properties so that the initial setting of the frequency variable need not be very close to the actual oscillating frequency.
While the steady-state performance of microwave oscillators is of considerable interest and enables oscillators to be synthesized, an oscillators transient performance is critical as such as attributes as chaotic behavior, start-up performance can only be predicted via transient simulation. In some cases it is possible to acquire some measure of performance of some timedomain based phenomena by approximating the transient behavior with a pulse sequence which has a Fourier description of manageable size. An example is the analysis of injection locking behavior in the recent work of number of a number of groups [38], [51], [53], [54], [55], [56], [52] [57]. Nonlinear simulation of power combiner is a particular case where transient analysis reigns supreme as prediction of transient and multifrequency behavior is generally of much greater importance than determining single frequency behavior. An harmonic balance simulator would to to assume that each oscillator is frequency locked, then searches for signal amplitudes, phases, and the system oscillation frequency. The derivation of this linear model should be physically based so that design parameters can be modified to improve power combiner performance.
York and Compton have studied weakly coupled oscillators [50]. Individual oscillators were characterized by their free running frequency, free running oscillation amplitude, and the oscillator external Q. Starting with an equation similar to Adler’s equation, an expression was derived from which the array oscillation frequency and phase distribution could be numerically determined. The stability of various phase distributions was investigated. A simplified example of four identical oscillators was analyzed and it was concluded, for free space coupling, that the oscillators must be separated by distances equal to multiples of one wavelength. This spacing is generally undesirable as it produces of radiation pattern grating lobes and it precludes a high area density of oscillators. Observations were made that array end effects are significant for small arrays, and that tight tolerances are necessary to achieve synchronization of weakly coupled oscillators. Experiments were conducted using
3

4 × 4 arrays of Gunn diodes and MESFET’s with rectangular microstrip patch antennas. Synchronization was achieved by use of an open, low Q, rectangular cavity. Modification of the theory is required to include the effects of the resonator on array behavior.
Stephan has used an analysis similar to Kurakawa’s to derive a set of time domain differential equations which describe the behavior of a network of lumped element coupled oscillators [68]. The nonlinear equations were solved numerically assuming that the active devices were purely conductive. The results were not extended to distributed circuits however. Stephan and Young [69] have modeled the mode selection and frequency pulling of two Gunn diode oscillators which were coupled through direct radiation.
The fundamental difficulty encountered in integrating distributed microwave element models in a transient circuit simulator arises because circuits containing nonlinear devices or time dependent characteristics must be characterized in the time domain while transmission lines and other distributed structures are best simulated in the frequency domain. Djordjevic, Sarkar and Harrington [25, 26, 29] and Schutt-Aine and Mittra [27, 28] were among the first to present general techniques for simulating such systems. These methods were based on developing the impulse response of the linear circuit and incorporating this characterization in a transient circuit simulation using a convolution operation. The subject described in these publications is concerned with the simulation of interconnects, particularly coupled lines, but it does form the basis for a more generally applicable analyses. More recently the asymptotic waveform evaluation (AWE) method has been developed which approximates the linear subcircuit by a limited order pole-zero description and so currently has the potential for handling essentially lumped microwave circuits. Currently much of this work is directed at the incorporation of transmission lines in high speed digital system simulation.
IV.A Asymptotic Waveform Evaluation
In the special case of no nonlinear devices the entire network can be modeled in the frequency-domain and the transient (albeit periodic) response of the circuit determined by inverse Fourier transforming the appropriate voltage and current phasors. However if an s-domain transfer function is first fitted to the frequency domain response then the response of an arbitrary network can be incorporated in transient analysis by using the asymptotic waveform evaluation technique pioneered by Rohrer et al. [66]. The extension of this technique to arbitrary interconnects was performed by Nakhla et al. [61] [63] [62] and has also been extended to more general distributed circuits [67]. Nakhla et al. show that given with the frequency-domain characteristics of the distributed circuit it is a simple matter to automatically generate the transient response function through inverse Laplace transformation [60]. The inverse Laplace transform operation is segmented in such a way that it can be incorporated into the modified nodal admittance matrix of the associated discrete circuit model. Thus very efficient simulation of networks for which only frequency domain information is available can be achieved. The accuracy of this method is limited in the first place by the modeling of the interconnect by a finite number of poles. Unfortunately the response of a distributed interconnect approximates a pole-zero response only poorly. The big advantage is that the technique should be marginally faster than other transient analysis techniques as the response of the distributed network is described by just a few poles and their residues rather than by 1000 or so discrete frequency values.
IV.B Convolution Techniques
In one of Djordjevic et al’s techniques [29], the frequency-domain admittance (Y) parameter description of a distributed network is converted to a time-domain description using a Fourier transform. This time-domain description is then the Dirac delta impulse response of the distributed system. Using the method of Green’s function the system response is found by convolving the impulse response with the transient response of the terminating nonlinear load. Normally this requires that the impulse response be extended in time to include many reflections. While this technique can handle lossy coupled networks, a difficulty arises as the y parameters of a typical multiconductor array have a wide dynamic range. For a low loss, closely matched, strongly coupled system, the Y parameters describing the coupling
4

mechanism approach zero at low frequencies and become very large at high frequencies. Conversely the transmission and self-admittance Y parameters approach infinity at DC and zero at resonance frequencies. Both numerical extremes are important in describing the physical process of reflections and crosstalk. The dynamic range of the time-domain solution is similarly large and values close to zero are significant in determining reflections and crosstalk. Consequently, aliasing in the frequency-domain to time-domain transformation can cause appreciable errors in the simulated transient response.
An alternative formulation that avoids the dynamic range problems is that of Schutt-Aine and Mittra [28] who use a scattering parameter formulation to consider a parallel coupled transmission line system. This approach offers good computational stability and efficiency. If a scattering (S) parameter is small in magnitude, it corresponds to a small coupling, transmission, or reflection component. If the S parameter is near one in magnitude it corresponds to a significant reflection, transmission, or coupling component. Correspondingly, in the time domain, components with a small magnitude are less significant. Consequently dynamic range errors in the frequency- to time-domain transformation are not so important.
The Laplace transform-based techniques, and the AWE techniques extend the range of transient analysis of distributed systems to the high 100’s of MHz. For systems operating at GHz frequencies and above, convolution-based techniques appear to be the only viable approaches as they use the frequency-domain characterization of the distributed network with no modification. Convolution-based techniques inverse transform (e.g., using Fourier, Hankel or Laplace transformation) the frequency domain characterization of the distributed network to obtain its impulse response of the network. Transient analysis is achieved using convolution of this impulse response with the voltages at the terminals interfacing the interconnect network with the remainder of the circuit. Generally low pass filtering of the frequency-domain parameters is required to ensure that they are band-limited or else aliasing effects may result in non-causal response and numerical convergence problems. As well, the impulse response derived must be time limited and this is achieved by terminating the linear network in the average characteristic impedance of the transmission lines in the network. Even with all of these considerations performance has not been completely accurate as non-causal effects, aliasing , and Gibb’s phenomenon often lead to numerical convergence problems. The most serious of these problems is non-causal effects which are due to the artificial filtering required to band-limiting frequency-domain parameters. This is a problem with all convolution-based techniques proposed to date. The various methods differ in the way in which the frequency-domain parameters are derived, the attention given to limiting the dynamic range of the response presented to the inverse Fourier transformation, and actions taken to limiting the time response and band-limit the frequency response of the network. These are required to minimize aliasing effects upon Fourier transformation or numerical inversion of the Laplace transform. All of the techniques can be adapted to function with a transient simulator. The techniques develop a green’s function relating the voltage or current response at one port of the distributed network to the voltages at past times and at all ports of the network. Djordjevic et al. [64] limit the time response by introducing a augmentation network to terminate transmission lines in their approximate characteristic impedance. The inverse Fourier transform of the y parameters of this augmented network is calculated and the effect of the augmented network removed in transient analysis.
The convolution-based technique presented here uses naturally band-limited parameters
IV.C Shooting Methods
Shooting methods are a specialized form of time-domain analysis applicable to circuits with strictly periodic excitation [70] [71], [72], [73], [74], [75]. They bypass the transient response altogether and are advantageous in situations that would require many iterations for the transient components to die out, if direct integration methods were used. It is assumed that the nonlinear circuit has a periodic solution and that the solution can be determined by finding an initial state such that transients are not excited. If x(t) is the set of state variables obtained by a time-domain analysis, the boundary value constraint for periodicity is that x(t) = x(t + T ) where T is the period. A series of iterations at time points between t and t + T can be performed for a given set of initial conditions, and the condition for
5

periodicity checked. Thus in the shooting method the problem of solving the state equations is converted into the two-point boundary value problem [71, 72]

x(0) = x(T )

T

x(T ) = f(x, τ )dτ + x(0)

(1)

0

where T is the period such that

x(t) = x(t + T ).

if x(t) is indeed the solution. If x(t) = x(t + T ) then a new set of initial conditions can then be determined using a gradient method based upon the error in achieving a periodic solution. Once the sensitivity of the circuit to the choice of initial conditions is established in this way, a set of initial conditions that establishes steady-state operation can be determined; this set is, of course, the desired solution. This iterative procedure can be implemented using the Newton’s method iteration

Xk+1(0) = Xk(0) − I − ∂Xk(T ) −1 [Xk(0) − Xk(T )] (2) ∂Xk(0)

where the superscripts refer to iteration numbers and Xk(T ) is found by integrating the
circuit equations over one period from the initial state Xk(0). To begin the analysis the period (T ) is determined and the initial state (X(0)) is esti-
mated. Using these values, the circuit equations are numerically integrated from t = 0 to t = T and the necessary derivatives calculated. Then, the estimate of the initial state is updated using the Newton iteration or alternative optimization scheme [74] may be used as calculation of the derivatives required for the Newton iteration can be time consuming for large problems. This process is repeated until X(0) = X(T ) is satisfied within a reasonable tolerance.
The computation becomes further complicated when transmission lines are present, because functional initial conditions are then required to establish the initial conditions at every point along the line (corresponding to the delayed instants in time seen at the ports of the line).
Multiple shooting algorithms have been developed which enable a circuit to be partitioned into two component networks and a different shooting algorithm applied to each part [76]. This algorithm is the basis of the FATE method — FATE being an acronym for Frequency and Time Evaluation of an oscillators periodic steady state [31], [37] In FATE the resonant tank circuit is partitioned from the rest of the circuit which now may comprise just the active device and a few linear, lumped device parasitics. An unrestricted standard shooting algorithm can be applied to the nonlinear lumped subcircuit but the tank circuit is assumed to have a single frequency response. The difference between the two sets of waveforms is then minimized using the Newton method.

IV.D Alternative Methods

There are many other transient analysis techniques that could lend themselves to oscillator simulation. Work that deserves particular mention is that of Sobhy et al. [33], [36], [32] who have developed state equation based analysis of distributed circuits. The analysis takes place using direct integration of the state equations and a method has been developed to handle distributed elements to some extent. Similarly a traditional time-marching scheme allows numerical device models to be used in conjunction with restricted distributed circuits. [24]. [49]

VI. Simulation of a Quasioptical Oscillator

6

SOURCE ARRAY
a^ y a^x a^ z
PLANAR REFLECTOR
d D

PARTIALLY TRANSPARENT SPHERICAL REFLECTOR

Figure 1: Quasi-Optical Power Combiner Configuration.

Ca

ia = f(ia,va)

V0a + va I0

Cd I0 ie V0d + vd

Figure 2: IMPATT diode model.
An example of the use of convolution-based transient simulation is the analysis of the quasioptical power combining oscillator shown in Fig. 1. The purpose of this work is to investigate strange and currently unpredictable phase-locking behavior. Current design and operation concerns of quasi-optical combining systems largely involve locking behavior. Much of the concept behind quasi-optical power combining relies on global locking of the various oscillating elements. However maximizing the power extracted from the system requires that the elements be placed relatively closely to each other. This increases the nearest element interactions so that local locking rather than global locking becomes dominant. However some local locking is required to ensure that the sources initially lock. To model these effects transient simulation is required. We have used an existing convolution-based transient simulator [58] and [65] to successfully analyze a simple model of a quasi-optical cavity oscillator with a single IMPATT diode source. The IMPATT diode model used was that of Goeller and Kaertner [45], see Fig. 2. For the simulation their parameters were used. The simulation strategy uses the impulse response of the linear circuitry to interface the frequency-dependent cavity model to the nonlinear oscillating element model. Interfacing is achieved through convolution. The major problem of using a convolution based approach is handling the the cavity. The unloaded Q of the cavity is about 300,000 so special techniques were developed to permit the cavity to be described by its impulse response. This is described in the next section.
VI.A Convolution Method Development
In convolution-based analyses the distributed network can be partitioned into a linear
7

LINEAR SUBCIRCUIT

iN v 1
iN v N

NONLINEAR SUBCIRCUIT

Figure 3: N port distributed network partitioned into linear and nonlinear subcircuits.

v1

LINEAR

i1

SUBCIRCUIT [Y] vN

AUGMENTED LINEAR i

SUBCIRCUIT

N

[Y ]

v1

v1

i 1

i1

vN

vN

iN

iN

COMPLEMENTARY NETWORK

NONLINEAR SUBCIRCUIT

Figure 4: Augmented network.

subcircuit and a nonlinear subcircuit, Fig. 3. Each time-domain y parameter (being the Fourier transform of the corresponding frequency-domain y parameter) is then the response of the current with respect to an impulse voltage at one of the ports of the network, i..e. yij(t) = ii(t)/vj(0). Other parameters such as S parameters have been used principally because these parameters intrinsicly have a limited range (between 0 and 1). Reflections in the circuit are reduced by terminating each port in its reference impedance so ensuring that there are no reflections at the ports of the distributed network. This is illustrated by considering the augmented network of Fig. 4 with ports terminated in Thevenin equivalent loads (or sources) each with reference impedance Zm. The procedure follows the devlopment in [58] and [65].
VI.B Quasioptical Cavity Model
The cavity model is based on an impedance matrix of a dipole array in a quasi-optical resonator [30]. This model develops an impedance matrix of the cavity with respect to the diodes in the resonator, see Fig. 5. The impedance matrix is developed by considering the fields in the cavity and a reflector using the field distribution shown in Fig. 6. The reflector is described by its reflectance R which is -1 for a perfectly reflecting surface. Typically there is very little loss at the reflector and the cavity has an extremely high unloaded Q. In order for us to perform an FFT of the frequency domain data as required in using the convolution-based transient analysis method we must develop a technique for reducing the Q of the cavity in the frequency domain description and reconstruct it during transient analysis. One mechanism for doing this is to introduce an additional port at the reflector as shown in Fig. 7. So that the previous work on cavity modeling can be reused we develop the model of the cavity with the additional port by placing loads (i.e. different R values) at the reflector port and recalculating the impedance matrix of the reflector. This was
implemented by deriving the equivalent impedances of the mirror (= µ/ (1 + R)/(1 − R)
8

DIODE 1 DIODE N

CAVITY

Figure 5: N diode cavity model.

a^ y a^ x
a^ z a Emn mn mn
- mn amnHm- n

PARTIALLY TRANSMITTING REFLECTOR

FILAMENTARY ANTENNA

- a E + (1 + ) mn mn mn mn + - mnamnHmn(1 + mn)

z=0

z = D

Figure 6: Cross-Section of Quasioptical Resonator Showing Fields, E2, Established by Radiating Current Elements.

CAVITY

LOAD

DIODE PORT

REFLECTOR PORT

Figure 7: Quasioptical cavity with a single IMPATT diode source.

9

FREQUENCY DOMAIN ANALYSIS SUBCIRCUIT

CAVITY

AUGMENT NETWORK

CORRECTION NETWORK

REFLECTOR SURFACE

DIODE PORT

REFLECTOR PORT
TRANSIENT ANALYSIS SUBCIRCUIT

Figure 8: Single diode analysis structure.

400 300 200 100
50

Mag. Admittance

50.2

50.4

50.6

50.8

Freq.(GHz) 51

Figure 9: Admittance calculated looking into the quasioptical cavity.

and then using a standard three load calibration procedure. The dynamic-range-limiting and band-limiting requirements for the frequency domain parameters necessitates using an augmentation network at the reflector terminals so that the circuit simulated is as shown in Fig. 8. The augmentation network is included in the frequency domain characterization of the cavity and so its effect appears in the impulse response of the circuit. Subsequently the effect of the augmenting circuit is removed during transient analysis through the use of a complementary network. Note that the complementary network and reflector surface element are linear and resistive so that it would seem that they could have been absorbed into the cavity model. However in effect what we have done is introduce an additional state variable in the transient analysis which enables the original low loss and high Q nature of the quasioptical cavity to be reconstructed during transient analysis.
The distance between the flat reflector and the parabolic reflector was 6.02 cm which was approximately equal to the focal length of the reflector. Simulation of the cavity yields a very rich set of resonances as shown in the admittance plot with a perfectly reflecting mirror shown in Fig. 9 Experimentally it has been found that multifrequency oscillation at several of the modes may occur at high bias levels. Since this is unforseen performing an harmonic balance analysis on this circuit would not predict this effect.
VI.D Results
The single diode quasioptical cavity was simulated using a simplified distributed model of the cavity owing to problems in correctly modeling the input impedance of the cavity at DC and at very low microwave frequencies. After thresholding the impulse response the impulse response of the cavity essentially consisted of a y11 free space admittance term at time 0 and a y12 impulse response at 0.1035 ns corresponding to the transit time of the cavity which after a round trip almost exactly canceled the response at time zero. The current response is given in Fig. 10.
VII. Conclusion

10
NetworkFrequencyCircuitAnalysisCavity