Zero-crossings of Nyquist Diagram

Hi,

I have been trying to extract all the zero-crossings of the Nyquist plot of my transfer function zp. So I want to obtain the real part value where the Nyquist plot crosses the Real Axis.

Here is an example of how the Nyquist diagram could look like (I keep changing specific parameters so the Nyquist diagram also keeps changing):

I already tried extracting the zero crossings with the folloing code,

[re,im] = nyquist(zp);
Re = squeeze(re);
Im = squeeze(im);

With the given values I searched where the imaginary part was close to zero. But unfortunately it isn't precise enough and what I need is very detailed data. Increasing the amount of values also didn't get me much further. This is the code I used to increase my amount of values:

upperbound = 2*pi*(1600);
lowerbound = 2*pi*(-1600);

increase_number_of_Values = linspace(lowerbound,upperbound,5000);
[re,im] = nyquist(zp, increase_number_of_Values);
Re = squeeze(re);
Im = squeeze(im);

This code doesn't give me every zero-crossing of the Nyquist plot. I already changed the frequency interval various times, yet it doesn't show me every zero-crossing. Then I tried to extract the zero crossings directly from the transfer function itself. With the following code:

fp = [-1600:0.1:1600];
zp_values = [];
for k = 1 : length(fp)
    w= 2*pi*fp(k);
     zp_values(k) = (my transfer function); %it calculates the value of zp for every frequency in my array fp
end  

However this also did not work. I chose the frequency interval -1600Hz <= fp <= 1600Hz because the values which I want to observe in the Nyquist plot have a frequency between -10000 rad/s to 10000 rad/s.

Has anybody an idea how I can extract the zero-crossings of my Nyquist plot, as accurately as possible?

Thanks very much for your help!

Answer 1

[Kiran Felix Robert]

Hi Leopold,

I understand that you need to find the zero-crossings of the imaginary part of the Nyquist diagram. I assume that you have already extracted the real-part and the imaginary part separately.

In that case you can refer the following answer thread for a solution,

https://in.mathworks.com/matlabcentral/answers/267222-easy-way-of-finding-zero-crossing-of-a-function

The following is a specific example

H = tf([2 5 1],[1 2 3]);
[re,im] = nyquist(H);
re = squeeze(re);
im = squeeze(im);
plot(re,im)
zci = @(v) find(v(:).*circshift(v(:), [-1 0]) <= 0); 
zind = zci(im);
figure(1)
plot(re,im,'r');
hold on
plot(re(zind),im(zind),'bp')
hold off
grid

You can also try interpolating your data for a more accurate values.

[Leo]

Hi Kiran,

thanks for the quick answer, it helped a lot. I have one more problem. The code you gave tells me where the zero-crossings of my Nyquist plot are. But it doesn't tell me how often this exact spot is being crossed, because it can be more than one time.

In the following plot you can see, that the Real-Axis gets crossed two times at Re= -11.752. However the given code tells me that there is only one crossing. For what I am trying to do it is also important how often the Real-Axis gets crossed.

I tried to solve it by comparing the frequencies of both crossings (at the same zero-crossing spot) because they have opposite signs, aus you can see in the second copy of the plot.

However exctracting the frequency from the Nyquist plot with the following code I only received frequencies (in rad/s) with a positive sign.

[re,im,w] = nyquist(zp);
Re = squeeze(re);
Im = squeeze(im);

I believe that is the reason why the crossings of the Real-Axis at the point Re=-11.752 are only counted once, because the Matlab function "[re,im,w]= nyquist(zp)" only checks on the positive frequencies. The reason why it only extracts the values for the positive frequencies could be the symmetry of my plot.

Do you know how I can also extract the values for the negative frequencies?

Or do you have another idea for how I can count the number of crossings at the same crossing point of the Real-Axis?

Thanks very much for your help!

[Paul]

Wouldn't it be true by symmetry that if the Nyquist plot crosses the real axis at frequency w it must also cross the real axis (at the same point) at -w, as shown in your plot?

Answer 2

[Paul]

By definition, the gain margin is defined at the point(s) where the Nyquist plot crosses the real axis. Use the allmargin() function to find the gain margin(s). The actual crossing point is at 1/GM.

>> h=tf([2 5 1],[1 2 3])*tf(1,[1 2 2 1])  % example transfer function

h =
 
             2 s^2 + 5 s + 1
  --------------------------------------
  s^5 + 4 s^4 + 9 s^3 + 11 s^2 + 8 s + 3
 
Continuous-time transfer function.

>> nyquist(h)

>> s=allmargin(h)

s = 

  struct with fields:

     GainMargin: 2.250221178446187e+00
    GMFrequency: 1.870895057257997e+00
    PhaseMargin: [1.675188645227028e+02 5.814149251422266e+01]
    PMFrequency: [5.894485378547998e-01 1.287131228874161e+00]
    DelayMargin: [4.960154377113239e+00 7.883892905655167e-01]
    DMFrequency: [5.894485378547998e-01 1.287131228874161e+00]
         Stable: 1

>> realcrossing1 = -1./s.GainMargin  % crossings on the negative real axis

realcrossing1 =

    -4.444007591691567e-01

>> s=allmargin(-h)

s = 

  struct with fields:

     GainMargin: [3.000000000000005e+00 1.186263434433761e+00]
    GMFrequency: [0 4.858495157107252e-01]
    PhaseMargin: [-1.248113547729714e+01 -1.218585074857773e+02]
    PMFrequency: [5.894485378547998e-01 1.287131228874161e+00]
    DelayMargin: [1.028986927474128e+01 3.229160350356917e+00]
    DMFrequency: [5.894485378547998e-01 1.287131228874161e+00]
         Stable: 1

>> realcrossing2 = 1./s.GainMargin  % crossings on the postive real axis

realcrossing2 =

     3.333333333333328e-01     8.429830769228168e-01

As shown, these crossings are for positive (and zero) frequencies. By symmetry, the same crossings should occur for negative frequencies (don't count w = 0 or w=inf twice). Presumably you don't care about crossings on the real axis at infinity due to poles on the imaginary axis.