On this page, we explain the details and implementation of the PisaProject function. This method of functional projection to perform stability analysis is the foundation of the pisa toolbox.

The concept is simple: to determine the stability of a system, project its frequency response on a basis of stable and unstable functions. when a significant part if projected onto the unstable basis functions, we have detected an instability.

Here, we consider the abstract notion of frequency response of a system. In the stability analysis of a circuit, this frequency response is either the impedance or admittance presented by the circuit to a small-signal excitation source. More information about the simulation set-up required to determine such frequency response can be found on the page on closed loop stability analysis.

If you want to see the functional projection technique applied to several examples, go to the examples page.


We are given a frequency response matrix Latex formula of a system with Latex formula outputs and Latex formula inputs on a discrete set of frequencies Latex formula. We assume that Latex formula, which means that we can write Latex formula as the sum of a function in the Hardy space Latex formula (the stable part) and its orthogonal complement Latex formula (the unstable part):

Latex formula

Computing the stable and unstable parts of a given frequency response is done by projecting onto the basis of Latex formula and Latex formula respectively:

Latex formula

In our algorithm, we map the frequency axis to the unit circle such that the complex right half-plane is mapped to the inside of the circle. After this mapping, Latex formula with Latex formula forms a basis for Latex formula and Latex formula with Latex formula forms a basis for Latex formula. The projection is therefore performed by computing the Fourier coefficients of our function on the unit circle.

On the unit circle, computing the projection boils down to computing the Fourier coefficients.

The simulated frequency response is only known for frequencies below Latex formula. Above that frequency, we assume that the frequency response is zero. This sudden jump to zero at Latex formula will have a strong unwanted influence on the obtained stable and unstable part. To suppress this influence, the frequency response is multiplied by a low-pass filter with its stop-band starting at Latex formula.

Algorithm steps

We can summarise the projection algorithm in the following steps:

  1. Perform a frequency normalisation on Latex formula
  2. Transform Latex formula to the unit circle
  3. Multiply by a filtering function Latex formula
  4. Compute the Fourier coefficients Latex formula of Latex formula
  5. Compute the stable and unstable parts by constructing the functions made of positive and negative Fourier coefficients respectively:

Latex formula

  1. Transform Latex formula and Latex formula back to the complex plane
  2. Perform a frequency denormalisation of Latex formula and Latex formula

Frequency normalisation

Internally, the data frequency response data is mapped to the frequency interval Latex formula for ease-of-computation. We provide two different ways to perform this mapping:

In a Bandpass normalisation, we start with data on a frequency interval Latex formula and map Latex formula and Latex formula. The data outside of the interval is assumed to be zero, so using a bandpass normalisation will cause the system to be non-hermitian Latex formula.

When a Lowpass normalisation is used, we map Latex formula and Latex formula. In the lowpass case, the system is assumed to be Hermitian, so the data for the negative frequencies is taken to be the complex conjugate from the data at the positive frequencies. Lowpass normalisation is best used when you know the frequency response for Latex formula. Otherwise, the method will perform an interpolation between Latex formula and Latex formula to estimate the value at zero. This interpolation can introduce a large error in some cases, you should use a bandpass normalisation then.

Mapping to the unit circle

After the frequency normalisation, we are ready to transform the data to the unit circle. To this end, we use Mobius transform.

Latex formula

where Latex formula is chosen such that the following mapping happens in the transform:

Latex formula

The complex right half-plane is mapped to the inside of the unit circle. We can visualise the transform like this:

Filtering to suppress edge effects

After frequency normalisation and mapping to the unit circle, we know the frequency response Latex formula on the interval Latex formula. Outside this interval, we don’t have any data, so we assume it is zero there. This sudden ending of the data will show up in the unstable part, so its influence needs to be suppressed. We do this by multiplying the data by a filtering function Latex formula:

Latex formula

We then compute the projection of Latex formula instead of Latex formula. This filtering will suppress unstable poles close to the interval boundary. It is therefore recommended to simulate the frequency response of the system on a very wide band to minimise the influence of the filtering on the projection result.

The filter we use is a FIR filter:

Latex formula

Because it only consists of positive Fourier coefficients, Latex formula is stable. Its influence will therefore mainly be visible in the stable part of the projection. The only influence the filter has on the unstable poles it that their residues are modified.

The coefficients of the filtering function are precomputed such that the out-of-band suppresion is maximised, for a fixed in-band ripple and transition region. Also, the filtering function is designed to have no zeros inside of the unit disc. That way, they cannot accidentally coincide with an unstable pole. To enforce smoothness in the derivative of the data at the edge of the interval, Two transmission zeros are placed at Latex formula.

The default filter we use in PisaProject has order Latex formula, a transistion region of 10% of the data interval and has an in-band ripple of 0.02:

Computing the Fourier coefficients

We provide two different options to compute the Fourier coefficients of the filtered impedance function. The first computes the classic Fourier integral:

Latex formula

using standard numeric integration techniques. The second option uses the Fast Fourier Transform (FFT) to perform the computation.

The number of Fourier coefficients should be large enough to capture the dynamics in the frequency response function. When a sharp resonance is present in the data, a high amount of Fourier coefficients is needed to capture this resonance.

As a small example, consider the unstable resonance shown below:

On the left, we used too little Fourier coefficients, which shows oscillations in the stable part around the unstable resonance. When we increase the amount of Fourier coefficients (shown on the right), these resonances disappear.

Choosing the correct amount of Fourier coefficients is not so difficult: more is always better. Performing a trade-off between computation time and accuracy is a little bit more difficult and depends on the data.


Both techniques integration techniques (quadrature and FFT) require to interpolate the function values on the unit circle. We either use a simple linear interpolation method, or a rational interpolation method.

Linear interpolation is the easiest to use and to understand. It should be used when the data is noisy because it does not impose continuous derivatives. When the data contains strong resonances, linear interpolation will perform poorly. This is why we also offer a rational interpolation method. The rational interpolation ensures a continuous first derivative, so it can only be used on noise-free data.

In the rational interpolation scheme, we use a rational interpolation function with two zeros and one pole:

Latex formula

where Latex formula, Latex formula,Latex formula and Latex formula are determined using the constraints:

Latex formula

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

Estimating the interpolation error

The interpolation introduces both stable and unstable artefacts into the projection method. To estimate the size of these artefacts, we compute an estimate of the interpolation error. The level of this interpolation error can be estimated by interpolating the frequency response 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.

Interpreting the result

Once the projection has been computed and the impedance function has been split into two parts, we can interpret the results to determine whether the circuit is stable or not. As an example, we show the analysis results for an MMIC power amplifier where a stabilisation resistor is tuned to guarantee circuit stability:

In the plot on the left, the unstable part shows a nice resonance around 10GHz, which indicates an instability. Comparing the size of the obtained unstable part to the level of the interpolation error allows us to affirm that the obtained unstable part is correctly computed and did not arise due to interpolation errors.

In the plot on the right, the unstable part is a lot smaller and more noisy than before. When we compare the unstable part to the interpolation error, we notice that they are both of similar magnitude. This allows us to conclude that the obtained unstable part is due to interpolation errors and not due to an instability in the circuit.

To see more examples on how a result of projection is interpreted, have a look at our different examples.


The projection algorithm described above is implemented in the PisaProject function. Before calling that function, we need to store the circuit impedance data in an FRM object:

Z = FRM(data,'Freq',freq);

where data is an Latex formula array which contains the circuit impedance function(s). and freq is the frequency vector of the data (of length Latex formula).

When the data is gathered in the FRM object, we can pass it to the PisaProject function:

PROJ = PisaProject(Z);

The result is a matlab struct with the following fields:

  • stable: FRM which contains Latex formula
  • unstable: FRM which contains Latex formula
  • coeffs_stable: Positive Fourier coefficients Latex formula (array of size Latex formula)
  • coeffs_unstable: Negative Fourier coefficients Latex formula (array of size Latex formula)
  • InterpErr: FRM which contains the estimated interpolation error made during the projection

Setting the normalisation

The frequency normalisation in pisa is a property of the FRM object.

Z = FRM(data,'Freq',freq,'Normalisation',LowpassNormalisation(freq));
PROJ = PisaProject(Z);

Setting the interpolation

The interpolation scheme used in the computation of the Fourier coefficients is controlled by setting the “Interpolation” parameter in PisaProject:

PROJ = PisaProject(Z,'Interpolation',RationalInterpolation);

Where the provided interpolation is an object which implements the abstract Interpolation class. We implemented two different interpolation methods: RationalInterpolation and LinearInterpolation. The default is linear interpolation, because of its speed and robustness.

Linear interpolation performs poorly when sharp resonances are present in the frequency response, a lot of samples are needed inside of the resonance to prevent the interpolation error to become too big. Linear interpolation should be used when the data is noisy because it does not impose a continuous derivative.

The rational interpolation method is perfectly suited for noiseless data with resonances. It is however slower than linear interpolation and it can encounter numerical issues for data which is constant or linear.

Computation of the Fourier coefficients

PisaProject can compute the Fourier coefficients in two different ways: By Fast Fourier Transform, or by quadrature integration. This is controlled by the “Integration” parameter for PisaProject:

PROJ = PisaProject(Z,'Integration','FFT');

The default is to use the FFT because of its speed. It can however introduce aliasing, but it is the best way to compute a large amount of Fourier coefficients. When only a small amount of Fourier coefficients is needed, it is better to use the quadrature integration.

Besides the computation style, one can control the number of Fourier coefficients by setting the “numberOfCoefficients” parameter in PisaProject:

PROJ = PisaProject(Z,'numberOfCoefficients',1e6);

When the FFT is used to compute the Fourier coefficients, setting this parameter correctly is important to avoid aliasing. When the data contains very sharp resonances, a high amount of Fourier coefficients will be needed to avoid this aliasing.


A. Cooman, F. Seyfert, M. Olivi, S. Chevillard, and L. Baratchart, “Model-free closed-loop stability analysis: a linear functional approach” IEEE transactions on microwave theory and techniques, vol. 66, no. 1, pp. 73-80, 2018.

Comments are closed.