This page contains the full text of the accepted version of the paper published at the 2018 International Microwave Symposium. The content was modified slightly to better fit the format of this website

The impedance of a microwave circuit has an infinite number of poles due to the distributed elements. This complicates locating those poles with a rational approximation. In this paper, we propose an algorithm to locate the unstable poles of a circuit with distributed elements. The proposed method exploits the fact that a realistic circuit can only have a finite number of unstable poles. We first determine the unstable part whose poles coincide with the unstable poles of the circuit. A rational approximation of the unstable part is used to estimate the unstable poles.

1. Introduction

Having the ability to trace a circuit’s poles as a function of the circuit parameters is a useful design tool. Pole tracking techniques have been used for the design of oscillators^{1 2}, to stabilise power amplifiers^{3 4 5} and during the design of frequency dividers^{1 6}. The core of a pole tracking tool is a robust automatic algorithm to estimate the poles of a circuit.

To determine the poles of a circuit, a two-step procedure is followed. First, an impedance $Z\!\left(j\omega\right)$ of the circuit is determined at discrete frequency points between $0$ and $f_{\mathrm{max}}$ with an AC simulation (Fig. 1). Then, the poles of $Z\!\left(j\omega\right)$ are determined in a post-processing step.

Figure 1: $Z\!\left(j\omega\right)$ is determined by connecting a small-signal current source $I_{ex}$ in parallel to a node in the circuit and by determining the voltage response $V\!\left(j\omega\right)$ of the circuit with an AC simulation.

In lumped circuits, $Z\!\left(j\omega\right)$ is a rational function, so a rational approximation can be used to determine the circuit poles. The impedance presented by a circuit with distributed elements, like transmission lines, is non-rational but still meromorphic ⁷. This makes estimating the poles of a microwave circuit more difficult. A good fit of $Z\!\left(j\omega\right)$ can be obtained with a high-order rational function, but the rational approximation will contain spurious poles that do not correspond to poles of the underlying function ⁸. To circumvent this problem, another approach proposed in ⁹, is to compute low-order rational approximants of the circuit’s response restricted to small frequency intervals. This local version of the rational approximation scheme, yields precise estimates of poles when these are close enough to the imaginary axis.

In this paper, we propose an algorithm that can estimate the unstable poles of a circuit without performing a rational approximation of a non-rational function. The proposed technique exploits the fact that the equilibrium solution of a realistic circuit ¹⁰ can only have a finite amount of poles in the right half-plane ¹¹.

When $Z\!\left(j\omega\right)$ is square integrable, it can be split uniquely into two meromorphic parts: one part in $H_{2}$, the Hardy space of the complex right half-plane, and another part in $\overline{H}_{2}$, the orthogonal complement of $H_{2}$. The projection of $Z\!\left(j\omega\right)$ onto a basis of $\overline{H}_{2}$ is defined as its “unstable part” $Z_{\mathrm{unstab}}\!\left(j\omega\right)$ because its poles coincide with the unstable poles of $Z\!\left(j\omega\right)$.

There are only a finite amount of unstable poles, so $Z_{\mathrm{unstab}}\!\left(j\omega\right)$ is rational. To reliably estimate the unstable poles of $Z\!\left(j\omega\right)$, we propose to compute its unstable part and to approximate this unstable part with a rational function.

In¹², it was shown that estimates of the unstable part can be obtained by computing the Fourier coefficients of a pre-filtered version of the impedance function. In this paper, the method is extended so as to yield estimates of the unstable poles. In order to do so:

• A Hankel-operator-based rational approximation technique is presented (Section 2.3).
• Fourier coefficients are obtained via a rational interpolation scheme that reduces the errors due to sampling (Section 2.4).
• The spectral range of the filtering function is reduced to improve the procedure’s robustness (Section 2.5).

After detailing the approach, we run it on two examples: a time-delayed Chua’s circuit (Section 3) and on simulations of a microwave oscillator (Section 4).

2. Algorithm

The response $Z\!\left(j\omega\right)$ is first transformed to the unit circle (Section 2.1) to allow for a simple determination of its unstable part by means of its Fourier series (Section 2.2 & 2.4). With its Fourier expansion, the unstable poles can be recovered by means of a singular value decomposition of a truncated version of the system’s Hankel matrix (Section 2.3). These steps assume that $Z \left(j\omega\right)$ is known everywhere while it is only known on a finite interval. To deal with finite data, the Fourier coefficients are convolved by the Fourier coefficients of the filtering function (Section 2.5).

2.1. Transform to the unit disc

A Mobius transform is used to transform $Z\!\left(j\omega\right)$ from the complex plane to the unit circle:

where $\alpha=\frac{2\pi f_{\mathrm{max}}}{\left(1+\sqrt{2}\right)}$. This transform maps $\omega=0$ to $\theta=\pi$ and $f_{\mathrm{max}}$ to $\theta=\frac{\pi}{4}$. The complex right half-plane is mapped to the inside of the unit circle (Fig. 2). We will use $\breve \cdot$ to indicate functions on the unit disc.

2.2 Determine the unstable part

On the unit circle, a basis for $\overline{H}_{2}$ is $z^{-k},k\in\mathbb{N}$ ¹² so we compute the Fourier series of $\breve Z\left(\theta\right)$ to obtain its unstable part

In ¹², the Fast Fourier Transform (FFT) is used to compute the $f_{k}$. This introduces an aliasing error in the method. To avoid aliasing, we compute the Fourier coefficients with the help of an interpolation function as explained in Section 2.4.

2.3. Estimating unstable poles

We use the theory introduced by Adamyan, Arov and Krein ¹³ to recover the location of the unstable poles on the unit disc. It uses the following infinite Hankel matrix which contains the negative Fourier coefficients of $\breve Z\left(\theta\right)$:

In practice, we truncate the infinite $\Phi_{\,}$ to $\Phi_{N}=\Phi_{1..N,1..N}$. $N>50$ has yielded good results in all examples we tried. By Kronecker’s theorem for a rational symbol with $P$ poles, $\Phi_{N}$, is of rank $P$. We compute the rank using the Singular Value Decomposition of $\Phi_{N}$

Ideally the number of non-zero singular values of $\Phi_{N}$ equals the number of poles of the function $Z_{\mathrm{unstab}}\!\left(j\omega\right)$. In practice due to numerical errors, the first $P$ singular values are large when compared to the others. This property is used to determine the number of unstable poles $P$ of $Z$.

Once $P$ has been determined, we can determine the circuit poles using the Principal Hankel Components (PHC) algorithm ¹⁴: We split the matrices $U$, $S$ and $W$ into parts

such that $S_{P}$ is a diagonal matrix of size $P\times P$ which contains the singular values, $U_{P}$ is of size $N\times P$ and $W_{P}^{\prime}$ is size $P\times N$. With $U_{p}$ and $S_{P}$, we compute an estimate of the observability matrix of the approximation:

From the observability matrix $\mathcal{O}$, we can compute the state matrix $A$ using

where $\mathcal{O}_{l}$ and $\mathcal{O}_{f}$ are formed from $\mathcal{O}$ by deleting the last and first row respectively. $\cdot^{\prime}$ indicates the complex conjugate. The eigenvalues of $A$ correspond to the location of the unstable poles.

2.4. Computing the Fourier coefficients

The Fourier coefficients needed to construct $\Phi_{N}$ are given by:

The function $\breve Z\left(\theta\right)$ is known on a set of $M$ sample points $\theta_{m}$. This calls for the use of an interpolation procedure for a precise evaluation of (4). When the impedance $\breve Z\left(\theta\right)$ contains strong resonances, classical spline interpolation fails to model the response’s stiffness. We therefore introduce a rational interpolation scheme by defining $\breve I_{m}\left(\theta\right)$ between the points $\theta_{m}$ and $\theta_{m+1}$ as:

where $a_{m}$, $b_{m}$, $c_{m}$ and $d_{m}$ are determined using the constraints:

The final constraint ensures that the interpolated function has a continuous first derivative at the interpolation points.

When the rational interpolation models $\breve I_{m}\left(\theta\right)$ are known, the integral (4) is computed with Matlab’s built-in quadrature method ¹⁵.

2.5. Filtering to suppress edge effects

The simulated impedance is only known for frequencies below $f_{\mathrm{max}}$. Above that frequency, we assume that $Z\!\left(j\omega\right)=0$. This sudden jump to zero at $f_{\mathrm{max}}$ will have a strong unwanted influence on the obtained stable and unstable part ¹². To suppress this influence, $\breve Z\left(\theta\right)$ is multiplied by a low-pass filter with its stop-band starting at $f_{\mathrm{max}}$:

In ¹², an elliptic filter on the complex plane is used. On the unit disc, this elliptic filter has an infinite number of Fourier coefficients, which increases the aliasing error in the method. Here, we use a FIR filter with a finitely supported spectrum, and defined by its $N_{H}$ Fourier coefficients of positive index:

The coefficients $h_{k}$ are computed with a Remes algorithm to obtain a filter with a pass-band ripple of $0.02$ and a stop-band rejection of $85\mathrm{dB}$ without zeroes inside the unit circle (Fig. 3). The filter has double transmission zeroes at $\theta=\pm\frac{\pi}{4}$ to guarantee a smooth transition to zero at the edge of the interval. Good results are obtained when using a filter of order $N_{H}=72$. With the $h_{k}$ known, the Fourier coefficients of $\breve Z_{\mathrm{filt}}$ are obtained by convolution:

The $f_{k}^{\mathrm{filt}}$ are then used in $\Phi_{N}$ to compute the unstable poles. With the elliptic filter used in ¹², the sum in the convolution is infinite. With the new filter, this sum becomes finite. We only need Fourier coefficients $f_{k-N_{H}}$ to $f_{k}$ to determine $f_{k}^{\mathrm{filt}}$ exactly. This is the second advantage of the new filter.

Figure 2: Mapping from plane to disc

Figure 3: Filtering function $\breve H\left(\theta\right)$

2.6. Summary

Unstable poles of $Z\!\left(j\omega\right)$ are obtained with the following steps:

1. Transform $Z\!\left(j\omega\right)$ to the unit circle (1)
2. Compute the local rational interpolation models $M_{i}$ (5)
3. Compute the Fourier coefficients of $\breve Z\left(\theta\right)$ (4)
4. Use convolution to apply the filtering function (6)
5. Construct $\Phi_{N}$ and compute its singular values to determine the number of unstable poles
6. Use the PHC algorithm to compute the location of the unstable poles

3. Time-delayed Chua’s circuit

As a first example, we determine the unstable poles of the equilibrium solution of the time-delayed Chua’s circuit shown in Fig. 4¹⁶. The impedance presented by the circuit was determined on 1500 points between 0 and $9\,\mathrm{GHz}$ (Fig. 5). A finer grid of points was placed around the strong resonances in the impedance. The singular values of $\Phi_{50}$ indicate that there are four unstable poles in the circuit (Fig. 5). Using the PHC algorithm, we estimate the following unstable poles:

The method finds the unstable pole pair around $2\mathrm{GHz}$, but it also finds two unstable poles on the real axis. To verify this result, we computed the exact location of the circuit poles using its analytical expression (Appendix 1). The poles obtained with the stability analysis lie very close to the exact poles. The maximum normalised error we obtain is

Figure 4: Time-delayed Chua’s circuit analysed here. The impedance is determined at the location of the red dot.

Figure 5: (left) Impedance presented by the time-delayed Chua’s circuit. (middle) zoom on the impedance on a 10MHz span around 2GHz (right) the first 10 singular values of $\Phi_{50}$.

4. MMIC oscillator

As a more realistic example, we estimate the unstable poles of an MMIC oscillator from the example library of ADS Fig. 6. $Z\!\left(j\omega\right)$ is determined by co-simulation. Lumped models were used for the passives and transistors in the oscillator, while a Momentum simulation was used for the transmission lines. $Z\!\left(j\omega\right)$ is determined at the source of the oscillating transistor on $600$ points between $1\mathrm{MHz}$ and $100\mathrm{GHz}$ (Fig. 7). Because $Z\!\left(j\omega\right)$ is not available close to $0\mathrm{Hz}$, we first use a bandpass transformation on the data. For this circuit, $\Phi_{50}$ has only one significant singular value, which indicates the presence of one unstable pole in the circuit. Its estimated location is at $s=2\pi\left(0.79+j20.94\right)10^{9}$. Here, the exact location of the unstable pole in the equilibrium solution is unknown, so we cannot compute the error. The oscillating solution found with Harmonic Balance has a frequency of $20.42\mathrm{GHz}$, which gives an indication that the estimated unstable pole is correctly estimated.

Figure 6: Lay-out of the MMIC oscillator. The circuit impedance is determined at the location indicated by the arrow.

Figure 7: (left) Impedance of the oscillator and (right) the first five normalised singular values of $\Phi_{50}$

5. Conclusion

In this paper, we have introduced an algorithm to estimate the unstable poles of the equilibrium solution of a microwave circuit. This is no easy task, because the impedance $Z\!\left(j\omega\right)$ of a circuit with distributed elements is not rational. The impedance however only has a finite number of unstable poles: computing its unstable part via a projection step yields a rational function, the unstable poles of which are recovered by the proposed procedure. The algorithm, which computes the singular values of a Hankel matrix based on the Fourier coefficients of the filtered impedance function is shown to yield very accurate pole estimates on two examples.

The presented technique only works to estimate the unstable poles of the equilibrium solution of a circuit. The stable poles remain infinite in number, so we cannot exploit the rationality of the stable part. The same problem is encountered in periodic solutions where both the stable and unstable part have an infinite amount of poles. Extending the pole estimating technique to these cases remains a challenging task.

Appendix 1. Exact poles of the time-delayed Chua’s circuit

We have the following analytical expression for the circuit impedance at the node indicated in red in Fig. 4:

with

To determine the amount of unstable poles, we use the fact that the roots of the denominator quasipolynomial vary continuously with respect to the delay $\tau$ ¹⁷. For $\tau=0$, the denominator is rational and has its root in the left half-plane. Because the circuit is realistic, its denominator will have all roots in the left half-plane for $\tau=0+\varepsilon$ ¹¹. So, to determine the amount of roots in the complex right half-plane at $\tau=0.25\mathrm{ns}$, we compute the amount of zero pairs that cross the imaginary axis in the interval $[0+\varepsilon,\tau]$.

A value $\tau$ for which a zero lies on the imaginary axis is called a critical delay. To compute the critical delays we construct the two variable polynomial $p\left(s,z\right)$ by replacing $z=e^{-s\tau}$ in the denominator quasipolynomial. The solutions of $p\left(s,z\right)=0$ where $s$ lies on the imaginary axis and $z$ on the unit circle are called critical pairs, and can be found by solving a set of non-linear polynomial equations ¹⁷. In this case, $p\left(s,z\right)$ has one critical pair. Its first three corresponding critical delays are $\tau=0.077\mathrm{ps}$, $0.16\mathrm{ns}$ and $0.33\mathrm{ns}$. Two zero pairs have crossed the imaginary axis for our delay of $\tau=0.25\,\mathrm{ns}$, so the denominator of the circuit has 4 unstable zeroes.

We observed two zeroes of the denominator quasipolynomial on the real axis, so we determined their location using dichotomy (fzero in Matlab). The location of the complex pole pair, close to $2\mathrm{GHz}$, was determined using Matlab’s fsolve, with its starting point at $s=j2\pi 2\cdot 10^{9}$.

Acknowledgements

We would like to thank department ELEC of the VUB, who provided the simulation tools for this paper through their strategic research program SRP-19.