Move lattp as a dedicated macros (it was documented as a single function) 39/11039/2
Sylvestre Ledru [Wed, 27 Mar 2013 09:40:09 +0000 (10:40 +0100)]
Change-Id: Ib8b6eb14d82c97c3741beb2e91f3f08077ef00c3

scilab/modules/signal_processing/macros/lattn.sci
scilab/modules/signal_processing/macros/lattp.sci [new file with mode: 0644]

index 37e9b3e..627e1b2 100644 (file)
-// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab\r
-// Copyright (C) INRIA - 1989 - G. Le Vey\r
-// Copyright (C) 2012 - Scilab Enterprises - Adeline CARNIS\r
-// \r
-// This file must be used under the terms of the CeCILL.\r
-// This source file is licensed as described in the file COPYING, which\r
-// you should have received as part of this distribution.  The terms\r
-// are also available at    \r
-// http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt\r
-function [la,lb]=lattn(n,p,cov)\r
-//[la,lb]=lattn(n,p,cov)\r
-//macro which solves recursively on n (p being fixed)\r
-//the following system (normal equations), i.e. identifies\r
-//the AR part(poles) of a vector ARMA(n,p) process\r
-//\r
-//                        |Rp+1 Rp+2 . . . . . Rp+n  |\r
-//                        |Rp   Rp+1 . . . . . Rp+n-1|\r
-//                        | .   Rp   . . . . .  .    |\r
-//                        |                          |\r
-//    |I -A1 -A2 . . .-An|| .    .   . . . . .  .    |=0\r
-//                        | .    .   . . . . .  .    |\r
-//                        | .    .   . . . . .  .    |\r
-//                        | .    .   . . . . . Rp+1  |\r
-//                        |Rp+1-n.   . . . . . Rp    |\r
-//\r
-//    where {Rk;k=1,nlag} is the sequence of empirical covariances\r
-//\r
-//   n   : maximum order of the filter\r
-//   p   : fixed dimension of the MA part. If p is equal to -1,\r
-//       : the algorithm reduces to the classical Levinson recursions.\r
-//   cov : matrix containing the Rk(d*d matrices for\r
-//       : a d-dimensional process).It must be given the\r
-//       : following way:\r
-//\r
-//                        |  R0 |\r
-//                        |  R1 |\r
-//                    cov=|  R2 |\r
-//                        |  .  |\r
-//                        |  .  |\r
-//                        |Rnlag|\r
-//\r
-//   la  : list-type variable, giving the successively calculated\r
-//       : polynomials (degree 1 to degree n),with coefficients Ak\r
-//!\r
-\r
-\r
-   if argn(2)<>3 then\r
-     error(msprintf(gettext("%s: Wrong number of input argument(s): %d expected.\n"),"lattn",3))\r
-   end\r
-   if type(n)<>1 then\r
-     error(msprintf(gettext("%s: Wrong type for input argument #%d: Array of floating point numbers expected.\n"),"lattn",1))\r
-   end\r
-   if  size(n,'*')<>1 then\r
-     error(msprintf(gettext("%s: Wrong size for input argument #%d: A scalar expected.\n"),"lattn",1))\r
-   end\r
-   if  n<0|n<>round(n) then\r
-     error(msprintf(gettext("%s: Wrong values for input argument #%d: Non-negative integers expected.\n"),"lattn",1))\r
-   end\r
-   if type(p)<>1 then\r
-     error(msprintf(gettext("%s: Wrong type for input argument #%d: Array of floating point numbers expected.\n"),"lattn",2))\r
-   end\r
-   if  size(p,'*')<>1 then\r
-     error(msprintf(gettext("%s: Wrong size for input argument #%d: A scalar expected.\n"),"lattn",2))\r
-   end\r
-   if  p<-1|p<>round(p) then\r
-     error(msprintf(gettext("%s: Wrong values for input argument #%d: Elements must be in the interval [%s, %s].\n"),"lattn",2,"-1","%inf"))\r
-   end\r
-   \r
-   if type(cov)<>1 then\r
-     error(msprintf(gettext("%s: Wrong type for input argument #%d: Array of floating point numbers expected.\n"),"lattn",3))\r
-   end\r
-\r
-   [l,d]=size(cov);\r
-   if d>l then\r
-     error(msprintf(gettext("%s: Wrong size for input argument #%d: A tall matrix expected.\n"),"lattn",3))\r
-   end\r
-   \r
-   \r
-   a=eye(d);b=eye(d);\r
-   z=poly(0,'z');la=list();lb=list();\r
-   no=p-n-1;cv=cov;\r
-   \r
-   if -no >= floor(l/d) then\r
-       error(msprintf(gettext("%s: Wrong values for input arguments #%d and #%d.\n"),"lattn",1, 2))\r
-   end\r
-   \r
-   if no<0,\r
-     for j=1:-no,cv=[cov(j*d+1:(j+1)*d,:)';cv];end;\r
-     p=p-no;\r
-   end\r
-   \r
-   if (p + 2 + (n-1)) * d >  size(cv,'r') then\r
-       error(msprintf(gettext("%s: Wrong values for input arguments #%d and #%d.\n"),"lattn",1, 2))\r
-   end\r
-   \r
-   for j=0:n-1,\r
-     jd=jmat(j+1,d);\r
-     r1=jd*cv((p+1)*d+1:(p+2+j)*d,:);\r
-     r2=jd*cv(p*d+1:(p+1+j)*d,:);\r
-     r3=jd*cv((p-1-j)*d+1:p*d,:);\r
-     r4=jd*cv((p-j)*d+1:(p+1)*d,:);\r
-     c1=coeff(a);c2=coeff(b);\r
-     k1=(c1*r1)*inv(c2*r2);\r
-     k2=(c2*r3)*inv(c1*r4);\r
-     a1=a-k1*z*b;\r
-     b=-k2*a+z*b;a=a1;\r
-     la(j+1)=a;lb(j+1)=b;\r
-   end;\r
-endfunction\r
-function [la,lb]=lattp(n,p,cov)\r
-// Author G Levey INRIA\r
-//!\r
-  [l,d]=size(cov);\r
-  id=eye(d);\r
-  [a,b]=lattn(n,0,cov);a=a(n);b=b(n);\r
-  z=poly(0,'z');la=list();lb=list();\r
-  jd=jmat(n+1,d);\r
-  for j=0:p-1,\r
-    r1=jd*cov((j+1)*d+1:(j+n+2)*d,:);\r
-    r2=jd*cov(j*d+1:(j+n+1)*d,:);\r
-    c1=coeff(a);c2=coeff(b);\r
-    k=(c1*r1)*inv(c2*r2);\r
-    hst=-inv(c1(:,n*d+1:(n+1)*d));\r
-    r=k*hst;\r
-    a=(id-r*z)*a-k*z*b;\r
-    b=-hst*a;\r
-    la(j+1)=a;\r
-  end;\r
-endfunction\r
-\r
-\r
-\r
-\r
-\r\r
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) INRIA - 1989 - G. Le Vey
+// Copyright (C) 2012 - Scilab Enterprises - Adeline CARNIS
+//
+// This file must be used under the terms of the CeCILL.
+// This source file is licensed as described in the file COPYING, which
+// you should have received as part of this distribution.  The terms
+// are also available at
+// http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
+
+function [la,lb]=lattn(n,p,cov)
+    //[la,lb]=lattn(n,p,cov)
+    //macro which solves recursively on n (p being fixed)
+    //the following system (normal equations), i.e. identifies
+    //the AR part(poles) of a vector ARMA(n,p) process
+    //
+    //                        |Rp+1 Rp+2 . . . . . Rp+n  |
+    //                        |Rp   Rp+1 . . . . . Rp+n-1|
+    //                        | .   Rp   . . . . .  .    |
+    //                        |                          |
+    //    |I -A1 -A2 . . .-An|| .    .   . . . . .  .    |=0
+    //                        | .    .   . . . . .  .    |
+    //                        | .    .   . . . . .  .    |
+    //                        | .    .   . . . . . Rp+1  |
+    //                        |Rp+1-n.   . . . . . Rp    |
+    //
+    //    where {Rk;k=1,nlag} is the sequence of empirical covariances
+    //
+    //   n   : maximum order of the filter
+    //   p   : fixed dimension of the MA part. If p is equal to -1,
+    //       : the algorithm reduces to the classical Levinson recursions.
+    //   cov : matrix containing the Rk(d*d matrices for
+    //       : a d-dimensional process).It must be given the
+    //       : following way:
+    //
+    //                        |  R0 |
+    //                        |  R1 |
+    //                    cov=|  R2 |
+    //                        |  .  |
+    //                        |  .  |
+    //                        |Rnlag|
+    //
+    //   la  : list-type variable, giving the successively calculated
+    //       : polynomials (degree 1 to degree n),with coefficients Ak
+    //!
+
+
+    if argn(2)<>3 then
+        error(msprintf(gettext("%s: Wrong number of input argument(s): %d expected.\n"),"lattn",3))
+    end
+    if type(n)<>1 then
+        error(msprintf(gettext("%s: Wrong type for input argument #%d: Array of floating point numbers expected.\n"),"lattn",1))
+    end
+    if  size(n,"*")<>1 then
+        error(msprintf(gettext("%s: Wrong size for input argument #%d: A scalar expected.\n"),"lattn",1))
+    end
+    if  n<0|n<>round(n) then
+        error(msprintf(gettext("%s: Wrong values for input argument #%d: Non-negative integers expected.\n"),"lattn",1))
+    end
+    if type(p)<>1 then
+        error(msprintf(gettext("%s: Wrong type for input argument #%d: Array of floating point numbers expected.\n"),"lattn",2))
+    end
+    if  size(p,"*")<>1 then
+        error(msprintf(gettext("%s: Wrong size for input argument #%d: A scalar expected.\n"),"lattn",2))
+    end
+    if  p<-1|p<>round(p) then
+        error(msprintf(gettext("%s: Wrong values for input argument #%d: Elements must be in the interval [%s, %s].\n"),"lattn",2,"-1","%inf"))
+    end
+
+    if type(cov)<>1 then
+        error(msprintf(gettext("%s: Wrong type for input argument #%d: Array of floating point numbers expected.\n"),"lattn",3))
+    end
+
+    [l,d]=size(cov);
+    if d>l then
+        error(msprintf(gettext("%s: Wrong size for input argument #%d: A tall matrix expected.\n"),"lattn",3))
+    end
+
+
+    a=eye(d);b=eye(d);
+    z=poly(0,"z");la=list();lb=list();
+    no=p-n-1;cv=cov;
+
+    if -no >= floor(l/d) then
+        error(msprintf(gettext("%s: Wrong values for input arguments #%d and #%d.\n"),"lattn",1, 2))
+    end
+
+    if no<0,
+        for j=1:-no,cv=[cov(j*d+1:(j+1)*d,:)';cv];end;
+        p=p-no;
+    end
+
+    if (p + 2 + (n-1)) * d >  size(cv,"r") then
+        error(msprintf(gettext("%s: Wrong values for input arguments #%d and #%d.\n"),"lattn",1, 2))
+    end
+
+    for j=0:n-1,
+        jd=jmat(j+1,d);
+        r1=jd*cv((p+1)*d+1:(p+2+j)*d,:);
+        r2=jd*cv(p*d+1:(p+1+j)*d,:);
+        r3=jd*cv((p-1-j)*d+1:p*d,:);
+        r4=jd*cv((p-j)*d+1:(p+1)*d,:);
+        c1=coeff(a);c2=coeff(b);
+        k1=(c1*r1)*inv(c2*r2);
+        k2=(c2*r3)*inv(c1*r4);
+        a1=a-k1*z*b;
+        b=-k2*a+z*b;a=a1;
+        la(j+1)=a;lb(j+1)=b;
+    end;
+endfunction
diff --git a/scilab/modules/signal_processing/macros/lattp.sci b/scilab/modules/signal_processing/macros/lattp.sci
new file mode 100644 (file)
index 0000000..117ad6b
--- /dev/null
@@ -0,0 +1,29 @@
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) INRIA - 1989 - G. Le Vey
+// Copyright (C) 2012 - Scilab Enterprises - Adeline CARNIS
+//
+// This file must be used under the terms of the CeCILL.
+// This source file is licensed as described in the file COPYING, which
+// you should have received as part of this distribution.  The terms
+// are also available at
+// http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
+
+function [la,lb]=lattp(n,p,cov)
+    // See lattn for more information
+    [l,d]=size(cov);
+    id=eye(d);
+    [a,b]=lattn(n,0,cov);a=a(n);b=b(n);
+    z=poly(0,"z");la=list();lb=list();
+    jd=jmat(n+1,d);
+    for j=0:p-1,
+        r1=jd*cov((j+1)*d+1:(j+n+2)*d,:);
+        r2=jd*cov(j*d+1:(j+n+1)*d,:);
+        c1=coeff(a);c2=coeff(b);
+        k=(c1*r1)*inv(c2*r2);
+        hst=-inv(c1(:,n*d+1:(n+1)*d));
+        r=k*hst;
+        a=(id-r*z)*a-k*z*b;
+        b=-hst*a;
+        la(j+1)=a;
+    end;
+endfunction