Plot the Bode, Black, Nichols and Nyquist diagrams of a transfer function from an XCOS diagram

[Pascal77]

This code can be used to plot the Bode, Black, Nichols and Nyquist diagrams of a transfer function from an XCOS diagram. It is not necessary to open the XCOS file before running the script.
The post processing “post_xcos_simulate” functionality is used for this operation.

function post_xcos_simulate(%cpr, scs_m, needcompile)
// Do Nothing
endfunction

//We load the XCOS file to be simulated
importXcosDiagram("C:\\...path...\Exemple diagramme de Bode depuis XCOS/Essai_BODE_plot.zcos")


// Clear any pre_xcos_simulation;
clear pre_xcos_simulate;

function post_xcos_simulate(%cpr, scs_m, needcompile)
    // Retrieve all objects
    objs = scs_m.objs;

    clrBlock = [];
    //Looking for CLR block
    for i=1:size(objs)
        if objs(i).gui == "CLR" then
            clrBlock = objs(i);
            break;
        end
    end


    // Retrieve exprs
    s = poly(0,'s');
    expr = clrBlock.graphics.exprs;
    num = evstr(expr(1));
    den = evstr(expr(2));
    h = syslin('c', num/den);
    
    // Open new figure then plot Bode, Black, Nyquist
    clf();//Clear all previous curves
    f = gcf();// Create a figure
    f.figure_size=[1200,1200];
    // f window background color
    f.background=color('lightsteelblue');
    f.figure_name='BODE BLACK NICHOLS NYQUIST DIAGRAMS FROM XCOS FILE'

    subplot(221);
    bode(h, 0.01, 100,"Mon filtre");
   
    subplot(222);
    black(h,0.01,10,"Mon filtre");
    nicholschart(colors=color('gray')*[1 1]);
  
    subplot(223);
    nyquist(h,0.01,10,"Mon filtre");
    hallchart(colors=color('blue')*[1 1])
endfunction

//We run the simulation of the Bode and other diagrams 
xcos_simulate(scs_m, 4);

XCOS diagram is the following one:

The plot is depicted below:

Of course, it is possible to split the code in order to keep only one kind of diagram.

Example:

importXcosDiagram("C:\\...path...\Essai_BODE_plot.zcos")

// Clear any pre_xcos_simulation;
clear pre_xcos_simulate;

function post_xcos_simulate(%cpr, scs_m, needcompile)
    // Retrieve all objects
    objs = scs_m.objs;

    clrBlock = [];
    //Looking for CLR block
    for i=1:size(objs)
        if objs(i).gui == "CLR" then
            clrBlock = objs(i);
            break;
        end
    end


    // Retrieve exprs
    s = poly(0,'s');
    expr = clrBlock.graphics.exprs;
    num = evstr(expr(1));
    den = evstr(expr(2));
    h = syslin('c', num/den);
    // Open new figure then plot Bode
    scf(max(winsid())+1);
    bode(h, 0.01, 100,"Mon filtre");
endfunction

//Bode simulation running
xcos_simulate(scs_m, 4);

[pascal]

actually one plot with Magnitude and Phase as second Y-Axis would be nice, so you have all the information at once.

[Pascal77]

Hello, I have prepared a first try (not perfect) with a common single axis for both plots : Magnitude and Phase.

I used gainplot and phaseplot instructions because it is very convenient but it is also possible to do this with the classical formula 20.log (|h|) and arg(h), but I did not explore this way to do.

Here is the code (I am started from the previous file) :

// Loading the XCOS diagram:
importXcosDiagram("C:\\..path...\Exemple diagramme de Bode depuis XCOS/Essai_BODE_plot.zcos")

// Clear any pre_xcos_simulation;
clear pre_xcos_simulate;

function post_xcos_simulate(%cpr, scs_m, needcompile)
    // Retrieve all objects
    objs = scs_m.objs;

    clrBlock = [];
    //Looking for CLR block
    for i=1:size(objs)
        if objs(i).gui == "CLR" then
            clrBlock = objs(i);
            break;
        end
    end


    // Retrieve exprs
    s = poly(0,'s');
    expr = clrBlock.graphics.exprs;
    num = evstr(expr(1));
    den = evstr(expr(2));
    h = syslin('c', num/den);
    // Open new figure then plot Bode
    //scf(max(winsid())+1);
    //bode(h, 0.01, 100,"Mon filtre");
    //----------------------------------------------------------------------
    //Plot magnitude and phase on 1 plot with 2 y axis
    clf;
    f=gcf();
    f.figure_size=[1000,800];
        a=gca();
    a.background=color('dimgrey');// or any other predefined color name (in color_list)
    a.foreground=color('blue');
    
    gainplot(h,0.01,100);
    phaseplot(h,0.01,100);
    
    //first curve
    poly1= a.children.children(1); //store polyline handle into poly1
    poly1.foreground = 2;  // change the color...
    poly1.thickness = 2;   // ...and the thickness of a curve.

    ylabel(".    Magnitude (dB)    &    Phase (deg)   . ");ylabel edgecolor b background grey95;  
   
    //second curve
    poly2= a.children.children(2); //store polyline handle into poly1
    poly2.foreground = 21;  // change the color...
    poly2.thickness = 2;   // ...and the thickness of a curve.    
    
    legend("Magnitude (dB)","Phase (deg)");
    t = gce();   // get the handle of the newly created object
    t.font_foreground=32; // change font properties
    t.font_size=2;
    t.font_style=9;  

    //title('Bode plot of h(s) = s^2+2*0.9*10*s+100)/(s^2+2*0.3*10.1*s+102.01');//identifying the curve 1 on the plot
    title(["$\frac{1000}{1+20*s+s^2}$"]);
    title font_foreground gold
    title edgecolor red
    title fontsize 4
    title backgroundcolor dimgrey
    title font_style 8
    
    //----------------------------------------------------------------------
endfunction

//Bode simulation running
xcos_simulate(scs_m, 4);

[pascal]

is it possible to do the second Y (Phase on the right side of the graph?)

[Pascal77]

Here is the code to plot magnitude and phase on the same plot with mag axis at letf and phase axis at right.

// plot the magnitude and phase versus frequency
//cf. repfreq help 

clf;
f1=gcf();//
f1.figure_size=[1000,1000];//figure size
f1.figure_position=[100,100];//figure position
f1.background=color('gray90');//color surrounding the plots

//Magnitude and phase computation


/*
f=linspace(fmin,fmax,n); // frequency vector (Hz) (fmin,fmax,number of points)
w=f*2*%pi;//w=pulsation in rad/sec
Pmag=zeros(1,n);Pph=zeros(1,n);//we define the two vectors magnitude and phase (filled with zeros at the start)
for k=1:n, // compute P(jw) at each frequency point
  P=1/(-10*w(k)^2+0.1*%i*w(k)+1); //transfer function in jω (complex expression)
  [Pmag(k),Pph(k)]=dbphi(P);//get magnitude in db and phase in degrees from P (we fill in the two vectors Pmag and Php defined above))
end
*/

s=poly(0,'s');//definition of Laplace variable
h=syslin('c',1/(10*s^2+0.1*s+1));//definition of transfer function with a polynome (rational fraction)

fmin=0.001;//Hz
fmax=10;//Hz
f=fmin:0.001:fmax;//Hz (fmin:step:fmax)
w=f*2*%pi;// conversion in degrees
[fr1,rep] =repfreq(h,f);//definition of frequency response (rep)
[db,phi]=dbphi(rep);//magnitude and phase extraction from rep 


//Plotting curves
//1st curve
ymin = -180;
ymax = 40;
plot2d(f,db,logflag = "ln",rect=[fmin,ymin,fmax,ymax]);  // magnitude plot (same scale than phase plot) rect=[xmin,ymin,xmax,ymax]
xlabel("Frequency (Hz)");
ylabel("Magnitude (dB)");
a = gca();//The left axis is defined (magnitude in dB)
aa = a.children(1).children(1);//this represents the magnitude track itself
aa.foreground=color('red');//we decided to paint it red
aa.thickness=2;//with a thickness of 2
xgrid();//plot the grid
legend("Magnitude (dB)",1);

title("$\frac{1}{10s^2+0.1s+1}$","FontSize",4');//We express the formula for the transfer function in s and jω (in Latex style)
title color blue
title edgecolor blue

//2nd curve
b=newaxes();//We define the second axis (phase in degrees) which will be on the right
plot2d(f,phi,logflag = "ln",rect=[fmin,ymin,fmax,ymax]);  // phase plot (same scale than magnitude plot)
b.filled="off";
b.axes_visible(1)="off";//make y2 axis transparent to avoid masking y1 axis
b.y_location="right";//put y2 axis on the right
ylabel("Phase (degree)");
bb = b.children(1).children(1);//this represents the phase track itself
bb.foreground=color('blue');//we decided to paint it blue
bb.thickness=2;//with a thickness of 2
legend("Phase (°)",4);

Below the code with bode instruction. Very simple compared to the previous one.

s = poly(0, 's');
h1 = syslin('c', 1/(10*s^2+0.1*s+1));


clf(); 
bode(h1, 0.01, 10);

And the corresponding plots: