The first step in the pole estimation algorithm is to perform a projection using the PisaProject function. Once the unstable part has been determined, you can estimate the location of the unstable poles with the PisaEstimate function.

Here, we detail the theory and implementation behind the PisaEstimate function. The function estimates the location of the unstable poles in a given frequency function. We assume that you are already familiar with the projection method. If not, please refer first to the projection page.

If you want to see the PisaEstimate function in action, you can check out our list of examples.


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 Latex formula:

Latex formula

In practice, we truncate the infinite Latex formula to Latex formula. Latex formula has yielded good results in all examples we tried. By Kronecker’s theorem for a rational symbol with Latex formula poles, Latex formula, is of rank Latex formula. We compute the rank using the Singular Value Decomposition of Latex formula. Ideally the number of non-zero singular values of Latex formula equals the number of poles of the function Latex formula. In practice due to numerical errors, the first Latex formula singular values are large when compared to the others. This property is used to determine the number of unstable poles Latex formula of Latex formula.

Principal Hankel Components

Once Latex formula has been determined, we can determine the circuit poles using the Principal Hankel Components (PHC) algorithm. Assume we have a state-space realisation of the unstable part of the impendace on the unit circle:

Latex formula

The eigenvalues of the Latex formula-matrix correspond to the poles of the system. The system’s, observability matrix Latex formula and controllability matrix Latex formula are defined as:

Latex formula

the Hankel matrix is equal to the product of Latex formula and Latex formula:

Latex formula

Using the singular value decomposition of the truncated Hankel matrix and the order of the unstable part, we obtain:

Latex formula

such that Latex formula is a diagonal matrix of size Latex formula which contains the singular values, Latex formula is of size Latex formula and Latex formula is size Latex formula. With Latex formula and Latex formula, we compute an estimate of the observability matrix of the approximation:

Latex formula

From the estimated observability matrix Latex formula, we create Latex formula and Latex formula by removing the last and first Latex formula columns respectively. Using the definition of Latex formula we find

Latex formula

We can use fact that Latex formula to estimate the A-matrix by computing the pseudo-inverse of Latex formula:

Latex formula

Latex formula indicates the complex conjugate. The eigenvalues of Latex formula correspond to the location of the unstable poles.

Algorithm steps

Estimating the unstable poles using this PHC method boils down to taking the following steps:

  1. Perform a projection
  2. Construct the Hankel matrix
  3. Compute the singular values of the Hankel matrix
  4. Determine the order from the singular values
  5. Determine the system poles
  6. Transform the poles to the complex plane


The PHC pole estimation algorithm is implemented in the PisaEstimate function:

poles = PisaEstimate(Z);

where Z is an FRM that contains the frequency response of the system. The function returns a vector which contains the estimated pole locations. When the original FRM is specified on the complex plane, the estimated poles will be in the complex right half-plane. The poles are in the laplace variable Latex formula, which is equal to Latex formula or Latex formula where Latex formula is the frequency. When Z is defined on the unit circle, the estimated poles will be inside of the unit circle.

Configuring the Projection

The first step of the PisaEstimate function is to perform a projection. You can pass the parameters for the PisaProject function to PisaEstimate easily. If you want, for example, change the interpolation used in the projection, you can do the following:

poles = PisaEstimate(Z,'Interpolation',RationalInterpolation);

Please refer to the page about projection to read more about the different parameters there.

If you have already performed a projection, you can avoid computing the projection twice by passing the projection results in the Projection parameter:

PROJ = PisaProject(Z);
poles = PisaEstimate(Z,'Projection',PROJ);

Choosing the order

There are two ways control the amount of unstable poles that is estimated in PisaProject. The first is to use the order parameter

poles = PisaEstimate(Z,'order',2);

When you do not specify an order, the PisaEstimate function will plot the singular values of the Hankel matrix and ask you to specify the order on the command line. An example of a singular values plot can be:

We see that the Hankel matrix has 4 significant singular values. This indicates that there are 4 unstable poles in the system. On the matlab command line, you will see the following message:

Enter the amount of unstable poles:

just type 4 and press enter to confirm the order.

Setting the Hankel size

The default number of Fourier coefficients used to construct the Hankel matrix is 100. This size gives good results for low-order systems. When the amount of unstable poles is large, a larger Hankel matrix is recommmended. You can control the number of Fourier coeffients used in the Hankel matrix as follows:

 poles = PisaEstimate(Z,'numberOfCoefficientsForEstimate',200)

Watch out that you are not asking more coefficients than the amount computed in the projection.

Poles estimated outside disc

During the computation of the system’s poles, there is no guarantee that they are estimated inside of the unit circle. Because we are looking for stable poles, they need to be inside the unit circle by definition, so when the function detects poles estimated outside of the unit circle, it will throw an error:

Error using PisaEstimate (line 158)
Estimated a pole outside of the unit circle

When you encounter this error, it has usually one of the following causes:

The order is not chosen correctly. When you ask for an order which is higher than the order of the underlying system, the extra poles will be influenced by the finite interval of the data. This will create poles close to Latex formula. Check the singular values of the Hankel matrix and verify whether you are providing the correct order.

The Hankel matrix is too small. For systems with a high amount of poles inside of the unit circle, a larger Hankel matrix is needed. Try increasing the size of the Hankel matrix to resolve the issue with poles estimated outside the circle.

The interpolation error is high. The influence of the interpolation error on the Fourier coefficients can cause issues with the estimated poles. You can check the interpolation error by setting the “EstimateInterpolationError” parameter in PisaProject to true and by plotting the resulting interpolation error. When it is too large, either change a different interpolation mehtod, or increase the amount of samples in the frequency response.


A. Cooman, F. Seyfert, and S. Amari, “Estimating unstable poles in simulations of microwave circuits” in 2018 ieee mtt-s international microwave symposium (ims), 2018.

Comments are closed.