3. Stable/unstable projection

The projection described here has been used before to perform stable interpolation and extrapolation of FRF data [21]. With slight modifications, it can be turned into a full-blown stability analysis. We first start with a brief introduction to the notion of Hardy spaces.

The Hardy space $H^{2}(\mathbb{C}^{+})$, is defined as the set of all functions $g$ defined on $\mathbb{C}^{+}$ such that:

• $\forall z\in\mathbb{C}^{+}$, $g$ is holomorphic at $z$
• $\sup_{x>0}\int_{-\infty}^{+\infty}|g(x+j\omega)|^{2}d\omega<\infty$

A classical result [12] [25] states that every function $g\in H^{2}$ admits a limiting function $G(j\omega)$ defined on the imaginary axis $j\mathbb{R}$. The latter is obtained by taking the limit of $g(z)$ when $z$ tends non tangentially toward $j\omega$. Moreover $\forall z\in\mathbb{C}^{+}$, $g(z)$ is equal to the poisson integral of $G$, that is:

The holomorphic nature of $g$ ensures that it is also the Cauchy integral of its boundary value $G$:

There is a one to one linear correspondence between the Hardy functions $g$ and its boundary value function $G$. Using this identification, $H^{2}(\mathbb{C}^{+})$ becomes a subspace of the Hilbert space $L^{2}\!\left(j\mathbb{R}\right)$ of square integrable functions on $j\mathbb{R}$. A direct consequence of (1) is that $H^{2}(\mathbb{C}^{+})$ is closed in $L^{2}\!\left(j\mathbb{R}\right)$ and therefore admits an orthogonal complement. An important result [12] asserts that, $\left(H^{2}(\mathbb{C}^{+})\right)^{\perp}=H^{2}(\mathbb{C}^{-})$ where $H^{2}(\mathbb{C}^{-})$ is defined exactly as $H^{2}(\mathbb{C}^{+})$ above by replacing $\mathbb{C}^{+}$ by $\mathbb{C}^{-}$ and taking the supremum over $x<0$. We therefore have that

This decomposition asserts that any square integrable function on $j\mathbb{R}$ decomposes uniquely as the sum of the traces on the imaginary axis, of an analytic function in the right half-plane and a function analytic in the left half-plane. The projection on $H^{2}\!\left(\mathbb{C}^{+}\right)$ defines the stable part of the function. The projection onto $H^{2}\!\left(\mathbb{C}^{-}\right)$ is the unstable part. As an example, consider $P/Q$ a strictly proper ($\deg(P)<\deg(Q)$) rational function devoid of poles on the imaginary axis. We write its partial fraction expansion as,

where the $\lambda_{i}'s$ with $i\in I^{-}$ are poles belonging to $\mathbb{C}^{-}$, and the ones with $i\in I^{+}$ belong to $\mathbb{C}^{+}$. The strict properness of $P/Q$ ensures its square integrability on $j\mathbb{R}$. By unicity its stable part obtained after projection on $H^{2}\!\left(\mathbb{C}^{+}\right)$ is found to be,

while its unstable part is

In the general case the projection boils down to calculating certain inner products of $Z_{mn}^{\left[ b \right]}\left(j\omega\right)$ with the basis functions $B_{k}$, which form an orthogonal basis of $L^{2}\!\left(j\mathbb{R}\right)$

The overbar $\overline{\bullet}$ indicates the complex conjugate. $\alpha$ is a positive constant used for scaling. All $B_{k}$ with $k\geq0$ create a basis for the stable part, while the $B_{k}$ with negative $k$ form a basis for $H^{2}\!\left(\mathbb{C}^{-}\right)$. Once the $c_{k}$ coefficients are calculated, the stable and unstable parts are easily recovered by calculating

The inner product in (2) runs over all frequencies while the impedance function $Z_{mn}^{\left[ b \right]}\left(j\omega\right)$ is only known over a frequency range $\mathbf{f}=\left[ f_{\mathrm{min}}, f_{\mathrm{max}}\right]$. To impose the finite frequency band on the data, the impedance is filtered before the analysis

The filter $H\left(j\omega\right)$ is a high-order elliptic lowpass filter with its first transmission zero placed at $f_{\mathrm{max}}$ that imposes band limitation. This filter will stabilise poles close to its cutoff frequency, so $f_{\mathrm{max}}$ should be chosen well beyond the maximum frequency at which the circuit can become unstable. $f_{\mathrm{min}}$ should be placed very close to DC. If $f_{\mathrm{min}}$ can’t be close to DC, a bandpass filter should be used for $H\left(j\omega\right)$. The filtering ensures that $Z_{\mathrm{f}}\left(j\omega\right)\in L^{2}\!\left(j\mathbb{R}\right)$ by suppressing anything outside of the frequency band of interest. The smooth decay to zero of $Z_{\mathrm{f}}\left(j\omega\right)$ at the edges of the frequency interval will avoid instabilities to pop up due to the discontinuity of $Z_{mn}^{\left[ b \right]}\left(j\omega\right)$.

In this paper, we use an elliptic filter of order 10 to filter the data. The filter has one transmission zero which is placed exactly at $f_{\mathrm{max}}$.

Figure 3.1 The Möbius transform used in the analysis maps the complex right-half plane to the inside of the unit disc. DC is mapped to $-1$. When $\alpha=2\pi f_{\mathrm{max}}\left(\sqrt{2}-1\right)$, $f_{\mathrm{max}}$ is mapped to $e^{j\frac{\pi}{4}}$.

Working with the basis functions defined in (3) is troublesome from a numeric point of view. Performing the projection when working on unit disc yields better results [21]. The mapping from the complex plane to the unit circle is performed with the Möbius mapping visualised in Figure 3.1:

Our mapping of choice converts square integrable functions on the frequency axis into square integrable functions of the same norm on the unit circle. Square integrable functions which are analytic on the right half-plane are mapped onto square integrable functions which are analytic inside the unit disc. Appendix 1 shows that the basis functions $B_{k}$ map onto powers of $z$

Projecting on this basis boils down to calculating the Fourier series of $Z_{\mathrm{f}}^{\mathrm{disc}}\!\left(z\right)$, with coefficients given by

This Fourier series can be calculated in a numerically efficient way using the Fast Fourier Transform (FFT). The FFT requires that the $\theta$-values are linearly spaced between $0$ and $2\pi$. Due to the mapping from the complex plane to the unit disc, the samples will not satisfy this constraint. A simple interpolation, can be used to obtain $Z_{\mathrm{f}}^{\mathrm{disc}}$ on the $\theta$-values required to perform the FFT.

The interpolation can introduce artefacts in the unstable part if the FRF is not sampled on a sufficiently dense frequency grid, which will be shown on an example later. The level of this interpolation error can be estimated by interpolating the FRF using only the data at the even data points. The difference between the original and the interpolated data for the odd data points will give an indication of the interpolation error encountered in the stability analysis. When the unstable part lies significantly above the interpolation error, we can conclude that the original impedance is unstable.

The threshold to determine when the unstable part lies significantly above the level of the interpolation error will depend on the simulation set-up for the circuit. In the examples that follow, a level of $20 \mathrm{dB}$ has been used as a threshold. A more strict threshold could be used at the cost of requiring a more dense frequency grid and an increased simulation time. When measured components are used in the simulation set-up, the noise level in those measurements should be taken into account to choose the correct threshold.

The stable and unstable parts of a FRF are determined using the following steps:

1. Multiply the FRF with a high-order filter as in (6)
2. Transform the filtered FRF to the unit disc using (7)
3. Interpolate the transformed FRF to a linear grid and use the FFT to calculate the $c_{k}$ coefficients
4. Reconstruct the stable and unstable part using (4)-(5)

A function is meromorphic on $\mathbb{C}$ if it is holomorphic on $\mathbb{C}$ but on a countable number of isolated poles. The function $\tanh(\omega)$ is for example meromorphic, having infinitely many isolated poles on the imaginary axis placed at $j\pi/2+jk\pi\,\,(\forall k\in\mathbb{Z})$ and being analytic elsewhere. In a small-signal stability analysis of a circuit composed of lumped elements, transmission lines and active devices modelled by negative resistors, the impedances can be shown to be meromorphic functions of the frequency. Under the additional realistic assumption that active elements can only deliver power over a finite bandwidth, the impedances are proven to possess only finitely many unstable poles in $\mathbb{C}^{+}$ [5]. Under the generic condition that $Z_{mn}^{\left[ b \right]}\left(j\omega\right)$ is devoid of poles on the imaginary axis, and that the filtering function $H\left(j\omega\right)$ decays strongly enough in order to render $Z_{mn}^{\left[ b \right]}\left(j\omega\right) H\left(j\omega\right)$ square integrable, we conclude that the unstable part of $Z_{mn}^{\left[ b \right]}\left(j\omega\right) H\left(j\omega\right)$ is a rational function. Its poles coincide with the unstable poles of $Z_{mn}^{\left[ b \right]}\left(j\omega\right)$: note here that the multiplication by $H\left(j\omega\right)$ does not add any unstable pole as $H\left(j\omega\right)$ is stable. This means that most of the complexity of the frequency response, like the delay, will be projected onto the stable part, while the unstable part can easily be approximated by a low-order rational model to recover the unstable poles¹.

Classic rational approximation tools can be used to approximate the unstable part and determine the unstable poles. When multiple frequency responses are analysed simultaneously, an approximation method suited for approximation of rational matrices is preferred [20]. In our current implementation, Kung’s method[14] [19] is used to estimate the poles of the unstable part of a single FRF at a time. Alternatively, more sophisticated rational approximation engines like RARL2[20] can be used to recover and track unstable poles.

Compared to working with a high-order rational approximation of the total impedance of the circuit, the “split-first, approximate later” approach proposed here could be a faster and easier method to recover the unstable poles. The post-processing will be faster, but the amount of points required to obtain a sufficiently low interpolation error might be higher than the amount of points required for a tool based on rational approximation.

When the circuit is stable, designers often require information about critical stable poles, to determine how far the circuit is from instability and to track the location of the poles as the circuit varies. In its current form, the projection-based analysis does not simplify finding the location of the stable poles. To perform such an analysis, methods based on local modelling may still be required.