function definition and help file made coherent + unit test added
Serge Steer [Wed, 16 Jul 2008 09:24:26 +0000 (09:24 +0000)]
scilab/modules/cacsd/help/en_US/calfrq.xml
scilab/modules/cacsd/macros/calfrq.sci
scilab/modules/cacsd/tests/unit_tests/calfrq.dia.ref [new file with mode: 0644]
scilab/modules/cacsd/tests/unit_tests/calfrq.tst [new file with mode: 0644]

index a79de11..4fd514c 100644 (file)
@@ -1,4 +1,4 @@
-<?xml version="1.0" encoding="ISO-8859-1"?>
+<?xml version="1.0" encoding="UTF-8"?>
 <!--
  * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
  * Copyright (C) INRIA - 
  * http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
  *
  -->
-<refentry xmlns="http://docbook.org/ns/docbook" xmlns:xlink="http://www.w3.org/1999/xlink" xmlns:svg="http://www.w3.org/2000/svg" xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:db="http://docbook.org/ns/docbook" version="5.0-subset Scilab" xml:lang="en" xml:id="calfrq">
+<refentry version="5.0-subset Scilab" xml:id="calfrq" xml:lang="en"
+          xmlns="http://docbook.org/ns/docbook"
+          xmlns:xlink="http://www.w3.org/1999/xlink"
+          xmlns:svg="http://www.w3.org/2000/svg"
+          xmlns:ns4="http://www.w3.org/1999/xhtml"
+          xmlns:mml="http://www.w3.org/1998/Math/MathML"
+          xmlns:db="http://docbook.org/ns/docbook">
   <info>
-    <pubdate>$LastChangedDate$</pubdate>
+    <pubdate>$LastChangedDate: 2008-03-26 09:50:39 +0100 (Wed, 26 Mar 2008)
+    $</pubdate>
   </info>
+
   <refnamediv>
     <refname>calfrq</refname>
-    <refpurpose> frequency response discretization</refpurpose>
+
+    <refpurpose>frequency response discretization</refpurpose>
   </refnamediv>
+
   <refsynopsisdiv>
     <title>Calling Sequence</title>
-    <synopsis>[frq,split]=calfrq(h,[fmin,fmax])</synopsis>
+
+    <synopsis>[frq,bnds,split]=calfrq(h,fmin,fmax)</synopsis>
   </refsynopsisdiv>
+
   <refsection>
     <title>Parameters</title>
+
     <variablelist>
       <varlistentry>
         <term>h</term>
+
         <listitem>
-          <para>SISO linear system (<literal>syslin</literal> list)</para>
+          <para>Linear system in state space or transfer representation
+          (<literal>see <link linkend="syslin">syslin</link></literal>)</para>
         </listitem>
       </varlistentry>
+
       <varlistentry>
         <term>fmin,fmax</term>
+
         <listitem>
-          <para>real scalars (min and max frequencies)</para>
+          <para>real scalars (min and max frequencies in Hz)</para>
         </listitem>
       </varlistentry>
+
       <varlistentry>
         <term>frq</term>
+
         <listitem>
-          <para>row vector (discretization of interval)</para>
+          <para>row vector (discretization of the frequency interval)</para>
         </listitem>
       </varlistentry>
+
+      <varlistentry>
+        <term>bnds</term>
+
+        <listitem>
+          <para>vector <literal>[Rmin Rmax Imin Imax]</literal> where
+          <literal>Rmin</literal> and <literal>Rmax</literal> are the lower
+          and upper bounds of the frequency response real part,
+          <literal>Imin</literal> and <literal>Imax</literal> are the lower
+          and upper bounds of the frequency response imaginary part,</para>
+        </listitem>
+      </varlistentry>
+
       <varlistentry>
         <term>split</term>
+
         <listitem>
           <para>vector of frq splitting points indexes</para>
         </listitem>
       </varlistentry>
     </variablelist>
   </refsection>
+
   <refsection>
     <title>Description</title>
-    <para>
-    frequency response discretization ; <literal>frq</literal> is the discretization of  
-    <literal>[fmin,fmax]</literal> 
-    such that the peaks in the frequency response are well represented.</para>
-    <para>
-    Default values for <literal>fmin</literal> and <literal>fmax</literal> are
-    <literal>1.d-3</literal>, <literal>1.d+3</literal> if <literal>h</literal> is continuous-time
-    or <literal>1.d-3</literal>, <literal>1/(2*h('dt'))</literal> if <literal>h</literal> is discrete-time.</para>
-    <para>
-    Singularities are located between <literal>frq(split(k))</literal> and <literal>frq(split(k)+1)</literal>
-    for <literal>k&gt;1</literal>.</para>
+
+    <para>frequency response discretization; <literal>frq</literal> is the
+    discretization of <literal>[fmin,fmax]</literal> such that the peaks in
+    the frequency response are well represented.</para>
+
+    <para>Singularities are located between <literal>frq(split(k)-1)</literal>
+    and <literal>frq(split(k))</literal> for <literal>k&gt;1</literal>.</para>
   </refsection>
+
   <refsection>
     <title>Examples</title>
-    <programlisting role="example"><![CDATA[
+
+    <programlisting role="example">
 
 s=poly(0,'s')
 h=syslin('c',(s^2+2*0.9*10*s+100)/(s^2+2*0.3*10.1*s+102.01))
 h1=h*syslin('c',(s^2+2*0.1*15.1*s+228.01)/(s^2+2*0.9*15*s+225)) 
-[f1,spl]=calfrq(h1,0.01,1000);
+[f1,bnds,spl]=calfrq(h1,0.01,1000);
 rf=repfreq(h1,f1);
 plot2d(real(rf)',imag(rf)')
  
-  ]]></programlisting>
+  </programlisting>
   </refsection>
+
   <refsection>
     <title>See Also</title>
+
     <simplelist type="inline">
-      <member>
-        <link linkend="bode">bode</link>
-      </member>
-      <member>
-        <link linkend="black">black</link>
-      </member>
-      <member>
-        <link linkend="nyquist">nyquist</link>
-      </member>
-      <member>
-        <link linkend="freq">freq</link>
-      </member>
-      <member>
-        <link linkend="repfreq">repfreq</link>
-      </member>
-      <member>
-        <link linkend="logspace">logspace</link>
-      </member>
+      <member><link linkend="bode">bode</link></member>
+
+      <member><link linkend="black">black</link></member>
+
+      <member><link linkend="nyquist">nyquist</link></member>
+
+      <member><link linkend="freq">freq</link></member>
+
+      <member><link linkend="repfreq">repfreq</link></member>
+
+      <member><link linkend="logspace">logspace</link></member>
     </simplelist>
   </refsection>
-</refentry>
+</refentry>
\ No newline at end of file
index 7d6ada4..736ce93 100644 (file)
@@ -14,7 +14,7 @@ function [frq,bnds,splitf]=calfrq(h,fmin,fmax)
   eps=1.d-14   //minimum absolute lower frequency
   k=0.001;     // Minimum relative prediction error in the nyquist plan 
   epss=0.002   // minimum frequency distance with a singularity
-  nptmax=5000 //maximum number of discretisation points
+  nptmax=5000  //maximum number of discretisation points
   tol=0.01     // Tolerance for testing pure imaginary numbers 
 
   // Check inputs
@@ -22,14 +22,11 @@ function [frq,bnds,splitf]=calfrq(h,fmin,fmax)
   if and(typeof(h)<>['state-space' 'rational'])
    error(msprintf(gettext("%s: Wrong type for input argument #%d: Linear state space or a transfer function expected.\n"),"calfrq",1))
   end
-  flag=h(1);
-  if flag(1)=='lss' then
+  if typeof(h)=='state-space' then
     h=ss2tf(h)
   end
+  
   [m,n]=size(h.num)
-  if  n<>1 then  
-    error(msprintf(gettext("%s: Wrong size for input argument #%d: Single input system expected.\n"),"calfrq",1))
-  end
   dom=h('dt')
   select dom
   case 'd' then 
diff --git a/scilab/modules/cacsd/tests/unit_tests/calfrq.dia.ref b/scilab/modules/cacsd/tests/unit_tests/calfrq.dia.ref
new file mode 100644 (file)
index 0000000..b1d1c69
--- /dev/null
@@ -0,0 +1,85 @@
+// =============================================================================
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) ????-2008 - rINRIA - Serge Steer
+//
+//  This file is distributed under the same license as the Scilab package.
+// =============================================================================
+K=0.001;     // Minimum relative prediction error in the nyquist plan
+Epss=0.002;   // minimum frequency distance with a singularity
+nptmax=5000;  //maximum number of discretisation points
+pas=100/(2*%pi);
+s=%s;
+h=syslin('c',1/%s);n=1;
+[f,bnds,split]=calfrq(h,0.01,100);
+if split<>1 then bugmes();quit;end
+rf=freq(h.num,h.den,2*%pi*%i*f);
+//finite difference derivative estimate
+ddf=diff(f)/pas;
+drf=(freq(h.num,h.den,2*%pi*%i*(f(1:$-1)+ddf))-rf(:,1:$-1));
+//error between computed and extrapolated value
+e=rf(:,2:$)-(rf(:,1:$-1)+drf*pas);
+if max([abs(real(e))/max(bnds(2)-bnds(1),1); abs(imag(e))/max(bnds(4)- bnds(3),1)])>K then bugmes();quit;end
+h=syslin('c',[1/(%s+0.5);1/(%s+0.3)]);n=2;
+[f,bnds,split]=calfrq(h,0.01,100);
+if split<>1 then bugmes();quit;end
+rf=freq(h.num,h.den,2*%pi*%i*f);
+//finite difference derivative estimate
+ddf=diff(f)/pas;
+drf=(freq(h.num,h.den,2*%pi*%i*(f(1:$-1)+ddf))-rf(:,1:$-1));
+//error between computed and extrapolated value
+e=rf(:,2:$)-(rf(:,1:$-1)+drf*pas);
+if max([abs(real(e))/max(bnds(2)-bnds(1),1); abs(imag(e))/max(bnds(4)- bnds(3),1)])>K then bugmes();quit;end
+h=syslin('c',(%s^2+2*0.9*10*%s+100)/(%s^2+2*0.3*10.1*%s+102.01));
+//h=h*syslin('c',(%s)/(%s^2+81)) ;
+n=1;
+[f,bnds,split]=calfrq(h,0.01,100);
+if size(split,1)<>1 then bugmes();quit;end
+rf=freq(h.num,h.den,2*%pi*%i*f);
+//finite difference derivative estimate
+ddf=diff(f)/100;
+drf=(freq(h.num,h.den,2*%pi*%i*(f(1:$-1)+ddf))-rf(:,1:$-1));
+//error between computed and extrapolated value
+e=rf(:,2:$)-(rf(:,1:$-1)+drf*100);
+if max([abs(real(e))/max(bnds(2)-bnds(1),1); abs(imag(e))/max(bnds(4)- bnds(3),1)])>K then bugmes();quit;end
+//case with singularity
+h=syslin('c',(%s)/(%s^2+81)) ;
+sing=9/(2*%pi);
+n=1;
+[f,bnds,split]=calfrq(h,0.01,100);
+if size(split,2)<>2 then bugmes();quit;end
+ks=split(2);
+if abs(f(ks-1)-f(ks))*2*%pi<Epss then bugmes();quit;end
+if f(ks-1)>sing|f(ks)<sing then bugmes();quit;end
+//finite difference derivative estimate
+f1=f(1:ks-15);//remove points near singularity for which pasmin constraint may be active
+rf=freq(h.num,h.den,2*%pi*%i*f1);
+ddf=diff(f1)/pas;
+drf=(freq(h.num,h.den,2*%pi*%i*(f1(1:$-1)+ddf))-rf(:,1:$-1));
+//error between computed and extrapolated value
+e=rf(:,2:$)-(rf(:,1:$-1)+drf*pas);
+if max([abs(real(e))/max(bnds(2)-bnds(1),1); abs(imag(e))/max(bnds(4)- bnds(3),1)])>K then bugmes();quit;end
+f1=f(ks+15:$);//remove points near singularity for which pasmin constraint may be active
+rf=freq(h.num,h.den,2*%pi*%i*f1);
+ddf=diff(f1)/pas;
+drf=(freq(h.num,h.den,2*%pi*%i*(f1(1:$-1)+ddf))-rf(:,1:$-1));
+//error between computed and extrapolated value
+e=rf(:,2:$)-(rf(:,1:$-1)+drf*pas);
+if max([abs(real(e))/max(bnds(2)-bnds(1),1); abs(imag(e))/max(bnds(4)- bnds(3),1)])>K then bugmes();quit;end
+//state space case
+sl=tf2ss(h);
+n=1;
+[f1,bnds1,split1]=calfrq(sl,0.01,100);
+if norm(f-f1)>1d-13 then bugmes();quit;end
+if or(split<>split1) then bugmes();quit;end
+//discrete case
+h=syslin('c',(s^2+2*0.9*10*s+100)/(s^2+2*0.3*10.1*s+102.01));
+hd=ss2tf(dscr(h,0.01));
+[f,bnds,split]=calfrq(hd,0.00001,3);
+if size(split,1)<>1 then bugmes();quit;end
+rf=freq(h.num,h.den,exp(2*%pi*%i*f));
+//finite difference derivative estimate
+ddf=diff(f)/pas;
+drf=(freq(h.num,h.den,exp(2*%pi*%i*(f(1:$-1)+ddf)))-rf(:,1:$-1));
+//error between computed and extrapolated value
+e=rf(:,2:$)-(rf(:,1:$-1)+drf*pas);
+if max([abs(real(e))/max(bnds(2)-bnds(1),1); abs(imag(e))/max(bnds(4)- bnds(3),1)])>K then bugmes();quit;end
diff --git a/scilab/modules/cacsd/tests/unit_tests/calfrq.tst b/scilab/modules/cacsd/tests/unit_tests/calfrq.tst
new file mode 100644 (file)
index 0000000..469bb94
--- /dev/null
@@ -0,0 +1,104 @@
+// =============================================================================
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) ????-2008 - rINRIA - Serge Steer
+//
+//  This file is distributed under the same license as the Scilab package.
+// =============================================================================
+K=0.001;     // Minimum relative prediction error in the nyquist plan 
+Epss=0.002;   // minimum frequency distance with a singularity
+nptmax=5000;  //maximum number of discretisation points
+pas=100/(2*%pi);
+s=%s;
+
+h=syslin('c',1/%s);n=1;
+
+
+[f,bnds,split]=calfrq(h,0.01,100);
+if split<>1 then pause,end
+rf=freq(h.num,h.den,2*%pi*%i*f);
+//finite difference derivative estimate
+ddf=diff(f)/pas;
+drf=(freq(h.num,h.den,2*%pi*%i*(f(1:$-1)+ddf))-rf(:,1:$-1));
+//error between computed and extrapolated value
+e=rf(:,2:$)-(rf(:,1:$-1)+drf*pas);
+if max([abs(real(e))/max(bnds(2)-bnds(1),1); abs(imag(e))/max(bnds(4)- bnds(3),1)])>K then pause,end
+
+h=syslin('c',[1/(%s+0.5);1/(%s+0.3)]);n=2;
+[f,bnds,split]=calfrq(h,0.01,100);
+if split<>1 then pause,end
+rf=freq(h.num,h.den,2*%pi*%i*f);
+//finite difference derivative estimate
+ddf=diff(f)/pas;
+drf=(freq(h.num,h.den,2*%pi*%i*(f(1:$-1)+ddf))-rf(:,1:$-1));
+//error between computed and extrapolated value
+e=rf(:,2:$)-(rf(:,1:$-1)+drf*pas);
+if max([abs(real(e))/max(bnds(2)-bnds(1),1); abs(imag(e))/max(bnds(4)- bnds(3),1)])>K then pause,end
+
+
+
+h=syslin('c',(%s^2+2*0.9*10*%s+100)/(%s^2+2*0.3*10.1*%s+102.01));
+//h=h*syslin('c',(%s)/(%s^2+81)) ;
+n=1;
+[f,bnds,split]=calfrq(h,0.01,100);
+if size(split,1)<>1 then pause,end
+
+rf=freq(h.num,h.den,2*%pi*%i*f);
+//finite difference derivative estimate
+ddf=diff(f)/100;
+drf=(freq(h.num,h.den,2*%pi*%i*(f(1:$-1)+ddf))-rf(:,1:$-1));
+//error between computed and extrapolated value
+e=rf(:,2:$)-(rf(:,1:$-1)+drf*100);
+if max([abs(real(e))/max(bnds(2)-bnds(1),1); abs(imag(e))/max(bnds(4)- bnds(3),1)])>K then pause,end
+
+
+
+
+
+
+//case with singularity
+h=syslin('c',(%s)/(%s^2+81)) ;
+sing=9/(2*%pi);
+n=1;
+[f,bnds,split]=calfrq(h,0.01,100);
+if size(split,2)<>2 then pause,end
+ks=split(2);
+if abs(f(ks-1)-f(ks))*2*%pi<Epss then pause,end
+if f(ks-1)>sing|f(ks)<sing then pause,end
+
+//finite difference derivative estimate
+f1=f(1:ks-15);//remove points near singularity for which pasmin constraint may be active
+rf=freq(h.num,h.den,2*%pi*%i*f1);
+ddf=diff(f1)/pas;
+drf=(freq(h.num,h.den,2*%pi*%i*(f1(1:$-1)+ddf))-rf(:,1:$-1));
+//error between computed and extrapolated value
+e=rf(:,2:$)-(rf(:,1:$-1)+drf*pas);
+if max([abs(real(e))/max(bnds(2)-bnds(1),1); abs(imag(e))/max(bnds(4)- bnds(3),1)])>K then pause,end
+f1=f(ks+15:$);//remove points near singularity for which pasmin constraint may be active
+rf=freq(h.num,h.den,2*%pi*%i*f1);
+ddf=diff(f1)/pas;
+drf=(freq(h.num,h.den,2*%pi*%i*(f1(1:$-1)+ddf))-rf(:,1:$-1));
+//error between computed and extrapolated value
+e=rf(:,2:$)-(rf(:,1:$-1)+drf*pas);
+if max([abs(real(e))/max(bnds(2)-bnds(1),1); abs(imag(e))/max(bnds(4)- bnds(3),1)])>K then pause,end
+
+
+//state space case
+sl=tf2ss(h);
+n=1;
+[f1,bnds1,split1]=calfrq(sl,0.01,100);
+if norm(f-f1)>1d-13 then pause,end
+if or(split<>split1) then pause,end
+
+//discrete case
+h=syslin('c',(s^2+2*0.9*10*s+100)/(s^2+2*0.3*10.1*s+102.01));
+hd=ss2tf(dscr(h,0.01));
+[f,bnds,split]=calfrq(hd,0.00001,3);
+if size(split,1)<>1 then pause,end
+
+rf=freq(h.num,h.den,exp(2*%pi*%i*f));
+//finite difference derivative estimate
+ddf=diff(f)/pas;
+drf=(freq(h.num,h.den,exp(2*%pi*%i*(f(1:$-1)+ddf)))-rf(:,1:$-1));
+//error between computed and extrapolated value
+e=rf(:,2:$)-(rf(:,1:$-1)+drf*pas);
+if max([abs(real(e))/max(bnds(2)-bnds(1),1); abs(imag(e))/max(bnds(4)- bnds(3),1)])>K then pause,end