CMSCOPE changed

147 files changed, 5377 insertions, 0 deletions

diff --git a/modules/signal_processing/macros/%k.bin b/modules/signal_processing/macros/%k.bin
new file mode 100755
index 000000000..8b3f2cc1a
--- /dev/null
+++ b/modules/signal_processing/macros/%k.bin
diff --git a/modules/signal_processing/macros/%k.sci b/modules/signal_processing/macros/%k.sci
new file mode 100755
index 000000000..6cc05a101
--- /dev/null
+++ b/modules/signal_processing/macros/%k.sci
@@ -0,0 +1,39 @@
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) INRIA - F.D
+//
+// 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.1-en.txt
+
+function [K]=%k(m)
+    //K=%k(m)
+    //Calculates Jacobi's complete elliptic integral
+    //of the first kind:
+    //  K = integral from 0 to 1 of
+    //      [(1-t**2)(1-m*t**2)]**(-1/2)
+    //m is allowed to be a vector
+    //Ref :Abramowitz and Stegun page 598
+    //  m :Parameter used in calculating the elliptic
+    //    :integral where 0<m<1.
+    //  K :Value of the elliptic integral from 0 to 1
+    //    :on the real axis.
+    //
+    //!
+
+    [n1,n2]=size(m);
+    un=ones(n1,n2);
+    a=un;
+    b=sqrt(un-m);
+    c=sqrt(m);
+    while max(abs(c)) > %eps,
+        an=0.5*(a+b);
+        bn=sqrt(a.*b);
+        cn=0.5*(a-b);
+        a=an;
+        b=bn;
+        c=cn;
+    end,
+    K=%pi*un./(2*a);
+endfunction
diff --git a/modules/signal_processing/macros/%sn.bin b/modules/signal_processing/macros/%sn.bin
new file mode 100755
index 000000000..66874e852
--- /dev/null
+++ b/modules/signal_processing/macros/%sn.bin
diff --git a/modules/signal_processing/macros/%sn.sci b/modules/signal_processing/macros/%sn.sci
new file mode 100755
index 000000000..94139476d
--- /dev/null
+++ b/modules/signal_processing/macros/%sn.sci
@@ -0,0 +1,36 @@
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) INRIA - F.D
+//
+// 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.1-en.txt
+
+function [y]=%sn(x,m)
+    //Jacobi 's elliptic function with parameter m
+    //which computes the inverse of the elliptic
+    //integral for the parameter m.
+    //x may be a vector.
+    //The amplitude is computed in fortran and apply
+    //the addition formulas for elliptic functions
+    //  x :A point inside the fundamental rectangle
+    //    :defined by the elliptic integral
+    //  m :Parameter of the elliptic integral (0<m<1)
+    //  y :Result
+    //
+    //!
+
+    [n1,n2]=size(x);
+    n=n1*n2;
+    a=amell(real(x),sqrt(m));
+    s=sin(a);
+    c=cos(a);
+    d=sqrt(ones(n1,n2)-m*s.*s);
+    m1=1-m;
+    a1=amell(imag(x),sqrt(m1));
+    s1=sin(a1);
+    c1=cos(a1);
+    d1=sqrt(ones(n1,n2)-m1*s1.*s1);
+    y=(s.*d1+%i*c.*d.*s1.*c1)./(c1.*c1+m*s.*s.*s1.*s1);
+endfunction
diff --git a/modules/signal_processing/macros/analpf.bin b/modules/signal_processing/macros/analpf.bin
new file mode 100755
index 000000000..a3f641f40
--- /dev/null
+++ b/modules/signal_processing/macros/analpf.bin
diff --git a/modules/signal_processing/macros/analpf.sci b/modules/signal_processing/macros/analpf.sci
new file mode 100755
index 000000000..6418b850e
--- /dev/null
+++ b/modules/signal_processing/macros/analpf.sci
@@ -0,0 +1,62 @@
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) INRIA - 1988 - C. Bunks
+//
+// 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.1-en.txt
+
+function [hs,pols,zers,gain]=analpf(n,fdesign,rp,omega)
+    //[hs,pols,zers,gain]=analpf(n,fdesign,rp,omega)
+    //Creates analog low-pass filter with cut-off frequency at omega
+    //  n       :filter order (pos. integer)
+    //  fdesign :filter design method
+    //          :    fdesign=('butt','cheb1','cheb2','ellip')
+    //  rp      :2-vector of error values for cheb1, cheb2, and
+    //          :ellip filters where only rp(1) is used for
+    //          :cheb1 case, only rp(2) is used for cheb2 case, and
+    //          :rp(1) and rp(2) are both used for ellip case.
+    //          :          0<rp(1),rp(2)<1
+    //          :for cheb1 filters:  1-rp(1)<ripple<1 in passband
+    //          :for cheb2 filters:  0<ripple<rp(2)   in stopband
+    //          :for ellip filters:  1-rp(1)<ripple<1 in passband
+    //          :                    0<ripple<rp(2)   in stopband
+    //  omega   :cut-off frequency of low-pass filter in rd/s
+    //  hs      :rational polynomial transfer function
+    //  pols    :poles of transfer function
+    //  zers    :zeros of transfer function
+    //  gain    :gain of transfer function
+    //
+    //    hs=gain*poly(zers,'s')/poly(pols,'s')
+    //
+    //!
+
+
+    select fdesign
+    case "butt" then
+        [pols,gain]=zpbutt(n,omega);
+        zers=[];
+        hs=gain/real(poly(pols,"s"));
+    case "cheb1" then
+        epsilon=sqrt(-1+1/(1-rp(1))**2);
+        [pols,gain]=zpch1(n,epsilon,omega);
+        zers=[];
+        hs=gain/real(poly(pols,"s"));
+    case "cheb2" then
+        att=1/rp(2);
+        [zers,pols,gain]=zpch2(n,att,omega);
+        hs=gain*real(poly(zers,"s"))./real(poly(pols,"s"));
+    case "ellip" then
+        epsilon=sqrt(-1+1/(1-rp(1))**2);
+        att=1/rp(2);
+        m=find_freq(epsilon,att,n);
+        omegar=omega/sqrt(m);
+        [zers,pols,gain]=zpell(epsilon,att,omega,omegar);
+        hs=gain*real(poly(zers,"s"))./real(poly(pols,"s"));
+    else
+        error(msprintf(gettext("%s: Wrong value for input argument #%d: Must be in the set {%s}.\n"),..
+        "analpf",2,"''butt'',''cheb1'',''cheb2'',''ellip''"));
+    end
+    hs.dt="c";
+endfunction
diff --git a/modules/signal_processing/macros/bilt.bin b/modules/signal_processing/macros/bilt.bin
new file mode 100755
index 000000000..26030534a
--- /dev/null
+++ b/modules/signal_processing/macros/bilt.bin
diff --git a/modules/signal_processing/macros/bilt.sci b/modules/signal_processing/macros/bilt.sci
new file mode 100755
index 000000000..dbf7dd0a8
--- /dev/null
+++ b/modules/signal_processing/macros/bilt.sci
@@ -0,0 +1,100 @@
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) INRIA - 1989 - C. Bunks
+// Copyright (C) INRIA - 1997 - C. Bunks
+//
+// 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.1-en.txt
+
+function [npl,nzr,ngn]=bilt(pl,zr,gn,num,den)
+    //[npl,nzr,ngn]=bilt(pl,zr,gn,num,den)
+    //macro for calculating the gain poles and zeros
+    //which result from a bilinear transform or from
+    //a biquadratic transform.  Used by the macros iir
+    //and trans
+    //Note: ***This macro is not intended for general use***
+    //  pl  :input poles
+    //  zr  :input zeros
+    //  gn  :input gain
+    //  num :numerator of transform
+    //  den :denominator of transform
+    //!
+    if type(pl)<>1 then
+        error(msprintf(gettext("%s: Wrong type for input argument #%d: Array of floating point numbers expected.\n") ,"bilt",1))
+    end
+    if type(zr)<>1 then
+        error(msprintf(gettext("%s: Wrong type for input argument #%d: Array of floating point numbers expected.\n") ,"bilt",2))
+    end
+    if type(gn)<>1 then
+        error(msprintf(gettext("%s: Wrong type for input argument #%d: Array of floating point numbers expected.\n") ,"bilt",3))
+    end
+    if size(gn,"*")<>1 then
+        error(msprintf(gettext("%s: Wrong size for input argument #%d: A scalar expected.\n") ,"bilt",3))
+    end
+    if type(num)<>2|size(num,"*")<>1 then
+        error(msprintf(gettext("%s: Wrong size for input argument #%d: A polynomial expected.\n") ,"bilt",4))
+    end
+    if type(den)<>2|size(den,"*")<>1 then
+        error(msprintf(gettext("%s: Wrong size for input argument #%d: A polynomial expected.\n") ,"bilt",5))
+    end
+
+    order=degree(den)
+
+    if and(order<>[1 2]) then
+        error(msprintf(gettext("%s: Wrong values for input argument #%d: degree must be in the set {%s}.\n"),"bilt",5,"1,2"))
+    end
+    if degree(num)<>order then
+        error(msprintf(gettext("%s: Incompatible input arguments #%d and #%d: Same degree expected.\n"),"bilt",4,5))
+    end
+    n=coeff(num);
+    d=coeff(den);
+    ms=max(size(pl));ns=max(size(zr));
+
+    select order
+    case 1  then
+        n0=n(1);n1=n(2);
+        d0=d(1);d1=d(2);
+        if zr == [] then
+            ngn=1/prod(n1*ones(pl)-d1*pl);
+        else
+            ngn=prod(n1*ones(zr)-d1*zr)/prod(n1*ones(pl)-d1*pl);
+        end
+        if ms<>ns then ngn=real(gn*d1**(ms-ns)*ngn);else ngn=real(gn*ngn);end
+        nzr=-(n0*ones(zr)-d0*zr)./(n1*ones(zr)-d1*zr);
+        npl=-(n0*ones(pl)-d0*pl)./(n1*ones(pl)-d1*pl);
+        if ms>ns then
+            nzr=[nzr';-(d0/d1)*ones(ms-ns,1)];
+        elseif ms<ns then
+            npl=[npl';-(d0/d1)*ones(ms-ns,1)];
+        end
+    case 2 then
+        n0=n(1);n1=n(2);n2=n(3);
+        d0=d(1);d1=d(2);d2=d(3);
+        a=n2*ones(zr)-d2*zr;
+        b=n1*ones(zr)-d1*zr;
+        c=n0*ones(zr)-d0*zr;
+        gz=a;
+        z1=-b./(2*a)+sqrt((b./(2*a)).^2-c./a);
+        z2=-b./(2*a)-sqrt((b./(2*a)).^2-c./a);
+        a=n2*ones(pl)-d2*pl;
+        b=n1*ones(pl)-d1*pl;
+        c=n0*ones(pl)-d0*pl;
+        gp=a;
+        p1=-b./(2*a)+sqrt((b./(2*a)).^2-c./a);
+        p2=-b./(2*a)-sqrt((b./(2*a)).^2-c./a);
+        gw=d2;
+        w1=-d1./(2*d2)+sqrt((d1./(2*d2))**2-d0./d2);
+        w2=-d1./(2*d2)-sqrt((d1./(2*d2))**2-d0./d2);
+        ngn=gn*prod(gz)/prod(gp);
+        nzr=[z1,z2];
+        npl=[p1,p2];
+        if ms>ns then
+            nzr=[nzr';-(d0/d1)*ones(ms-ns,1)];
+        elseif ms<ns then
+            npl=[npl';-(d0/d1)*ones(ms-ns,1)];
+        end
+        ngn=real(ngn*(gw**(ms-ns)));
+    end
+endfunction
diff --git a/modules/signal_processing/macros/buildmacros.bat b/modules/signal_processing/macros/buildmacros.bat
new file mode 100755
index 000000000..c4e35ec40
--- /dev/null
+++ b/modules/signal_processing/macros/buildmacros.bat
@@ -0,0 +1 @@
+@..\..\..\bin\scilex -nwni -ns -e exec('buildmacros.sce');quit;
\ No newline at end of file
diff --git a/modules/signal_processing/macros/buildmacros.sce b/modules/signal_processing/macros/buildmacros.sce
new file mode 100755
index 000000000..8a24332a6
--- /dev/null
+++ b/modules/signal_processing/macros/buildmacros.sce
@@ -0,0 +1,9 @@
+//------------------------------------
+// Allan CORNET INRIA 2005
+//------------------------------------
+if (isdef("genlib") == %f) then
+    exec(SCI+"/modules/functions/scripts/buildmacros/loadgenlib.sce");
+end
+//------------------------------------
+genlib("signal_processinglib","SCI/modules/signal_processing/macros",%f,%t);
+//------------------------------------
diff --git a/modules/signal_processing/macros/buttmag.bin b/modules/signal_processing/macros/buttmag.bin
new file mode 100755
index 000000000..ff7c1b3dd
--- /dev/null
+++ b/modules/signal_processing/macros/buttmag.bin
diff --git a/modules/signal_processing/macros/buttmag.sci b/modules/signal_processing/macros/buttmag.sci
new file mode 100755
index 000000000..0d83af242
--- /dev/null
+++ b/modules/signal_processing/macros/buttmag.sci
@@ -0,0 +1,31 @@
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) INRIA - F.D
+//
+// 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.1-en.txt
+
+function [h]=buttmag(order,omegac,sample_vector)
+    //<h>=buttmag(order,omegac,sample_vector)
+    //Squared Magnitude response of a Butterworth filter
+    //omegac = cutoff frequency ; sample_vector = sample of frequencies
+    //  order  :Filter order
+    //  omegac :Cut-off frequency in Hertz
+    //  sample_vector :Vector of frequency where buttmag is evaluated
+    //  h      :Butterworth filter values at sample points
+    //
+    //!
+
+    // For ascendant compatibility before bug 4618 fix
+    // http://bugzilla.scilab.org/show_bug.cgi?id=4618
+    // In case the users calls buttmag with named arguments
+    if exists("sample","local")==1 then
+        sample_vector = sample;
+    end
+
+    [n1,n2]=size(sample_vector);
+    un=ones(n1,n2);
+    h=un./(un+(sample_vector/omegac).^(2*order));
+endfunction
diff --git a/modules/signal_processing/macros/casc.bin b/modules/signal_processing/macros/casc.bin
new file mode 100755
index 000000000..3e67a0e5a
--- /dev/null
+++ b/modules/signal_processing/macros/casc.bin
diff --git a/modules/signal_processing/macros/casc.sci b/modules/signal_processing/macros/casc.sci
new file mode 100755
index 000000000..b038e3199
--- /dev/null
+++ b/modules/signal_processing/macros/casc.sci
@@ -0,0 +1,44 @@
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) INRIA - 1988 - F.D
+//
+// 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.1-en.txt
+
+function cels=casc(x,z)
+    //cels=casc(x,z)
+    //Creates cascade realization of filter
+    //from a matrix of coefficients
+    //  x    :(4xN)-Matrix where each column is a cascade
+    //       :element, the first two column entries being
+    //       :the numerator coefficients and the second two
+    //       :column entries being the denominator coefficients
+    //  z    :Character string representing the cascade variable
+    //  cels :Resulting cascade representation
+    //
+    //EXAMPLE:
+    //  x=[ 1.     2.     3.  ;
+    //      4.     5.     6.  ;
+    //      7.     8.     9.  ;
+    //      10.    11.    12. ]
+    //
+    //  cels=casc(x,'z')
+    //  cels      =
+    //
+    //  !             2               2               2  !
+    //  !   1 + 4z + z      2 + 5z + z      3 + 6z + z   !
+    //  !  ------------    ------------    ------------  !
+    //  !              2               2               2 !
+    //  !   7 + 10z + z     8 + 11z + z     9 + 12z + z  !
+    //!
+
+    cels=[];
+    for col=x,
+        nf=[col(1:2);1];
+        nd=[col(3:4);1];
+        cels=[cels,syslin([],poly(nf,"z","c"),poly(nd,"z","c"))];
+    end,
+
+endfunction
diff --git a/modules/signal_processing/macros/cepstrum.bin b/modules/signal_processing/macros/cepstrum.bin
new file mode 100755
index 000000000..4c0511d34
--- /dev/null
+++ b/modules/signal_processing/macros/cepstrum.bin
diff --git a/modules/signal_processing/macros/cepstrum.sci b/modules/signal_processing/macros/cepstrum.sci
new file mode 100755
index 000000000..d7cb3f662
--- /dev/null
+++ b/modules/signal_processing/macros/cepstrum.sci
@@ -0,0 +1,89 @@
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) INRIA
+//
+// 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.1-en.txt
+
+function fresp = cepstrum(w,mag)
+    // Uses the complex-cepstrum (Oppenheim & Schafer, Digital
+    //Signal Processing, p. 501) to generate, at the frequencies
+    //w, a complex frequency response fresp whose magnitude is
+    //equal to magnitude data mag and whose phase corresponds
+    //to a stable, minimum phase transfer function.
+
+    if ~isreal(w) then
+        error(msprintf(gettext("%s: Input argument #%d must be real.\n"),"cepstrum",1));
+    end
+    if ~isreal(mag) then
+        error(msprintf(gettext("%s: Input argument #%d must be real.\n"),"cepstrum",2));
+    end
+
+    [w,ind]=gsort(-w);w=-w;mag=mag(ind);
+    dnum=length(w);
+
+    pw=sqrt(w(dnum)*w(1));
+    if w(1)<=0 then
+        error(msprintf(gettext("%s: Wrong values for input argument #%d: Elements must be positive.\n"),"cepstrum",1));
+    end
+    if pw<1.d-8 then
+        error(msprintf(gettext("%s: Wrong values for input argument #%d: Distance between Max and Min elements must be greater than %g.\n"),..
+        "cepstrum",1,1d-8));
+    end
+    if size(w,"*")<>size(mag,"*") then
+        error(msprintf(gettext("%s: Incompatible input arguments #%d and #%d: Same numbers of elements expected.\n"),"cepstrum",1,2));
+    end
+
+    tcomes=(w.*w)/pw/pw;
+
+    dw=acos((1-tcomes)./(1+tcomes));
+    lowf=dw(1);
+    highf=dw(length(dw));
+    if lowf <= 0 | highf >=%pi
+        error(msprintf(gettext("%s: Wrong values for input argument #%d: Elements must be positive.\n"),"cepstrum",1));
+    end
+
+    nn=ceil((log(2*%pi/min(lowf,%pi-highf)))/log(2));
+    npts=2*(2^nn);hnpts =2^nn;
+    if npts<4096 then npts=4096;hnpts=2048;end
+
+    lmagin =(1/log(10))*log(mag);lindw =(%pi/hnpts)*(0:hnpts-1);
+    lmag=zeros(1,hnpts);topval=length(dw);
+    p=find(lindw<dw(1));lmag(p)=lmagin(1)*ones(1,length(p));
+    p=find(dw(topval)<=lindw);
+    lmag(p)=lmagin(topval)*ones(1,length(p));
+    for i=2:topval
+        p=find(lindw>=dw(i-1) & lindw<dw(i));
+        wrat=lindw(p) - dw(i-1)*ones(1,length(p));
+        wrat=(1/(dw(i)-dw(i-1)))*wrat;
+        lmag(p)=(ones(1,length(p))-wrat)*lmagin(i-1)+wrat*lmagin(i);
+    end
+    linmag=exp(log(10)*lmag);
+
+    dome=[lindw,(2*%pi)-fliplr(lindw)];
+    mag=[linmag,fliplr(linmag)];
+
+    ymag=log(mag.*mag);ycc=fft(ymag,1);nptso2=npts/2;xcc=ycc(1:nptso2);
+    xcc(1)=xcc(1)/2;xhat=exp(fft(xcc,-1));domeg=dome(1:2:nptso2-1);
+    xhat = xhat(1:nptso2/2);nptsslo=length(dw);nptsfst = length(domeg);
+
+    if domeg(1)<=dw(1) & domeg(nptsfst)>=dw(nptsslo)
+        fresp=zeros(1,nptsslo);
+        for i=1:nptsslo
+            p=min(find(domeg>=dw(i)));
+            wrat=(dw(i)-domeg(p-1))/(domeg(p)-domeg(p-1));
+            fresp(i)=wrat*xhat(p) + (1-wrat)*xhat(p-1);
+        end
+    else
+        error(msprintf(gettext("%s: Wrong values for input argument #%d: Not sampled high enough.\n"),"cepstrum",1));
+    end
+    fresp=fresp(:);
+
+endfunction
+function m=fliplr(m)
+    //utility fct.
+    [p,q]=size(m);
+    m=m(:,q:-1:1);
+endfunction
diff --git a/modules/signal_processing/macros/cheb1mag.bin b/modules/signal_processing/macros/cheb1mag.bin
new file mode 100755
index 000000000..4f1812d6f
--- /dev/null
+++ b/modules/signal_processing/macros/cheb1mag.bin
diff --git a/modules/signal_processing/macros/cheb1mag.sci b/modules/signal_processing/macros/cheb1mag.sci
new file mode 100755
index 000000000..ea2b312de
--- /dev/null
+++ b/modules/signal_processing/macros/cheb1mag.sci
@@ -0,0 +1,30 @@
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) INRIA - F.D
+//
+// 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.1-en.txt
+function [h2]=cheb1mag(n,omegac,epsilon,sampleFreq)
+    //<h2>=cheb1mag(n,omegac,epsilon,sample)
+    //Square magnitude response of a type 1 Chebyshev filter
+    //omegac=passband edge
+    //epsilon such that 1/(1+epsilon**2)=passband ripple
+    //sample vector of frequencies where the square magnitude
+    //is desired.
+    //  n       :Filter order
+    //  omegac  :Cut-off frequency
+    //  epsilon :Ripple in pass band
+    //  sample  :Vector of frequency where cheb1mag is evaluated
+    //  h2      :Chebyshev I filter values at sample points
+    //
+    //!
+
+
+    [n1,n2]=size(sampleFreq);
+    un=ones(n1,n2);
+    Tn=chepol(n,"x");  //n-th Chebyshev polynomial
+    fr=freq(Tn,1,sampleFreq/omegac);   //fr=Tn(sample/omegac)
+    h2=un./(un+epsilon*epsilon*fr.*fr)   //magnitude
+endfunction
diff --git a/modules/signal_processing/macros/cheb2mag.bin b/modules/signal_processing/macros/cheb2mag.bin
new file mode 100755
index 000000000..d3831bf99
--- /dev/null
+++ b/modules/signal_processing/macros/cheb2mag.bin
diff --git a/modules/signal_processing/macros/cheb2mag.sci b/modules/signal_processing/macros/cheb2mag.sci
new file mode 100755
index 000000000..89c3a8003
--- /dev/null
+++ b/modules/signal_processing/macros/cheb2mag.sci
@@ -0,0 +1,30 @@
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) INRIA - F.D
+//
+// 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.1-en.txt
+function [h2]=cheb2mag(n,omegar,A,samplefreq)
+    //<h2>=cheb2mag(n,omegar,A,samplefreq)
+    //Square magnitude response of a type 1 Chebyshev filter
+    //omegar = stopband edge
+    //samplefreq = vector of frequencies where the square magnitude
+    //h2 is desired.
+    //  n       :Filter order
+    //  omegar  :Cut-off frequency
+    //  A       :Attenuation in stop band
+    //  samplefreq  :Vector of frequency where cheb2mag is evaluated
+    //  h2      :Chebyshev II filter values at sample points
+    //
+    //!
+
+
+    [n1,n2]=size(samplefreq);
+    un=ones(n1,n2);
+    Tn=chepol(n,"x");             //n-th Chebyshev polynomial
+    frd=freq(Tn,1,omegar*un./samplefreq);   //frd=Tn(omegar/samplefreq)
+    h2=un./(un+(A*A-1)*un./real(frd.*frd))
+endfunction
+
diff --git a/modules/signal_processing/macros/cleanmacros.bat b/modules/signal_processing/macros/cleanmacros.bat
new file mode 100755
index 000000000..5079dfd71
--- /dev/null
+++ b/modules/signal_processing/macros/cleanmacros.bat
@@ -0,0 +1,3 @@
+@del *.bin 2>NUL
+@del lib 2>NUL
+@del names 2>NUL
\ No newline at end of file
diff --git a/modules/signal_processing/macros/conv.bin b/modules/signal_processing/macros/conv.bin
new file mode 100755
index 000000000..f61e13e83
--- /dev/null
+++ b/modules/signal_processing/macros/conv.bin
diff --git a/modules/signal_processing/macros/conv.sci b/modules/signal_processing/macros/conv.sci
new file mode 100755
index 000000000..50de64887
--- /dev/null
+++ b/modules/signal_processing/macros/conv.sci
@@ -0,0 +1,34 @@
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2012 - INRIA - Serge STEER
+//
+// 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.1-en.txt
+function h=conv(u,v,Shape)
+    if and(size(u)>1) then
+        error( msprintf(_("%s: Wrong size for argument #%d: Vector expected.\n"),"conv",1))
+    end
+    if and(size(v)>1) then
+        error( msprintf(_("%s: Wrong size for argument #%d: Vector expected.\n"),"conv",2))
+    end
+    if argn(2)==2 then
+        Shape="full",
+    elseif and(Shape<>["full","same","valid"]) then
+        error(msprintf(_("%s: Wrong value for input argument #%d: ""%s"" or ""%s"" expected.\n"),"conv",3, """full"", ""same""","""valid"""));
+    end
+
+    h=conv2(u(:),v(:),Shape);
+    //set result orientation
+    if Shape=="full" then
+        if size(u,"*")>size(v,"*") then
+            if size(u,1)==1 then h=matrix(h,1,-1);end
+        else
+            if size(v,1)==1 then h=matrix(h,1,-1);end
+        end
+    else
+        if size(u,1)==1 then h=matrix(h,1,-1);end
+    end
+endfunction
+
diff --git a/modules/signal_processing/macros/convol.bin b/modules/signal_processing/macros/convol.bin
new file mode 100755
index 000000000..1d55e4c83
--- /dev/null
+++ b/modules/signal_processing/macros/convol.bin
diff --git a/modules/signal_processing/macros/convol.sci b/modules/signal_processing/macros/convol.sci
new file mode 100755
index 000000000..d5a1c9945
--- /dev/null
+++ b/modules/signal_processing/macros/convol.sci
@@ -0,0 +1,57 @@
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) INRIA
+//
+// 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.1-en.txt
+function [y,e1]=convol(h,x,e0)
+    //  convol - convolution
+    //%CALLING SEQUENCE
+    //  [y]=convol(h,x)
+    //  [y,e1]=convol(h,x,e0)     (for use with overlap add method)
+    //%PARAMETERS
+    //  x,h       :input sequences (h is a "short" sequence, x a "long" one)
+    //  e0        : old tail to overlap add (not used in first call)
+    //  y         : output of convolution
+    //  e1        : new tail to overlap add (not used in last call)
+    //%DESCRIPTION
+    //  calculates the convolution y= h*x of two discrete sequences by
+    //  using the fft.  overlap add method can be used.
+    //%USE OF OVERLAP ADD METHOD
+    //  For x=[x1,x2,...,xNm1,xN]
+    //  First call : [y1,e1]=convol(h,x1)
+    //  Subsequent calls : [yk,ek]=convol(h,xk,ekm1)
+    //  Final call : [yN]=convol(h,xN,eNm1)
+    //  Finally y=[y1,y2,...,yNm1,yN]
+    //!
+
+    [lhs,rhs]=argn(0)
+    n=prod(size(x))
+    m=prod(size(h))
+    m1=2^(int(log(n+m-1)/log(2))+1)
+    x(m1)=0;h(m1)=0
+    if norm(imag(x))==0&norm(imag(h))==0 then
+        y=real(fft(fft(matrix(x,1,m1),-1).*fft(matrix(h,1,m1),-1),1))
+    else
+        y=fft(fft(matrix(x,1,m1),-1).*fft(matrix(h,1,m1),-1),1)
+    end
+    if lhs+rhs==5 then,
+        e0(n)=0;//update carried from left to right
+        e1=y(n+1:n+m-1)
+        y=y(1:n)+e0
+
+    elseif lhs+rhs==4 then
+        if rhs==2 then
+            e1=y(n+1:n+m-1)
+            y=y(1:n) //initial update
+        else
+            e0(n+m-1)=0 //final update
+            y=y(1:n+m-1)+e0
+        end,
+
+    else
+        y=y(1:n+m-1) //no update
+    end
+endfunction
diff --git a/modules/signal_processing/macros/convol2d.bin b/modules/signal_processing/macros/convol2d.bin
new file mode 100755
index 000000000..5b9d19bad
--- /dev/null
+++ b/modules/signal_processing/macros/convol2d.bin
diff --git a/modules/signal_processing/macros/convol2d.sci b/modules/signal_processing/macros/convol2d.sci
new file mode 100755
index 000000000..7f322ac3f
--- /dev/null
+++ b/modules/signal_processing/macros/convol2d.sci
@@ -0,0 +1,56 @@
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2012 - INRIA - Serge STEER
+//
+// 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.1-en.txt
+function y=convol2d(h,x)
+    //  convol2d - 2-D convolution
+    //%CALLING SEQUENCE
+    //  y=convol2d(h,x)
+    //%PARAMETERS
+    //  x,h       :input matrices
+    //  y         : output of convolution
+    //%DESCRIPTION
+    //  calculates the 2-D convolution y= h*x of two discrete sequences by
+    //  using the fft.
+    if argn(2)<2 then
+        error(msprintf(_("%s: Wrong number of input arguments: %d expected.\n"),"convol2d",2))
+    end
+
+    if type(h)<>1 then
+        error(msprintf(_("%s: Wrong type for argument #%d: Real or complex matrix expected.\n"),"convol2d",1))
+    end
+    if type(x)<>1 then
+        error(msprintf(_("%s: Wrong type for argument #%d: Real or complex matrix expected.\n"),"convol2d",2))
+    end
+    if isempty(h) then
+        y=zeros(x);
+        return
+    end
+    if isempty(x) then
+        y=zeros(h);
+        return
+    end
+    //inline fft2d function definition (fft2 does not manage inverse fft)
+    function y=fft2d(x,d)
+        [mx,nx]=size(x)
+        y=fft(fft(x,d,mx,1),d,nx,mx)
+    endfunction
+    [mx,nx]=size(x);
+    [mh,nh]=size(h);
+    //use  power of 2 dimensions for efficiency
+    m1=2^(int(log(mx+mh-1)/log(2))+1);
+    n1=2^(int(log(nx+nh-1)/log(2))+1);
+    //m1=mx+mh-1;
+    //n1=nx+nh-1;
+    x(m1,n1)=0;
+    h(m1,n1)=0;
+    y=fft2d(fft2d(x,-1).*fft2d(h,-1),1);
+    if isreal(h,0)&isreal(x,0) then
+        y=real(y);
+    end
+    y=y(1:(mx+mh-1),1:(nx+nh-1));
+endfunction
diff --git a/modules/signal_processing/macros/cspect.bin b/modules/signal_processing/macros/cspect.bin
new file mode 100755
index 000000000..27e60c0ab
--- /dev/null
+++ b/modules/signal_processing/macros/cspect.bin
diff --git a/modules/signal_processing/macros/cspect.sci b/modules/signal_processing/macros/cspect.sci
new file mode 100755
index 000000000..17a910dd7
--- /dev/null
+++ b/modules/signal_processing/macros/cspect.sci
@@ -0,0 +1,110 @@
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) INRIA - 1988 - C. Bunks
+//
+// 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.1-en.txt
+
+function [sm,cwp]=cspect(nlags,ntp,wtype,x,y,wpar)
+    //[sm,cwp]=cspect(nlags,ntp,wtype,x,y,wpar)
+    //Spectral estimation using the correlation method.
+    //Cross-spectral estimate of x and y is calculated when both
+    //x and y are given.  Auto-spectral estimate of x is calculated
+    //if y is not given.
+    //  x        :data if vector, amount of input data if scalar
+    //  y        :data if vector, amount of input data if scalar
+    //  nlags    :number of correlation lags (pos. integer)
+    //  ntp      :number of transform points (pos. integer)
+    //  wtype    :window type ('re','tr','hm','hn','kr','ch')
+    //  wpar     :optional window parameters
+    //           :         for wtype='kr', wpar>0
+    //           :         for wtype='ch', 0<wpar(1)<.5, wpar(2)>0
+    //  sm       :power spectral estimate in the interval [0,1]
+    //  cwp      :calculated value of unspecified Chebyshev
+    //           :window parameter
+    //
+    //!
+
+    [lhs,rhs]=argn(0);
+
+    if and(wtype<>["re","tr","hm","hn","kr","ch"]) then
+        error(msprintf(gettext("%s: Wrong value for input argument #%d: Must be in the set {%s}.\n"),..
+        "cspect",3,"''re'',''tr'',''hm'',''hn'',''kr'',''ch''"));
+    end
+
+    //Analyze calling sequence and construct window
+    if rhs==4 then,
+        //cspect(nlags,ntp,wtype,x)
+        [w,cwp]=window(wtype,2*nlags-1);
+        crossFlag=%f; //autocorrelation
+    elseif rhs==5 then
+        //cspect(nlags,ntp,wtype,x,wpar) or cspect(nlags,ntp,wtype,x,y)
+        if wtype=="kr" then //cspect(nlags,ntp,'kr',x,wpar)
+            wpar=y;
+            [w,cwp]=window(wtype,2*nlags-1,wpar);
+            crossFlag=%f; //autocorrelation
+        elseif wtype=="ch" then  //cspect(nlags,ntp,'ch',x,wpar)
+            wpar=y;
+            [w,cwp]=window(wtype,2*nlags-1,wpar);
+            crossFlag=%f; //autocorrelation
+        else,//cspect(nlags,ntp,wtype,x,y)
+            crossFlag=%t;//cross correlation
+            [w,cwp]=window(wtype,2*nlags-1);
+        end,
+    else,//cspect(nlags,ntp,wtype,x,y,wpar)
+        [w,cwp]=window(wtype,2*nlags-1,wpar);
+        crossFlag=%t;//cross correlation
+    end,
+
+
+    //Make x and y row vectors
+    x=matrix(x,1,-1);
+    if crossFlag then
+        y=matrix(y,1,-1);
+        if size(x,"*")<>size(y,"*") then
+            error(msprintf(gettext("%s: Arguments #%d and #%d must have the same sizes.\n"),"pspect",4,5));
+        end
+    end
+
+    //Estimate correlations
+    if size(x,"*")==1 then  //Batch processing of x and y data
+        nsects=int(x/(3*nlags));
+        xlen=int(x/nsects);
+        ss=zeros(1,2*nlags);
+        if crossFlag then,
+            for k=1:nsects
+                xk=getx(xlen,1+(k-1)*xlen);
+                yk=gety(xlen,1+(k-1)*xlen);
+                ss=corr("update",xk,yk,ss);
+            end
+            re=fft(ss,1)/x;
+            re=[re(nlags:-1:1) re(2:nlags)];
+        else
+            for k=1:nsects
+                xk=getx(xlen,1+(k-1)*xlen);
+                ss=corr("update",xk,ss);
+            end
+            re=fft(ss,1)/x;
+            re=[re(nlags:-1:1) re(2:nlags)];
+        end
+    else // Signal data given in x and y if cross-correlation
+        if crossFlag then
+            [re1,me]=corr(x,y,nlags);
+            [re2,me]=corr(y,x,nlags);
+            re=[re1(nlags:-1:1) re2(2:nlags)];
+        else
+            [re,me]=corr(x,nlags);
+            re=[re(nlags:-1:1) re(2:nlags)];
+        end
+    end
+
+    //Window correlation estimate
+    wre=w.*re;
+
+    //Fourier transform to obtain spectral estimate
+    wree=[wre zeros(1,ntp-2*nlags+1)];
+    sm=fft(wree,-1)
+    //sm=abs(fft(wree,-1));
+endfunction
diff --git a/modules/signal_processing/macros/czt.bin b/modules/signal_processing/macros/czt.bin
new file mode 100755
index 000000000..de3d1fb9a
--- /dev/null
+++ b/modules/signal_processing/macros/czt.bin
diff --git a/modules/signal_processing/macros/czt.sci b/modules/signal_processing/macros/czt.sci
new file mode 100755
index 000000000..e61e602d3
--- /dev/null
+++ b/modules/signal_processing/macros/czt.sci
@@ -0,0 +1,57 @@
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) INRIA - 1988 - C. Bunks
+//
+// 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.1-en.txt
+
+function [czx]=czt(x,m,w,phi,a,theta)
+    //czx=czt(x,m,w,phi,a,theta)
+    //chirp z-transform algorithm which calcultes the
+    //z-transform on a spiral in the z-plane at the points
+    //[a*exp(j*theta)][w**kexp(j*k*phi)] for k=0,1,...,m-1.
+    //  x     :Input data sequence
+    //  m     :czt is evaluated at m points in z-plane
+    //  w     :Magnitude multiplier
+    //  phi   :Phase increment
+    //  a     :Initial magnitude
+    //  theta :Initial phase
+    //  czx   :Chirp z-transform output
+    //!
+
+    //get the size of x and find the maximum of (n,m)
+
+    n=max(size(x));
+    nm=max([n,m]);
+
+    //create sequence h(n)=[w*exp(-j*phi)]**(-n*n/2)
+
+    w=w*exp(-%i*phi);
+    idx=(-nm+1:0);
+    idx=-idx.*idx/2;
+    h=exp(idx*log(w));
+    h(nm+1:2*nm-1)=h(nm-1:-1:1);
+
+    //create g(n)
+
+    a=a*exp(%i*theta);
+    idx=(0:n-1);
+    g=exp(-idx*log(a))./h(nm:nm+n-1);
+    g=x.*g;
+
+    //convolve h(n) and g(n)
+
+    hc=h(nm:nm+m-1);
+    hc(m+1:m+n-1)=h(nm-n+1:nm-1);
+    gc=0*ones(hc);
+    gc(1:n)=g;
+    hf=fft(hc,-1);
+    gf=fft(gc,-1);
+    hcg=fft(hf.*gf,1);
+
+    //preserve m points and divide by h(n)
+
+    czx=hcg(1:m)./h(nm:nm+m-1);
+endfunction
diff --git a/modules/signal_processing/macros/detrend.bin b/modules/signal_processing/macros/detrend.bin
new file mode 100755
index 000000000..88fa52672
--- /dev/null
+++ b/modules/signal_processing/macros/detrend.bin
diff --git a/modules/signal_processing/macros/detrend.sci b/modules/signal_processing/macros/detrend.sci
new file mode 100755
index 000000000..4c299c899
--- /dev/null
+++ b/modules/signal_processing/macros/detrend.sci
@@ -0,0 +1,97 @@
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) Bruno Pincon
+//
+// 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.1-en.txt
+//
+function [y] = detrend(x, flag, bp)
+    //
+    // this function removes the constant or linear or
+    // piecewise linear trend from a vector x. If x is
+    // matrix this function removes the trend of each
+    // column of x.
+    //
+    // flag = "constant" or "c" to removes the constant trend
+    //        (simply the mean of the signal)
+    // flag = "linear" or "l" to remove the linear or piecewise
+    //        linear trend.
+    //
+    // To define the piecewise linear trend the function needs the
+    // breakpoints and these must be given as the third argument bp.
+    //
+    // The "instants" of the signal x are in fact from 0 to m-1
+    // (m = length(x) if x is a vector and m = size(x,1) in case
+    // x is a matrix). So bp must be reals in [0 m-1].
+    //
+
+    rhs = argn(2)
+    if rhs < 1 | rhs > 3 then
+        error(msprintf(gettext("%s: Wrong number of input arguments: %d to %d expected.\n"),"detrend",1,3));
+    elseif rhs == 1
+        flag = "linear"; bp = []
+    elseif rhs == 2
+        bp = []
+    end
+
+    if type(x)~=1 then
+        error(msprintf(gettext("%s: Wrong type for input argument #%d: Array of floating point numbers expected.\n"),"detrend",1));
+    end
+    if type(flag)~=10 then
+        error(msprintf(gettext("%s: Wrong type for input argument #%d: Array of floating point numbers expected.\n"),"detrend",2));
+    end
+    if ~(type(bp)==1 & isreal(bp)) then
+        error(msprintf(gettext("%s: Wrong type for input argument #%d: Array of floating point numbers expected.\n"),"detrend",3));
+    end
+
+    [mx,nx] = size(x)
+    if mx==1 | nx==1 then
+        x_is_vector = %t; x = x(:); m = mx*nx; n = 1
+    elseif mx*nx == 0 then
+        y = []; return
+    else
+        x_is_vector = %f; m = mx; n = nx
+    end
+
+
+    if flag == "constant" | flag == "c" then
+        y = x - ones(m,1)*mean(x,"r")
+    elseif flag == "linear" | flag == "l"
+        bp = unique([0 ; bp(:) ; m-1])
+        // delete all the breakpoints outside [0,m-1]
+        while bp(1) < 0, bp(1)=[], end
+        while bp($) > m-1, bp($)=[], end
+        // breakpoints are 0-based so add one to
+        // compare them with signal vector indices (1-based)
+        bp = bp + 1;
+        nbp = length(bp);
+        // build the least square matrix with hat functions
+        // (as a basis for continuous piecewise linear functions)
+        A = zeros(m, nbp)
+        k1 = 1
+        delta_bp = diff(bp)
+        for j = 2:nbp-1
+            k2 = ceil(bp(j))-1
+            ind = (k1:k2)'
+            A(ind,j-1) = (bp(j) - ind)/delta_bp(j-1)
+            A(ind,j) = (ind - bp(j-1))/delta_bp(j-1)
+            k1 = k2+1
+        end
+        ind = (k1:m)'
+        A(ind,nbp-1) = (m - ind)/delta_bp(nbp-1)
+        A(ind,nbp) = (ind - bp(nbp-1))/delta_bp(nbp-1)
+        // solve the least square pb and retrieve the fitted
+        // piecewise linear func off the signal
+        y = x - A*(A\x)
+    else
+        error(msprintf(gettext("%s: Wrong value for input argument #%d: Must be in the set {%s}.\n") ,..
+        "detrend",2,"''constant'',''c'',''linear'',''l''"));
+    end
+
+    if x_is_vector then
+        y = matrix(y,mx,nx)
+    end
+
+endfunction
diff --git a/modules/signal_processing/macros/ell1mag.bin b/modules/signal_processing/macros/ell1mag.bin
new file mode 100755
index 000000000..b74c849f9
--- /dev/null
+++ b/modules/signal_processing/macros/ell1mag.bin
diff --git a/modules/signal_processing/macros/ell1mag.sci b/modules/signal_processing/macros/ell1mag.sci
new file mode 100755
index 000000000..15cbfc7c7
--- /dev/null
+++ b/modules/signal_processing/macros/ell1mag.sci
@@ -0,0 +1,23 @@
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) INRIA
+//
+// 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.1-en.txt
+
+function [v]=ell1mag(eps,m1,z)
+    //Function used for squared magnitude of an elliptic filter
+    //Usually m1=eps*eps/(a*a-1);
+    //  eps     :Passband ripple=1/(1+eps**2)
+    //  m1      :Stopband ripple=1/(1+(eps**2)/m1)
+    //  z       :Sample vector of values in the complex plane
+    //  v       :Elliptic filter values at sample points
+    //
+    //!
+    s=%sn(z,m1);un=ones(z);
+    v=real(un./(un+eps*eps*s.*s))
+
+
+endfunction
diff --git a/modules/signal_processing/macros/eqfir.bin b/modules/signal_processing/macros/eqfir.bin
new file mode 100755
index 000000000..4ef4ee94e
--- /dev/null
+++ b/modules/signal_processing/macros/eqfir.bin
diff --git a/modules/signal_processing/macros/eqfir.sci b/modules/signal_processing/macros/eqfir.sci
new file mode 100755
index 000000000..d26a9132c
--- /dev/null
+++ b/modules/signal_processing/macros/eqfir.sci
@@ -0,0 +1,79 @@
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) INRIA - 1988 - C. Bunks
+//
+// 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.1-en.txt
+
+function [hn]=eqfir(nf,bedge,des,wate)
+    //<hn>=eqfir(nf,bedge,des,wate)
+    //Minimax approximation of multi-band, linear phase, FIR filter
+    //  nf    :Number of output filter points desired
+    //  bedge :Mx2 matrix giving a pair of edges for each band
+    //  des   :M-vector giving desired magnitude for each band
+    //  wate  :M-vector giving relative weight of error in each band
+    //  hn    :Output of linear-phase FIR filter coefficients
+    //!
+
+    //get number of cosines
+
+    nc=int(nf/2);
+    if nf-2*nc<>0 then,
+        flag=0;
+        nc=nc+1;
+    else,
+        flag=1;
+    end,
+
+    //make frequency grid, desired function, and weight function
+
+    [nb,c2]=size(bedge);
+    ngp=nc*16;
+    b1=bedge(:,1);
+    b2=bedge(:,2);
+    sb=sum(b2-b1);
+    delf=sb/ngp;
+    bp=round((b2-b1)/delf);
+    bsum=0;
+    for k=1:nb,
+        bpk=bp(k);
+        et=b2(k)-b1(k);
+        fg(bsum+1:bsum+bpk)=b1(k)*ones(1:bpk)+(0:bpk-1)*et/(bpk-1);
+        ds(bsum+1:bsum+bpk)=des(k)*ones(1:bpk);
+        wt(bsum+1:bsum+bpk)=wate(k)*ones(1:bpk);
+        bsum=bsum+bpk;
+    end,
+
+    //adjust values of ds and wt if filter is of even length
+
+    if flag==1 then,
+        fgs=max(size(fg));
+        if fg(fgs)>.5-%eps then,
+            fg=fg(1:fgs-1);
+            ds=ds(1:fgs-1);
+            wt=wt(1:fgs-1);
+        end,
+        cf=cos(%pi*fg);
+        ds=ds./cf;
+        wt=wt.*cf;
+    end,
+
+    //call remez
+
+    [an]=remezb(nc,fg,ds,wt);
+
+    //obtain other half of filter coefficients (by symmetry)
+
+    if flag==1 then,
+        hn(1)=.25*an(nc);
+        hn(2:nc-1)=.25*(an(nc:-1:3)+an(nc-1:-1:2));
+        hn(nc)=.5*an(1)+.25*an(2);
+        hn(nc+1:2*nc)=hn(nc:-1:1);
+    else,
+        hn=an(nc:-1:2)/2;
+        hn(nc)=an(1);
+        hn(nc+1:2*nc-1)=hn(nc-1:-1:1);
+    end,
+endfunction
diff --git a/modules/signal_processing/macros/eqiir.bin b/modules/signal_processing/macros/eqiir.bin
new file mode 100755
index 000000000..e9907a291
--- /dev/null
+++ b/modules/signal_processing/macros/eqiir.bin
diff --git a/modules/signal_processing/macros/eqiir.sci b/modules/signal_processing/macros/eqiir.sci
new file mode 100755
index 000000000..82504e33e
--- /dev/null
+++ b/modules/signal_processing/macros/eqiir.sci
@@ -0,0 +1,81 @@
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) INRIA
+//
+// 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.1-en.txt
+
+function [cells,fact,zzeros,zpoles]=eqiir(ftype,approx,om,deltap,deltas)
+    //[cells,fact,zzeros,zpoles]=eqiir(ftype,approx,om,deltap,deltas)
+    //Design of iir filter :interface with eqiir (syredi)
+    // ftype  :filter type ('lp','hp','sb','bp')
+    // approx :design approximation ('butt','cheb1','cheb2','ellip')
+    // om     :4-vector of cutoff frequencies (in radians)
+    //        :      om=<om1,om2,om3,om4>
+    //        :      0 <= om1 <= om2 <= om3 <= om4 <= pi
+    //        :When ftype=3 or 4, om3 and om4 are not used
+    //        :and may be set to 0.
+    // deltap :ripple in the passband.  0<= deltap <=1.
+    // deltas :ripple in the stopband.  0<= deltas <=1.
+    //Outputs :
+    // cells  :realization of the filter as second order cells
+    // fact   :normalization constant
+    // zzeros :zeros in the z-domain
+    // zpoles :poles in the z-domain
+    //The filter obtained is h(z)=fact*product of the elements of
+    //cells.  That is
+    //
+    //     hz=fact*prod(cells(2))./prod(cells(3))
+    //
+    //!
+    select part(approx,1);
+    case "b"
+        iapro=1
+    case "e"
+        iapro=4
+    case "c"
+        last=part(approx,length(approx));
+        if last=="1" then iapro=2,end;
+        if last=="2" then iapro=3,end
+    else
+        error(msprintf(gettext("%s: Wrong value for input argument #%d: Must be in the set {%s}.\n"),..
+        "eqiir",2,"''b[utt]'',''e[llip]'',''c[heb]1'',''c[heb]2''"))
+    end
+    select ftype;
+    case "lp"
+        ityp=1
+    case "hp"
+        ityp=2
+    case "bp"
+        ityp=3
+    case "sb"
+        ityp=4
+    else
+        error(msprintf(gettext("%s: Wrong value for input argument #%d: Must be in the set {%s}.\n"),"eqiir",1,"''lp'',''hp'',''bp'',''sb''"))
+    end
+    if max(size(om))==2 then
+        om=matrix([matrix(om,1,2),0,0],1,4),
+    end
+    [fact,b2,b1,b0,c1,c0,zzeros,zpoles]=syredi(ityp,iapro,om,deltap,deltas);
+    nb=max(size(b0));
+    coeffs=[b0;b1;b2;c0;c1];
+    coeffs=coeffs(:,1:nb);
+    coeffs=[coeffs;ones(1,nb)];
+    cells=[];
+    for col=coeffs,
+        nf=col(1:3);nd=col(4:6);
+        [n1,d1]=simp(poly(nf,"z","c"),poly(nd,"z","c"));
+        cells=[cells,syslin([],n1,d1)];
+    end
+    //crapaud...
+    if iapro==1| iapro==2  then
+        zzeros=[];
+        [k,j]=size(cells);
+        w=cells(2);
+        for k=w;
+            zzeros=[zzeros,roots(k)'];
+        end
+    end
+endfunction
diff --git a/modules/signal_processing/macros/faurre.bin b/modules/signal_processing/macros/faurre.bin
new file mode 100755
index 000000000..41724bfe8
--- /dev/null
+++ b/modules/signal_processing/macros/faurre.bin
diff --git a/modules/signal_processing/macros/faurre.sci b/modules/signal_processing/macros/faurre.sci
new file mode 100755
index 000000000..5a36c3858
--- /dev/null
+++ b/modules/signal_processing/macros/faurre.sci
@@ -0,0 +1,32 @@
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) INRIA - 1989 - G. Le Vey
+//
+// 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.1-en.txt
+
+function [P,R,T]=faurre(n,H,F,G,R0)
+    //[P,R,T]=faurre(n,H,F,G,R0)
+    //macro which computes iteratively the minimal solution of the algebraic
+    //Riccati equation and gives the matrices Rt and Tt of the filter model.
+    //   n     : number of iterations.
+    //   H,F,G : estimated triple from the covariance sequence of y.
+    //   R0    : E(yk*yk')
+    //   P    : solution of the Riccati equation after n iterations.
+    //   R,T  : gain matrix of the filter.
+    //!
+
+    //initialization
+    Pn=G*inv(R0)*G'
+    //recursion
+    for k=1:n,
+        A=(G-F*Pn*H');
+        Pn=F*Pn*F'+A/(R0-H*Pn*H')*A',
+    end;
+    P=Pn
+    //gain matrices of the filter.
+    R=R0-H*P*H';
+    T=(G-F*P*H')/R;
+endfunction
diff --git a/modules/signal_processing/macros/ffilt.bin b/modules/signal_processing/macros/ffilt.bin
new file mode 100755
index 000000000..c43d48a93
--- /dev/null
+++ b/modules/signal_processing/macros/ffilt.bin
diff --git a/modules/signal_processing/macros/ffilt.sci b/modules/signal_processing/macros/ffilt.sci
new file mode 100755
index 000000000..6bdc31a0c
--- /dev/null
+++ b/modules/signal_processing/macros/ffilt.sci
@@ -0,0 +1,74 @@
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) INRIA - 1988 - C. Bunks
+// Copyright (C) INRIA - 2002 - S Steer
+//
+// 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.1-en.txt
+
+function x=ffilt(ft,n,fl,fh)
+
+    //x=ffilt(ft,n,fl,fh)
+    //Get n coefficients of an FIR low-pass,
+    //high-pass, band-pass, or stop-band filter
+    //  ft :Filter type where ft can take the values
+    //     :       'lp' for low-pass filter
+    //     :       'hp' for high-pass filter
+    //     :       'bp' for band-pass filter
+    //     :       'sb' for stop-band filter
+    //  n  :Number of filter samples desired
+    //  fl :Low frequency cut-off
+    //  fh :High frequency cut-off
+    //     :For low and high-pass filters one cut-off
+    //     :frequency must be specified whose value is
+    //     :given in fl.  For band-pass and stop-band
+    //     :filters two cut-off frequencies must be
+    //     :specified for which the lower value is in
+    //     :fl and the higher value is in fh.
+    //  x  :Filter coefficients
+
+    //Pre-calculation
+
+    no2  = (n-1)/2;
+    ino2 = int(no2);
+
+    //Calculate n samples of the sinc function
+
+    //Low pass filter
+
+    if ft=="lp" then
+        [x]=filt_sinc(n,fl)
+    end
+
+    //High pass filter
+
+    if ft=="hp" then
+        x=filt_sinc(n,fl)
+        x=-x;
+        x(no2+1)=1+x(no2+1);
+    end
+
+    //Band pass filter
+
+    if ft=="bp" then
+        wc=%pi*(fh+fl);
+        fl=(fh-fl)/2;
+        x=filt_sinc(n,fl)
+        y=2*cos(wc*(-no2:no2));
+        x=x.*y;
+    end
+
+    //Stop band filter
+
+    if ft=="sb" then
+        wc=%pi*(fh+fl);
+        fl=(fh-fl)/2;
+        x=filt_sinc(n,fl)
+        y=2*cos(wc*(-no2:no2));
+        x=-x.*y;
+        x(no2+1)=1+x(no2+1);
+    end
+
+endfunction
diff --git a/modules/signal_processing/macros/fft2.bin b/modules/signal_processing/macros/fft2.bin
new file mode 100755
index 000000000..d8914e951
--- /dev/null
+++ b/modules/signal_processing/macros/fft2.bin
diff --git a/modules/signal_processing/macros/fft2.sci b/modules/signal_processing/macros/fft2.sci
new file mode 100755
index 000000000..c5c0ff906
--- /dev/null
+++ b/modules/signal_processing/macros/fft2.sci
@@ -0,0 +1,59 @@
+
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) INRIA - F.B
+//
+// 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.1-en.txt
+//
+function x=fft2(varargin)
+    // Two-dimensional fast Fourier transform
+    // Syntax : y = fft2(x) or y = fft2(x,m,n)
+    // Inputs :
+    // x : scalar, vector, matrix, array (real or complex)
+    // m, n : numbers (respectively) of rows and colums of x which used for the perfom of DFT,if the rows number respectively (columns) of x is more than m then x is truncated in order to have m rows, else if the rows (respectively columns) number of x is less than m then x rows are completed by 0 to have m rows.
+    //
+    // Outputs :
+    // y : scalar, vector, matrix, array (real or complex), if there is one input argument x then y and x have the same size, else if there are 3 inputs arguments then the sizes of the first and second dimension of y are equal to m and n, the others dimension sizes are equal to the size of x
+
+
+    if size(varargin) == 1 then
+        a = varargin(1);
+        if type(a) == 1 then
+            x = fft(a);
+        elseif typeof(a) == "hypermat" then
+            dims = matrix(x.dims,-1,1);
+            v = matrix(x.entries,-1,1);
+            incr = 1;
+            for k=1:2
+                dk = double(dims(k));
+                v = fft(v ,-1,dk,incr);
+                incr = incr*dk;
+            end
+            x = matrix(v,double(dims));
+        end
+    elseif size(varargin) == 3 then
+        a = varargin(1);
+        m = varargin(2);
+        n = varargin(3);
+        asize1 = size(a,1);
+        asize2 = size(a,2);
+        if type(a) == 1 then
+            x(1:min(m,asize1),1:asize2)=a(1:min(m,asize1),1:asize2);
+        elseif typeof(a) == "hypermat" then
+            dims = matrix(a.dims,-1,1);
+            dims = double(dims);
+            dims(1) = m;
+            dims(2) = n;
+            x=hypermat(dims)
+            for i=1:prod(dims(3:$))
+                x(1:min(m,asize1),1:min(n,asize2),i)=a(1:min(m,asize1),1:min(n,asize2),i);
+            end
+        end
+        x=fft2(x);
+    else
+        error(msprintf(gettext("%s: Wrong number of input arguments: %d or %d expected.\n"),"fft2",1,3));
+    end
+endfunction
diff --git a/modules/signal_processing/macros/fftshift.bin b/modules/signal_processing/macros/fftshift.bin
new file mode 100755
index 000000000..ed0c3d52f
--- /dev/null
+++ b/modules/signal_processing/macros/fftshift.bin
diff --git a/modules/signal_processing/macros/fftshift.sci b/modules/signal_processing/macros/fftshift.sci
new file mode 100755
index 000000000..b10dbc023
--- /dev/null
+++ b/modules/signal_processing/macros/fftshift.sci
@@ -0,0 +1,22 @@
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) INRIA
+//
+// 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.1-en.txt
+
+function x = fftshift(x,job)
+    if argn(2)<2 then job="all",end
+    deff("sel=fun(sk)","c=ceil(sk/2);sel=[c+1:sk,1:c]")
+    if job=="r" then job=1,elseif job=="c" then job=2,end
+    ind=list()
+    if job=="all" then
+        for sk=size(x),ind($+1)=fun(sk),end
+    else
+        for sk=size(x),ind($+1)=:,end;
+        ind(job)=fun(size(x,job))
+    end
+    x=x(ind(:))
+endfunction
diff --git a/modules/signal_processing/macros/filt_sinc.bin b/modules/signal_processing/macros/filt_sinc.bin
new file mode 100755
index 000000000..2cd2f79d9
--- /dev/null
+++ b/modules/signal_processing/macros/filt_sinc.bin
diff --git a/modules/signal_processing/macros/filt_sinc.sci b/modules/signal_processing/macros/filt_sinc.sci
new file mode 100755
index 000000000..4a40601a2
--- /dev/null
+++ b/modules/signal_processing/macros/filt_sinc.sci
@@ -0,0 +1,30 @@
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) INRIA - 1988 - C. Bunks
+//
+// 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.1-en.txt
+
+function [x]=filt_sinc(n,fl)
+
+    //x=sinc(n,fl)
+    //Calculate n samples of the function sin(2*pi*fl*t)/(pi*t)
+    //for t=-n/2:n/2 (i.e. centered around the origin).
+    //  n  :Number of samples
+    //  fl :Cut-off freq. of assoc. low-pass filter in Hertz
+    //  x  :Samples of the sinc function
+
+    no2  = (n-1)/2;
+    ino2 = int(no2);
+    wl   = fl*2*%pi;
+    xn   = sin(wl*(-no2:no2));
+    xd   = %pi*(-no2:no2);
+    if ino2==no2 then
+        xn(no2+1) = 2*fl;
+        xd(no2+1) = 1;
+    end
+    x=xn./xd;
+
+endfunction
diff --git a/modules/signal_processing/macros/filter.bin b/modules/signal_processing/macros/filter.bin
new file mode 100755
index 000000000..e78ec4238
--- /dev/null
+++ b/modules/signal_processing/macros/filter.bin
diff --git a/modules/signal_processing/macros/filter.sci b/modules/signal_processing/macros/filter.sci
new file mode 100755
index 000000000..40257a0b7
--- /dev/null
+++ b/modules/signal_processing/macros/filter.sci
@@ -0,0 +1,148 @@
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) INRIA - Serge STEER <serge.steer@inria.fr>
+//
+// 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.1-en.txt
+
+function [y, z] = filter(b, a, x, z)
+    //Implements a direct form II transposed implementation of the standard
+    //difference equation
+
+    fname = "filter"
+    [lhs, rhs] = argn(0)
+
+    if rhs < 3 | rhs > 4
+        error(msprintf(_("%s: Wrong number of input arguments: %d to %d expected.\n"), fname, 3, 4));
+    end
+
+    if rhs == 3
+        z = 0;
+    end
+
+    type_b = typeof(b);
+    type_a = typeof(a);
+    type_x = typeof(x);
+    type_z = typeof(z);
+
+    if type_b <> "constant" & type_b <> "polynomial"
+        error(msprintf(_("%s: Wrong type for input argument #%d: Real matrix or polynomial expected.\n"), fname, 1));
+    end
+
+    if type_a <> "constant" & type_a <> "polynomial"
+        error(msprintf(_("%s: Wrong type for input argument #%d: Real matrix or polynomial expected.\n"), fname, 2));
+    end
+
+    if type_x <> "constant"
+        error(msprintf(_("%s: Wrong type for input argument #%d: Real matrix expected.\n"), fname, 3));
+    end
+
+    if type_z <> "constant"
+        error(msprintf(_("%s: Wrong type for input argument #%d: Real matrix expected.\n"), fname, 4));
+    end
+
+    if ~isreal(b)
+        error(msprintf(_("%s: Wrong type for input argument #%d: Real matrix or polynomial expected.\n"), fname, 1));
+    end
+
+    if ~isreal(a)
+        error(msprintf(_("%s: Wrong type for input argument #%d: Real matrix or polynomial expected.\n"), fname, 2));
+    end
+
+    if ~isreal(x)
+        error(msprintf(_("%s: Wrong type for input argument #%d: Real matrix expected.\n"), fname, 3));
+    end
+
+    if ~isreal(z)
+        error(msprintf(_("%s: Wrong type for input argument #%d: Real matrix expected.\n"), fname, 4));
+    end
+
+    if (size(a, "c") <> 1) & (size(a, "r") <> 1)
+        error(msprintf(_("%s: Wrong size for input argument #%d: Vector expected.\n"), fname, 1));
+    end
+
+    if (size(b, "c") <> 1) & (size(b, "r") <> 1)
+        error(msprintf(_("%s: Wrong size for input argument #%d: Vector expected.\n"), fname, 2));
+    end
+
+    if (size(x, "c") <> 1) & (size(x, "r") <> 1)
+        error(msprintf(_("%s: Wrong size for input argument #%d: Vector expected.\n"), fname, 3));
+    end
+
+    if (size(z, "c") <> 1) & (size(z, "r") <> 1)
+        error(msprintf(_("%s: Wrong size for input argument #%d: Vector expected.\n"), fname, 4));
+    end
+
+    // User mixes polynomial and vector notation
+    if type_b == "polynomial" & size(a, "*") <> 1
+        error(msprintf(_("%s: Incompatible input arguments #%d and #%d: a polynomial and 1-by-1 matrix or two polynomials expected.\n"), fname, 1, 2));
+    end
+
+    // User mixes polynomial and vector notation
+    if type_a == "polynomial" & size(b, "*") <> 1
+        error(msprintf(_("%s: Incompatible input arguments #%d and #%d: a polynomial and 1-by-1 matrix or two polynomials expected.\n"), fname, 1, 2));
+    end
+
+    if type_b == "polynomial" | type_a == "polynomial"
+        c = b/a;
+        b = numer(c);
+        a = denom(c);
+        deg_b = degree(b);
+        deg_a = degree(a);
+        deg = max(deg_b, deg_a);
+        b = coeff(b, deg:-1:0);
+        a = coeff(a, deg:-1:0);
+    end
+
+    ////remove high order coefficients equal to zero
+    //i = 0; while b($ - i) == 0, i = i + 1; end;
+    //b = b(1:$ - i);
+
+    ////remove high order coefficients equal to zero
+    //i = 1; while a(i) == 0, i = i + 1; end
+    //a = a(i:$);
+
+    if a(1) == 0
+        error(msprintf(_("%s: Wrong value for input argument #%d: First element must not be %s.\n"), fname, 2, "0"));
+    end
+
+    //force vector orientation
+    b   = matrix(b, -1, 1);
+    a   = matrix(a, -1, 1);
+    mnx = size(x);
+    x   = matrix(x, 1, -1);
+
+    //normalize
+    b = b / a(1);
+    a = a / a(1);
+
+    n = max(size(b, "*"), size(a, "*"))-1;
+    if n > 0 then
+        if argn(2) < 4 then
+            z = zeros(n, 1);
+        else
+            z = matrix(z, n, 1);
+        end
+
+        //pad the numerator and denominator if necessary
+        a($ + 1:(n + 1)) = 0;
+        b($ + 1:(n + 1)) = 0;
+
+        //form state space representation
+        A     = [-a(2:$), [eye(n - 1, n - 1); zeros(1, n - 1)] ];
+        B     = b(2:$) - a(2:$) * b(1); //C = eye(1, n); D = b(1);
+
+        [z, X] = ltitr(A, B, x, z);
+        y     = X(1, :) + b(1) * x;
+
+    else
+        y     = b(1) * x;
+        z     = [];
+    end
+    //make y orientation similar to the x one
+    y = matrix(y, mnx);
+
+endfunction
+
diff --git a/modules/signal_processing/macros/find_freq.bin b/modules/signal_processing/macros/find_freq.bin
new file mode 100755
index 000000000..664879d79
--- /dev/null
+++ b/modules/signal_processing/macros/find_freq.bin
diff --git a/modules/signal_processing/macros/find_freq.sci b/modules/signal_processing/macros/find_freq.sci
new file mode 100755
index 000000000..f0fdd3417
--- /dev/null
+++ b/modules/signal_processing/macros/find_freq.sci
@@ -0,0 +1,28 @@
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) INRIA - F.D
+//
+// 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.1-en.txt
+
+function [m]=find_freq(epsilon,A,n)
+    //Search for m such that n=K(1-m1)K(m)/(K(m1)K(1-m))
+    //with m1=(epsilon*epsilon)/(A*A-1);
+    //If  m = omegar^2/omegac^2,the parameters
+    //epsilon,A,omegac,omegar and n are then
+    //compatible for defining a prototype elliptic filter.
+    //  epsilon :Passband ripple
+    //  A       :Stopband attenuation
+    //  n       :filter order
+    //  m       :Frequency needed for construction of
+    //          :elliptic filter
+    //
+    //!
+
+    m1=(epsilon*epsilon)/(A*A-1);
+    chi1=%k(1-m1)/%k(m1);
+    m=findm(chi1/n);
+
+endfunction
diff --git a/modules/signal_processing/macros/findm.bin b/modules/signal_processing/macros/findm.bin
new file mode 100755
index 000000000..b441e57ac
--- /dev/null
+++ b/modules/signal_processing/macros/findm.bin
diff --git a/modules/signal_processing/macros/findm.sci b/modules/signal_processing/macros/findm.sci
new file mode 100755
index 000000000..f5b5d0d59
--- /dev/null
+++ b/modules/signal_processing/macros/findm.sci
@@ -0,0 +1,47 @@
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) INRIA
+//
+// 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.1-en.txt
+
+function m=findm(chi)
+    //Search for m such that chi = %k(1-m)/%k(m)
+    //!
+
+    if chi < 1 then
+        t=1;
+        tn=2;
+        m=0.99999;
+        mn=2;
+        v=16*exp(-%pi/chi);
+        while abs(t-tn) > 10.*%eps
+            t=tn;
+            lln=log(16/(1-m));
+            k1=delip(1,sqrt(1-m));
+            k=delip(1,sqrt(m));
+            y=(k1*lln/%pi)-k;
+            mn=m;
+            m=1-v*exp((-%pi*y)/k1);
+            tn=m+mn;
+        end
+    else
+        t=1;
+        tn=2;
+        m=0.00001;
+        mn=0.1;
+        v=16*exp(-%pi*chi);
+        while abs(t-tn) > 10.*%eps
+            t=tn;
+            lln=log(16/m);
+            k1=delip(1,sqrt(1-m));
+            k=delip(1,sqrt(m));
+            y=(k*lln/%pi)-k1;
+            mn=m;
+            m=v*exp((-%pi*y)/k);
+            tn=m+mn;
+        end
+    end
+endfunction
diff --git a/modules/signal_processing/macros/frfit.bin b/modules/signal_processing/macros/frfit.bin
new file mode 100755
index 000000000..3223146aa
--- /dev/null
+++ b/modules/signal_processing/macros/frfit.bin
diff --git a/modules/signal_processing/macros/frfit.sci b/modules/signal_processing/macros/frfit.sci
new file mode 100755
index 000000000..86816d46a
--- /dev/null
+++ b/modules/signal_processing/macros/frfit.sci
@@ -0,0 +1,220 @@
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) INRIA
+//
+// 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.1-en.txt
+
+function [num,den]=frfit(w,fresp,order,weight)
+    //Calling sequence:
+    //[num,den]=frfit(w,fresp,r,weight)
+    //sys=frfit(w,fresp,r,weight)
+    //
+    //  w: vector of frequencies in Hz
+    //  fresp: vector of frequency responses at these frequencies.
+    //  weight: vector of weights given to each point
+    //
+    //  Fits frequency response data points by a bi-stable transfer
+    //  function
+    //                G(s) = num(s)/den(s)
+    //  of order r.
+    //  freq(num,den,%i*w) should be close to fresp
+    //
+    // changing frequencies to rad/s
+    w=2*%pi*w;
+    [LHS,RHS]=argn(0);
+    if RHS==3
+        weight=ones(w);
+    end
+    w=w(:);fresp=fresp(:);weight=weight(:);
+
+    Mean=sum(fresp)/size(fresp,"*");
+    if max(abs(fresp-Mean)) < .1*max(abs(fresp))
+        num=real(Mean); den=1; return
+    end
+
+    order1=order+1;npts=length(w);wmed=1;
+
+    if RHS < 4 then
+        weight=ones(npts,1);
+    end
+
+    M0=zeros(order+1,1); M1=M0;
+    M0(1)=1; M1(2)=1; M=[M0 M1];
+    for i=1:order-1,
+        Mt=2*[0;M1(1:order)]-M0;
+        M=[M Mt];
+        M0=M1; M1=Mt;
+    end
+
+    k=2;km=1;
+    sl0=round((log10(abs(fresp(k)))-log10(abs(fresp(km))))/(log10(w(k))-log10(w(km))));
+    k=npts;km=k-1;
+    slinf=round((log10(abs(fresp(k)))-log10(abs(fresp(km))))/(log10(w(k))-log10(w(km))));
+
+    if slinf>0&slinf<20
+        w=[w;[10;15]*w(npts)];
+        fresp=[fresp;[1;1]*10.^slinf*abs(fresp(npts))];
+        weight=[weight;1;1];
+        npts=npts+2;
+        slinf=0;
+    end
+
+
+    if sl0>0
+        mindeg=max(abs(sl0),-slinf+abs(sl0));
+    else
+        mindeg=max(abs(sl0),-slinf);
+    end
+
+    if mindeg > order
+        warning(msprintf(gettext("%s: Filter order too small.\n"),"frfit"));
+        sl0=sign(sl0)*min(order,sl0);
+        if sl0>0 then
+            slinf=-(order-abs(sl0));
+        else
+        end
+        slinf=-order;
+    end
+
+
+    jw=%i*w;mag=abs(fresp);t0=ones(npts,1); t1=jw;A=[ones(npts,1),jw];
+
+    for i=1:order-1,
+        t=2*jw.*t1-t0;
+        A=[A,t];
+        t0=t1; t1=t;
+    end
+
+    Aom=A;
+
+    //A=[A -diag(fresp)*A];AA=A;om0=w;
+    A=[A, -fresp*ones(1,size(A,2)).*A];AA=A;om0=w;
+
+    Acons=[]; ycons=[];
+    if sl0<=0
+        Acons=[Acons; zeros(-sl0,order1) M(1:-sl0,:);...
+        M(1,:) -(mag(1)*w(1)^(-sl0))*M(1-sl0,:)];
+        ycons=[ycons; zeros(-sl0+1,1)];
+    elseif sl0>0
+        Acons=[Acons; M(1:sl0,:) zeros(sl0,order1);...
+        M(1+sl0,:)*w(1)^sl0 -mag(1)*M(1,:)];
+        ycons=[ycons; zeros(sl0+1,1)];
+    end
+
+    Acons=[Acons;zeros(1,order1) M(order1,:)];ycons=[ycons;1];
+
+    if slinf<=0
+        Acons=[Acons; M(order1+slinf:order1,:) zeros(-slinf+1,order1)];
+        ycons=[ycons; mag(npts)*w(npts)^(-slinf); zeros(-slinf,1)];
+    end
+
+    [nc,nv]=size(Acons);
+
+    if nc >= nv
+        x=Acons\ycons;
+    else
+        indin=1:nv; indout=[];
+        for i=1:nc,
+            if i<nc,
+                [m,ix]=max(abs(Acons(i:nc,i:nv)));ix=ix(2);
+            else
+                [m,ix]=max(abs(Acons(i:nc,i:nv)));
+            end
+            Acons(:,i:nv)=[Acons(:,i+ix-1) Acons(:,[i:i+ix-2,i+ix:nv])];
+            indout=[indout,indin(ix)];indin=indin([1:ix-1,ix+1:nv-i+1]);
+
+            x=Acons(i:nc,i); Atmp=Acons(i:nc,i+1:nv); ytmp=ycons(i:nc);
+            if x(1)>=0
+                sgn=1;
+            else
+                sgn=-1;
+            end
+            aux=sgn*sqrt(x'*x);x(1)=x(1)+aux;nx2=.5*x'*x;
+            Atmp=[[-aux;zeros(nc-i,1)] Atmp-x*((x'*Atmp)/nx2)];
+            ytmp=ytmp-x*((x'*ytmp)/nx2);
+            Acons(i:nc,i:nv)=Atmp;ycons(i:nc)=ytmp;
+        end
+        perm=[indout indin];
+        Ac1=Acons(:,1:nc); Ac2=Acons(:,nc+1:nv);
+        A=A(1:npts,:);
+        A=[real(A);imag(A)];
+        A=A(:,perm);
+        A1=A(:,1:nc); A2=A(:,nc+1:nv);
+        A=A2-A1*(Ac1\Ac2);
+        y=-A1*(Ac1\ycons);
+
+        fweight=ones(npts,1);
+        ind=find(w > 10);
+
+        fweight(ind)=((1) ./(w(ind).^order));
+        ind=find(w < .01);
+
+        fweight(ind)=(1./(w(ind)^min(0,sl0)))';
+        fweight=weight.*fweight;
+        //Wt=diag([fweight;fweight]);
+        //x=pinv(Wt*A)*(Wt*y);
+        //x=pinv([fweight;fweight]*ones(1,size(A,2)).*A)*(vvv.*y);
+        vvv=[fweight;fweight];
+        x=(vvv*ones(1,size(A,2)).*A)\(vvv.*y);
+        x=[Ac1\(ycons-Ac2*x);x];
+        [s,perm]=gsort(-perm);
+        s=-s;
+        x=x(perm);
+
+        nresp=Aom*x(1:order1); dresp=Aom*x(order1+1:2*order1);
+        relerr=abs(abs(nresp./dresp)./mag-1)';
+        for kk=1:prod(size(relerr));relerr(kk)=min(1,relerr(kk));end
+        relerr=relerr(:);
+        ind=find(mag < .01); relerr(ind)=relerr(ind)*.3;
+        ind=find(relerr>.5);
+        if ind<>[]
+            fweight(ind)=( exp(relerr(ind)*log(10)) ).*fweight(ind);
+        end
+        //Wt=diag([fweight;fweight]);
+        //x=pinv(Wt*A)*(Wt*y);
+        //x=pinv([fweight;fweight]*ones(1,size(A,2)).*A)*(vvv.*y);
+        x=([fweight;fweight]*ones(1,size(A,2)).*A)\(vvv.*y);
+        x=[Ac1\(ycons-Ac2*x);x];x=x(perm);
+    end
+
+    a=x(1:order1); b=x(order1+1:2*order1);
+    num=fliplr((M*a)');den=fliplr((M*b)');
+
+    err=abs(AA*x);
+    num=num(-slinf+1:order1);
+
+    junk=poly(fliplr(num),"s","c")
+    rn=roots(junk); rn=-abs(real(rn))+%i*imag(rn);
+    l=length(rn);
+    for k=1:l
+        if real(rn(k))>-1e-3
+            rn(k)=-1.e-3+%i*imag(rn(k))
+        end
+    end
+
+    junk=poly(fliplr(den),"s","c");
+    dn=roots(junk); dn=-abs(real(dn))+%i*imag(dn);
+    l=length(dn);
+    for k=1:l
+        if real(dn(k))>-1e-3
+            dn(k)=-1.e-3+%i*imag(dn(k))
+        end
+    end
+
+    polyrn=poly(fliplr(rn),"s");
+    polydn=poly(fliplr(dn),"s");
+    num=real(num(1)*polyrn);
+    den=real(den(1)*polydn);
+    if LHS==1
+        num=syslin("c",num/den)
+    end
+
+endfunction
+function m=fliplr(m)
+    //Utility fct
+    [p,q]=size(m);
+    m=m(:,q:-1:1);
+endfunction
diff --git a/modules/signal_processing/macros/frmag.bin b/modules/signal_processing/macros/frmag.bin
new file mode 100755
index 000000000..0e01c7eff
--- /dev/null
+++ b/modules/signal_processing/macros/frmag.bin
diff --git a/modules/signal_processing/macros/frmag.sci b/modules/signal_processing/macros/frmag.sci
new file mode 100755
index 000000000..4ace029b8
--- /dev/null
+++ b/modules/signal_processing/macros/frmag.sci
@@ -0,0 +1,79 @@
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) INRIA - 1988 - C. Bunks
+//
+// 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.1-en.txt
+
+function [xm,fr]=frmag(num,den,npts)
+    //[xm,fr]=frmag(num[,den],npts)
+    //Calculates the magnitude of the frequency respones of
+    //FIR and IIR filters.  The filter description can be
+    //one or two vectors of coefficients, one or two polynomials,
+    //or a rational polynomial.
+    //Case 1 (When den is not given):
+    //  num  :Vector coefficients/Polynomial/Rational
+    //       :polynomial of filter
+    //Case 2 (When den is given):
+    //  num  :Vector coefficients/Polynomial of filter numerator
+    //  den  :Vector coefficients/Polynomial of filter denominator
+    //Case 1 and 2:
+    //  npts :Number of points desired in frequency response
+    //  xm   :Magnitude of frequency response at the points fr
+    //  fr   :Points in the frequency domain where
+    //       :magnitude is evaluated
+    //!
+
+    select argn(2)
+    case 2 then //frmag(sys,npts)
+        npts=den;
+        if typeof(num)=="rational" then
+            h=num;
+            if size(h,"*")<>1 then
+                error(msprintf(_("%s: Wrong size for input argument #%d: A single input, single output system expected.\n"),"frmag",1))
+            end
+            num=h.num;den=h.den;
+        elseif typeof(num)=="state-space" then
+            h=ss2tf(num);
+            if size(h,"*")<>1 then
+                error(msprintf(_("%s: Wrong size for input argument #%d: A single input, single output system expected.\n"),"frmag",1))
+            end
+            num=h.num;den=h.den;
+        elseif type(num)==2 then,
+            if size(num,"*")<>1 then
+                error(msprintf(_("%s: Wrong size for input argument #%d: A polynomial expected.\n"),"frmag",1));
+            end
+            den=poly(1,"z","c");
+        elseif type(num)==1 then,
+            num=poly(num,"z","c");
+            den=1
+        else
+            error(msprintf(_("%s: Wrong type for input argument #%d:  Linear dynamical system or row vector of floats expected.\n"),"frmag",1))
+        end,
+    case 3 then,//frmag(num,den,npts)
+        if type(num)==2 then,
+            if size(num,"*")<>1 then
+                error(msprintf(_("%s: Wrong size for input argument #%d: A polynomial expected.\n"),"frmag",1));
+            end
+        elseif type(num)==1 then,
+            num=poly(num,"z","c");
+        else,
+            error(msprintf(_("%s: Wrong size for input argument #%d: A polynomial expected.\n"),"frmag",1));
+        end,
+        if type(den)==2 then,
+            if size(den,"*")<>1 then
+                error(msprintf(_("%s: Wrong size for input argument #%d: A polynomial expected.\n"),"frmag",2));
+            end
+        elseif type(den)==1 then,
+            den=poly(den,"z","c");
+        else,
+            error(msprintf(_("%s: Wrong size for input argument #%d: A polynomial expected.\n"),"frmag",2));
+        end,
+    end
+    fr=linspace(0,1/2,npts+1);
+    fr($)=[];
+    dfr=exp(2*%i*%pi*fr);
+    xm=abs(freq(num,den,dfr));
+endfunction
diff --git a/modules/signal_processing/macros/fsfirlin.bin b/modules/signal_processing/macros/fsfirlin.bin
new file mode 100755
index 000000000..8afb05d06
--- /dev/null
+++ b/modules/signal_processing/macros/fsfirlin.bin
diff --git a/modules/signal_processing/macros/fsfirlin.sci b/modules/signal_processing/macros/fsfirlin.sci
new file mode 100755
index 000000000..bc2936478
--- /dev/null
+++ b/modules/signal_processing/macros/fsfirlin.sci
@@ -0,0 +1,28 @@
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) INRIA - 1989 - G. Le Vey
+//
+// 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.1-en.txt
+function [hst]=fsfirlin(hd,flag)
+    //<hst>=fsfirlin(hd,flag)
+    //macro for the design of FIR, linear phase filters
+    //using the frequency sampling technique
+    //  hd   : vector of desired frequency response samples
+    //  flag : is equal to 1 or 2,
+    //         according to the choice of type 1 or type 2 design
+    //  hst  : vector giving the approximated continuous response
+    //         on a dense grid of frequencies
+    //!
+
+
+    n1=prod(size(hd));//filter length
+    if int(n1/2)==n1/2,n=2*n1;else,n=2*n1+1;end;//even or odd length
+    scd=sincd(n,flag);//calculates the function Sin(N*x)/Sin(x)
+    hst=hd(1)*scd(4*n+1:6*n+1);
+    eps=(-1)**(n-1);
+    for j=1:n1-1,hst=hst+hd(j+1)*[scd(-4*j+4*n+1:-4*j+6*n+1)+..
+    eps*scd(4*j+1:4*j+2*n+1)];end;
+endfunction
diff --git a/modules/signal_processing/macros/group.bin b/modules/signal_processing/macros/group.bin
new file mode 100755
index 000000000..ca36c1920
--- /dev/null
+++ b/modules/signal_processing/macros/group.bin
diff --git a/modules/signal_processing/macros/group.sci b/modules/signal_processing/macros/group.sci
new file mode 100755
index 000000000..ac06aef8b
--- /dev/null
+++ b/modules/signal_processing/macros/group.sci
@@ -0,0 +1,140 @@
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) INRIA - 1988 - C. Bunks
+//
+// 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.1-en.txt
+
+function [tg,fr]=group(npts,a1i,a2i,b1i,b2i)
+    //Calculate the group delay of a digital filter
+    //with transfer function h(z).
+    //The filter specification can be in coefficient form,
+    //polynomial form, rational polynomial form, cascade
+    //polynomial form, or in coefficient polynomial form.
+    //  npts :Number of points desired in calculation of group delay
+    //  a1i  :In coefficient, polynomial, rational polynomial, or
+    //       :cascade polynomial form this variable is the transfer
+    //       :function of the filter.  In coefficient polynomial
+    //       :form this is a vector of coefficients (see below).
+    //  a2i  :In coeff poly form this is a vector of coeffs
+    //  b1i  :In coeff poly form this is a vector of coeffs
+    //  b2i  :In coeff poly form this is a vector of coeffs
+    //  tg   :Values of group delay evaluated on the grid fr
+    //  fr   :Grid of frequency values where group delay is evaluated
+    //
+    //In the coefficient polynomial form the tranfer function is
+    //formulated by the following expression:
+    //
+    //       h(z)=prod(a1i+a2i*z+z**2)/prod(b1i+b2i*z+z^2)
+    //
+    //!
+
+
+    //get frequency grid values in [0,.5)
+
+    fr=(0:.5/npts:.5*(npts-1)/npts);
+
+    //get the of arguments used to called group
+
+    [ns,ne]=argn(0);
+    if and(ne <> [2 5]) then
+        error(sprintf(_("%s: Wrong number of input argument(s): %d or %d expected.\n"), "group", 2, 5));
+    end
+
+    //if the number of arguments is 2 then
+    //the case may be non-cascade
+
+    hcs=1;
+    if ne==2 then
+
+        //get the type of h and the variable name
+
+        h=a1i;
+        ht=type(h);
+        if and(ht <> [1 2 15 16]) then
+            error(sprintf(_("%s: Wrong type for input argument #%d: A real or polynomial matrix or a rational expected.\n"), "group", 2));
+        end
+        if ht == 16 & or(h.dt == "c" | h.dt == []) then
+            error(sprintf(_("%s: Wrong type for input argument #%d: A discrete system expected.\n"), "group", 2));
+        end
+
+        //if ht==1 then h is a vector containing filter coefficients
+
+        if ht==1 then
+
+            //make h a rational polynomial
+
+            hcs=max(size(h));
+            z=poly(0,"z");
+            h=poly(h,"z","c");
+            h=gtild(h,"d")*(1/z^(hcs-1));
+        elseif ht == 2 then
+            z = poly(0, varn(h));
+        else  //ht==15|ht==16
+            z=poly(0, varn(h(3)));
+            hcs=max(size(h(2)));
+        end,
+
+        //if the rational polynomial is not in cascade form then
+
+        if hcs==1 then
+            //get the derivative of h(z)
+
+            hzd=derivat(h);
+
+            //get the group delay of h(z)
+            tgz=-z*hzd/h;
+
+            //evaluate tg
+
+            rfr=exp(2*%pi*%i*fr);
+
+
+            tg=real(freq(tgz(2),tgz(3),rfr));
+
+            //done with non-cascade calculation of group delay
+
+        end,
+        //re-organize if h is in cascade form
+
+        if hcs>1 then
+            xc=[coeff(h(2)),coeff(h(3))];
+            a2i=xc(1:hcs);
+            a1i=xc(hcs+1:2*hcs);
+            b2i=xc(3*hcs+1:4*hcs);
+            b1i=xc(4*hcs+1:5*hcs);
+            ne=5;
+        end,
+    end,
+
+    //if implementation is in cascade form
+
+    if ne==5 then
+
+        //a1i,a2i,b1i,b2i are the coefficients of a
+        //second order decomposition of the filter
+        //(i.e. in cascade form, see Deczky)
+
+        phi=2*%pi*fr;
+        z=poly(0,"z");
+        ejw=freq(z,1,exp(%i*phi));
+        emjw=freq(z,1,exp(-%i*phi));
+
+        a1=a1i'*ones(phi);
+        b1=b1i'*ones(phi);
+        a2=a2i'*ones(phi);
+        b2=b2i'*ones(phi);
+        ejw=ones(a1i)'*ejw;
+        emjw=ones(a1i)'*emjw;
+
+        aterm=(a1+2*ejw)./(a1+ejw+a2.*emjw);
+        bterm=(b1+2*ejw)./(b1+ejw+b2.*emjw);
+
+        tgi=real(bterm-aterm);
+        tg=ones(a1i)*tgi;
+
+        //done with cascade calculation of group delay
+    end
+endfunction
diff --git a/modules/signal_processing/macros/hank.bin b/modules/signal_processing/macros/hank.bin
new file mode 100755
index 000000000..ea246ff42
--- /dev/null
+++ b/modules/signal_processing/macros/hank.bin
diff --git a/modules/signal_processing/macros/hank.sci b/modules/signal_processing/macros/hank.sci
new file mode 100755
index 000000000..9e6dc075e
--- /dev/null
+++ b/modules/signal_processing/macros/hank.sci
@@ -0,0 +1,65 @@
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) INRIA - 1989 - G. Le Vey
+// Copyright (C) Scilab Enterprises - 2011 - Calixte DENIZET
+// Copyright (C) DIGITEO - 2012 - Allan CORNET
+//
+// 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.1-en.txt
+
+function hk = hank(m, n, seq_cov)
+    //hk = hank(m,n,cov)
+    //this macro builds the hankel matrix of size (m*d,n*d)
+    //from the covariance sequence of a vector process
+    //   m  : number of bloc-rows
+    //   n  : number of bloc-columns
+    //   seq_cov: sequence of covariances; it must be given as :[R0 R1 R2...Rk]
+    //   hk : computed hankel matrix
+    //
+
+    [lhs, rhs] = argn(0);
+    if rhs <> 3 then
+        error(msprintf(gettext("%s: Wrong number of input argument(s): %d expected.\n"), "hank", 3));
+    end
+
+    if and(type(m) <> [1 8]) then
+        error(msprintf(gettext("%s: Wrong type for input argument #%d: A positive integer expected.\n"), "hank", 1));
+    end
+
+    if (size(m, "*") <> 1) then
+        error(msprintf(gettext("%s: Wrong size for input argument #%d: A positive integer expected.\n"), "hank", 1));
+    end
+
+    if m <= 0 | m <> fix(m) then
+        error(msprintf(gettext("%s: Wrong value for input argument #%d: A positive integer expected.\n"), "hank", 1));
+    end
+
+    if and(type(n) <> [1 8]) then
+        error(msprintf(gettext("%s: Wrong type for input argument #%d: A positive integer expected.\n"), "hank", 2));
+    end
+
+    if (size(n, "*") <> 1) then
+        error(msprintf(gettext("%s: Wrong size for input argument #%d: A positive integer expected.\n"), "hank", 2));
+    end
+
+    if n <= 0 | n <> fix(n) then
+        error(msprintf(gettext("%s: Wrong value for input argument #%d: A positive integer expected.\n"), "hank", 2));
+    end
+
+    t = type(seq_cov);
+    if (t > 10) | (t <= 0) then
+        error(msprintf(gettext("%s: Wrong values for input argument #%d: Unsupported ''%s'' type.\n"), "hank", 3, typeof(t)));
+    end
+
+    [r, c] = size(seq_cov);
+    mr = m * r;
+    nr = n * r;
+    if (mr + nr - r > c) then
+        error(msprintf(gettext("%s: Wrong size for input arguments: Incompatible sizes.\n"), "hank"));
+    end
+
+    index = ones(1, nr) .*. (1:r:mr)' + (0:(nr - 1)) .*. ones(m, 1);
+    hk = matrix(seq_cov(:, index), mr, -1);
+endfunction
diff --git a/modules/signal_processing/macros/hilb.bin b/modules/signal_processing/macros/hilb.bin
new file mode 100755
index 000000000..0a1c053c3
--- /dev/null
+++ b/modules/signal_processing/macros/hilb.bin
diff --git a/modules/signal_processing/macros/hilb.sci b/modules/signal_processing/macros/hilb.sci
new file mode 100755
index 000000000..c4f1dc8f3
--- /dev/null
+++ b/modules/signal_processing/macros/hilb.sci
@@ -0,0 +1,61 @@
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) INRIA
+//
+// 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.1-en.txt
+
+function xh=hilb(n,wtype,par)
+    //xh=hilb(n [,wtype [,par]])
+    //returns the first n points of the
+    //Hilbert IIR filter  centered around the origin.
+    //That is, xh=(2/(n*pi))*(sin(n*pi/2))**2.
+    //Window type and window parameters are optional.
+    // n     :Number of points in filter
+    // wtype :window type ('re','tr','hn','hm','kr','ch')
+    //       :     default wtype='re'
+    // par   :window parameter for wtype='kr' or 'ch'
+    //       :     default par=[0 0]
+    //       :see the macro window for more help
+    // xh    :Hilbert transform
+    //
+    //!
+    // References:
+    // http://ieeexplore.ieee.org/iel4/78/7823/00330385.pdf?tp=&arnumber=330385&isnumber=7823
+    //  A. Reilly, G. Frazer, and B. Boashash, "Analytic signal generation
+    //  Tips and traps,¡ IEEE Trans. Signal Processing, vol. 42,
+    //  pp.3241-3245, Nov. 1994.
+    [lhs,rhs]=argn(0);
+    if rhs==1 then,
+        wtype="re";
+        par=[0 0];
+    elseif rhs==2 then,
+        par=[0 0];
+    end,
+
+    if type(n)<>1|size(n,"*")<>1 then
+        error(msprintf(gettext("%s: Wrong size for input argument #%d: A scalar expected.\n"),"hilb",1));
+    end
+    if (n - floor(n) <> 0) | (floor(n/2) == n/2) then
+        error(msprintf(gettext("%s: Wrong type for input argument #%d: An odd integer expected.\n"),"hilb",1));
+    end
+    if and(wtype<>["re","tr","hn","hm","kr","ch"]) then
+        error(msprintf(gettext("%s: Wrong value for input argument #%d: Must be in the set {%s}.\n") ,"hilb",2,"''re'',''tr'',''hn'',''hm'',''kr'',''ch''"));
+
+    end
+    no2=int(n/2)
+    if no2==n/2 then,
+        xh0=[]
+    else
+        xh0=0
+    end,
+    th=zeros(1,no2);
+    th(1:2:no2)=ones(1:2:no2)./(1:2:no2);
+    xh=[-th(no2:-1:1),xh0,th];
+    xh=2*xh/%pi;
+
+    [win_l,cwp]=window(wtype,n,par);
+    xh=xh.*win_l;
+endfunction
diff --git a/modules/signal_processing/macros/hilbert.bin b/modules/signal_processing/macros/hilbert.bin
new file mode 100755
index 000000000..1345ca9fb
--- /dev/null
+++ b/modules/signal_processing/macros/hilbert.bin
diff --git a/modules/signal_processing/macros/hilbert.sci b/modules/signal_processing/macros/hilbert.sci
new file mode 100755
index 000000000..1a838dacf
--- /dev/null
+++ b/modules/signal_processing/macros/hilbert.sci
@@ -0,0 +1,46 @@
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2008 - INRIA
+//
+// 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.1-en.txt
+//
+//
+
+function x = hilbert(xr)
+    // Marple, S.L., "Computing the discrete-time analytic signal via FFT,"
+    // IEEE Transactions on Signal Processing, Vol. 47, No.9 (September
+    // 1999), pp.2600-2603
+    // http://ieeexplore.ieee.org/iel5/78/16975/00782222.pdf?arnumber=782222
+
+    if  type(xr)<>1 then
+        error(msprintf(gettext("%s: Wrong type for input argument #%d: Array of floating point numbers expected.\n"),"hilbert",1))
+    end
+
+    n=size(xr,"*");
+
+    if n==0 then x=[],return,end
+    if ~isreal(xr,0) then
+        error(msprintf(gettext("%s: Input argument #%d must be real.\n"),"hilbert",1));
+    end
+
+    no2 = int(n/2);
+
+    x = fft(real(xr),-1);
+
+    h  = zeros(xr);
+
+    if ((2*no2) == n) then  // n is even
+        h([1,no2+1]) = 1;
+        h(2:no2) = 2;
+    else // n is odd
+        h(1) = 1;
+        h(2:(n+1)/2) = 2;
+    end
+
+    x = fft(x.*h,1);
+
+endfunction
+
diff --git a/modules/signal_processing/macros/idct.bin b/modules/signal_processing/macros/idct.bin
new file mode 100755
index 000000000..1ee8ec790
--- /dev/null
+++ b/modules/signal_processing/macros/idct.bin
diff --git a/modules/signal_processing/macros/idct.sci b/modules/signal_processing/macros/idct.sci
new file mode 100755
index 000000000..70aa3ae15
--- /dev/null
+++ b/modules/signal_processing/macros/idct.sci
@@ -0,0 +1,15 @@
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) INRIA - 2012 - S. Steer
+//
+// 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.1-en.txt
+function x=idct(a,varargin)
+    if type(a)==1|(typeof(a)=="hypermat"&type(a.entries)==1) then
+        x=dct(a,1,varargin(:))
+    else
+        error(msprintf(_("%s: Wrong type for input argument #%d: Array of double expected.\n"),"idct", 1));
+    end
+endfunction
diff --git a/modules/signal_processing/macros/idst.bin b/modules/signal_processing/macros/idst.bin
new file mode 100755
index 000000000..9adab5074
--- /dev/null
+++ b/modules/signal_processing/macros/idst.bin
diff --git a/modules/signal_processing/macros/idst.sci b/modules/signal_processing/macros/idst.sci
new file mode 100755
index 000000000..7045c8147
--- /dev/null
+++ b/modules/signal_processing/macros/idst.sci
@@ -0,0 +1,15 @@
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) INRIA - 2012 - S. Steer
+//
+// 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.1-en.txt
+function x=idst(a,varargin)
+    if type(a)==1|(typeof(a)=="hypermat"&type(a.entries)==1) then
+        x=dst(a,1,varargin(:))
+    else
+        error(msprintf(_("%s: Wrong type for input argument #%d: Array of double expected.\n"),"idst", 1));
+    end
+endfunction
diff --git a/modules/signal_processing/macros/ifft.bin b/modules/signal_processing/macros/ifft.bin
new file mode 100755
index 000000000..7b04d217d
--- /dev/null
+++ b/modules/signal_processing/macros/ifft.bin
diff --git a/modules/signal_processing/macros/ifft.sci b/modules/signal_processing/macros/ifft.sci
new file mode 100755
index 000000000..8eb491892
--- /dev/null
+++ b/modules/signal_processing/macros/ifft.sci
@@ -0,0 +1,15 @@
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) INRIA - 1988 - C. Bunks
+//
+// 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.1-en.txt
+function x=ifft(a,varargin)
+    if type(a)==1|(typeof(a)=="hypermat"&type(a.entries)==1) then
+        x=fft(a,1,varargin(:))
+    else
+        error(msprintf(_("%s: Wrong type for input argument #%d: Array of double expected.\n"),"ifft", 1));
+    end
+endfunction
diff --git a/modules/signal_processing/macros/ifftshift.bin b/modules/signal_processing/macros/ifftshift.bin
new file mode 100755
index 000000000..8bbf09968
--- /dev/null
+++ b/modules/signal_processing/macros/ifftshift.bin
diff --git a/modules/signal_processing/macros/ifftshift.sci b/modules/signal_processing/macros/ifftshift.sci
new file mode 100755
index 000000000..8d030ac52
--- /dev/null
+++ b/modules/signal_processing/macros/ifftshift.sci
@@ -0,0 +1,27 @@
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2011 - Scilab Enterprises - Allan Cornet
+//
+// 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.1-en.txt
+
+function y = ifftshift(x)
+
+    if argn(2) < 1 then
+        error(sprintf(_("%s: Wrong number of input argument(s): %d expected.\n"), "ifftshift", 1));
+    end
+
+    numDims = ndims(x);
+    idx = list(1, numDims);
+
+    for k = 1:numDims
+        m = size(x, k);
+        p = floor(m/2);
+        idx(k) = [p+1:m 1:p];
+    end
+
+    y = x(idx(:));
+
+endfunction
diff --git a/modules/signal_processing/macros/iir.bin b/modules/signal_processing/macros/iir.bin
new file mode 100755
index 000000000..7e689ce74
--- /dev/null
+++ b/modules/signal_processing/macros/iir.bin
diff --git a/modules/signal_processing/macros/iir.sci b/modules/signal_processing/macros/iir.sci
new file mode 100755
index 000000000..ce549bd06
--- /dev/null
+++ b/modules/signal_processing/macros/iir.sci
@@ -0,0 +1,81 @@
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) INRIA - 1988 - C. Bunks
+//
+// 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.1-en.txt
+
+function [hz,zz,gz]=iir(n,ftype,fdesign,frq,delta)
+    // hz=iir(n,ftype,fdesign,frq,delta)
+    //macro which designs an iir digital filter
+    //using analog filter designs.
+    //the arguments of the filter are:
+    //  n        :filter order (pos. integer)
+    //  ftype    :specification of the filter type
+    //           :   ftype=('lp','hp','bp','sb')
+    //  fdesign  :specifiation of the analog filter design
+    //           :   fdesign=('butt','cheb1','cheb2','ellip')
+    //  frq      :2-vector of discrete cut-off frequencies
+    //           :(i.e., 0<frq<.5). For lp and hp filters only
+    //           :frq(1) is used.  For bp and sb filters frq(1)
+    //           :is the lower cut-off frequency and frq(2) is
+    //           :the upper cut-off frequency
+    //  delta    :2-vector of error values for cheb1, cheb2, and
+    //           :ellip filters where only delta(1) is used for
+    //           :cheb1 case, only delta(2) is used for cheb2 case,
+    //           :and delta(1) and delta(2) are both used for
+    //           :ellip case.
+    //           :          0<delta(1),delta(2)<1
+    //           :for cheb1 filters:  1-delta(1)<ripple<1 in passband
+    //           :for cheb2 filters:  0<ripple<delta(2)   in stopband
+    //           :for ellip filters:  1-delta(1)<ripple<1 in passband
+    //           :                    0<ripple<delta(2)   in stopband
+    //
+    //!
+    if and(argn(1)<>[1 3]) then
+        error(msprintf(gettext("%s: Wrong number of output arguments: %d or %d expected.\n"),"iir",1,3))
+    end
+    //select analog filter design for low-pass filter with fc=.25
+    if type(n)<>1 then
+        error(msprintf(gettext("%s: Wrong type for input argument #%d: Array of floating point numbers expected.\n"),"iir",1))
+    end
+    if  size(n,"*")<>1 then
+        error(msprintf(gettext("%s: Wrong size for input argument #%d: A scalar expected.\n"),"iir",1))
+    end
+    if  n<0|n<>round(n) then
+        error(msprintf(gettext("%s: Wrong values for input argument #%d: Non-negative integers expected.\n"),"iir",1))
+    end
+
+    if and(ftype<>["lp","hp","bp","sb"]) then
+        error(msprintf(gettext("%s: Wrong value for input argument #%d: Must be in the set {%s}.\n"),"iir",2,"''lp'',''hp'',''bp'',''sb''"))
+    end
+    if and(fdesign<>["butt","cheb1","cheb2","ellip"]) then
+        error(msprintf(gettext("%s: Wrong value for input argument #%d: Must be in the set {%s}.\n"),"iir",2,"''butt'',''cheb1'',''cheb2'',''ellip''"))
+    end
+    if type(frq)<>1 then
+        error(msprintf(gettext("%s: Wrong type for input argument #%d: Array of floating point numbers expected.\n"),"iir",4))
+    end
+    if type(delta)<>1 then
+        error(msprintf(gettext("%s: Wrong type for input argument #%d: Array of floating point numbers expected.\n"),"iir",5))
+    end
+
+    if max(abs(frq))>0.5 then
+        error(msprintf(gettext("%s: Wrong values for input argument #%d: Elements must be in the interval [%s, %s].\n"),"iir",4,"0","0.5"));
+    end
+    if delta(1)<0|delta(2)>1 then
+        error(msprintf(gettext("%s: Wrong values for input argument #%d: Elements must be in the interval [%s, %s].\n"),"iir",4,"0","1"));
+    end
+
+    [hs,pc,zc,gc]=analpf(n,fdesign,delta,2);
+    //make digital low-pass filter from analog low-pass filter
+    z=poly(0,"z");[pd,zd,gd]=bilt(pc,zc,gc,2*(z-1),(z+1));
+    //do change of variables to obtain general digital filter
+    if argn(1)==1 then
+        hz=trans(pd,zd,gd,ftype,frq);
+    else
+        [pz,zz,gz]=trans(pd,zd,gd,ftype,frq);
+        hz=pz
+    end
+endfunction
diff --git a/modules/signal_processing/macros/iirgroup.bin b/modules/signal_processing/macros/iirgroup.bin
new file mode 100755
index 000000000..b67ee0eb9
--- /dev/null
+++ b/modules/signal_processing/macros/iirgroup.bin
diff --git a/modules/signal_processing/macros/iirgroup.sci b/modules/signal_processing/macros/iirgroup.sci
new file mode 100755
index 000000000..dfe64a74f
--- /dev/null
+++ b/modules/signal_processing/macros/iirgroup.sci
@@ -0,0 +1,52 @@
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) INRIA
+//
+// 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.1-en.txt
+
+function [lt,grad]=iirgroup(p,r,theta,omega,wt,td)
+    //    p===> Lp
+    //    r=magnitude of poles and zeros of filters
+    //theta=phase   "    "   "    "    "
+    //omega=frequencies corresponding to filters specs.
+    //wt=weighting function for group delay
+    //!
+    [m,n]=size(td);if m>n,td=td';end;
+    [m,n]=size(omega);if m>n,omega=omega';end;
+    [m,n]=size(r);if n>m,r=r';m=n;end;
+    [m,n]=size(theta);if n>m,theta=theta';m=n;end;
+    m=m/2;
+    //
+    //VARIABLES LOCALES
+    unr=ones(r);unom=ones(omega);un=unr(1:m,:).*.unom;
+    n1=r.*.unom;n2=unr.*.omega;n3=theta.*.unom;n4=(unr+r.*r).*.unom;
+    cw=cos(omega);sw=sin(omega);ct=cos(theta);st=sin(theta);
+    c=cw.*.ct;s=sw.*.st;cos1=c+s;cos2=c-s;
+    c=sw.*.ct;s=cw.*.st;sin1=c-s;sin2=c+s;
+    //
+    //FORMES QUADRATIQUES
+    k1=n4-2*n1.*cos1;k2=n4-2*n1*cos2;
+    //
+    //RETARD DE GROUPE
+    grp=[];k3=k1(1:m,:);k4=k1(m+1:2*m,:);k5=k2(1:m,:);k6=k2(m+1:2*m,:);
+    c1=cos1(1:m,:);c2=cos1(m+1:2*m,:);c3=cos2(1:m,:);c4=cos2(m+1:2*m,:);
+    r1=n1(1:m,:);r2=n1(m+1:2*m,:);
+    t=(un-r2.*c2)./k4+(un-r2.*c4)./k6-(un-r1.*c1)./k3-(un-r1.*c3)./k5;
+    for k=t,grp=[grp sum(k)];end;
+    //
+    //GRADIENT DU RETARD DE GROUPE
+    k1=k1.*k1;k2=k2.*k2;
+    grrgrp=(n4.*cos1-2*n1)./k1+(n4.*cos2-2*n1)./k2;
+    grrgrp(1:m,:)=-grrgrp(1:m,:);
+    n5=n1.*(-n4+2*unr.*.unom);
+    grtgrp=(n5.*sin1)./k1-(n5.*sin2)./k2;
+    grtgrp(1:m,:)=-grtgrp(1:m,:);
+    //
+    //CRITERE D'ERREUR EN LE RETARD DE GROUPE ET SON GRADIENT
+    t=grp-td;t1=t^(2*p);t1=t1.*wt;lt=sum(t1);
+    t1=(t1./t)*2*p;
+    grad=[grrgrp*t1' grtgrp*t1'];
+endfunction
diff --git a/modules/signal_processing/macros/iirlp.bin b/modules/signal_processing/macros/iirlp.bin
new file mode 100755
index 000000000..74c69041c
--- /dev/null
+++ b/modules/signal_processing/macros/iirlp.bin
diff --git a/modules/signal_processing/macros/iirlp.sci b/modules/signal_processing/macros/iirlp.sci
new file mode 100755
index 000000000..5b6f04f9c
--- /dev/null
+++ b/modules/signal_processing/macros/iirlp.sci
@@ -0,0 +1,49 @@
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) INRIA
+//
+// 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.1-en.txt
+
+function [cout,grad,ind]=iirlp(x,ind,p,flag,lambda,omega,ad,wa,td,wt)
+    //
+    //optimization of IIR filters IIR with  Lp criterium for magnitude
+    //                                      and/or group delay
+    //     -cf Rabiner & Gold pp270-273-
+    //
+    //auteur : G. Le Vey
+    //
+    //   p ===> critere Lp
+    //
+    //r=module des poles et zeros des filtres
+    //theta=argument des  "    "   "    "    "
+    //omega=frequences ou sont donnees les specifications des filtres
+    //wa,wt=fonctions de ponderation pour l'amplitude et le
+    //retard de groupe ad,td=amplitudes et retard de groupe desires
+    //!
+    r=x(:,1);theta=x(:,2);
+    [m,n]=size(ad);if m>n,ad=ad';end
+    [m,n]=size(td);if m>n,td=td';end
+    [m,n]=size(omega);if m>n,omega=omega';end;
+    [m,n]=size(r);if n>m,r=r';m=n;end;
+    [m,n]=size(theta);if n>m,theta=theta';m=n;end;
+    //
+    select flag
+    case "a"
+        //AMPLITUDE
+        [cout,grad]=iirmod(p,r,theta,omega,wa,ad);
+        //
+    case "gd"
+        //RETARD DE GROUPE
+        [cout,grad]=iirgroup(p,r,theta,omega,wt,td);
+        //
+    else
+        //AMPLITUDE ET RETARD DE GROUPE
+        [la,ga]=iirmod(p,r,theta,omega,wa,ad);
+        [lt,gt]=iirgroup(p,r,theta,omega,wt,td);
+        cout=lambda*la+(1-lambda)*lt;
+        grad=lambda*ga+(1-lambda)*gt;
+    end;
+endfunction
diff --git a/modules/signal_processing/macros/iirmod.bin b/modules/signal_processing/macros/iirmod.bin
new file mode 100755
index 000000000..2c01d2fcf
--- /dev/null
+++ b/modules/signal_processing/macros/iirmod.bin
diff --git a/modules/signal_processing/macros/iirmod.sci b/modules/signal_processing/macros/iirmod.sci
new file mode 100755
index 000000000..74f3d81ab
--- /dev/null
+++ b/modules/signal_processing/macros/iirmod.sci
@@ -0,0 +1,49 @@
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) INRIA
+//
+// 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.1-en.txt
+
+function [la,grad,ind]=iirmod(x,ind)
+    //    p===> critere Lp
+    //    r=module des poles et zeros des filtres
+    //theta=argument des  "    "   "    "    "
+    //omega=frequences ou sont donnees les specifications des filtres
+    //wa=fonction de ponderation pour l'amplitude
+    //!
+    r=x(:,1);theta=x(:,2);
+    [m,n]=size(ad);if m>n,ad=ad';end;
+    [m,n]=size(omega);if m>n,omega=omega';end;
+    [m,n]=size(r);if n>m,r=r';m=n;end;
+    [m,n]=size(theta);if n>m,theta=theta';m=n;end;
+    m=m/2;
+    //
+    //VARIABLES LOCALES
+    unr=ones(r);unom=ones(omega);un=unr(1:m,:).*.unom;
+    n1=r.*.unom;n2=unr.*.omega;n3=theta.*.unom;n4=(unr+r.*r).*.unom;
+    cos1=[];cos2=[];sin1=[];sin2=[];
+    for k=(n2-n3),cos1=[cos1 cos(k)];sin1=[sin1 sin(k)];end;
+    for k=(n2+n3),cos2=[cos2 cos(k)];sin2=[sin2 sin(k)];end;
+    //
+    //FORMES QUADRATIQUES
+    k1=n4-2*n1.*cos1;k2=n4-2*n1.*cos2;
+    //
+    //AMPLITUDE
+    ampl=[];k3=k1(1:m,:);k4=k1(m+1:2*m,:);k5=k2(1:m,:);k6=k2(m+1:2*m,:);
+    a1=(k3.*k5)./(k4.*k6);
+    for k=a1,ampl=[ampl sqrt(prod(k))];end;
+    //
+    //GRADIENT DE L'AMPLITUDE
+    grrampl=(n1-cos1)./k1+(n1-cos2)./k2;
+    grrampl(m+1:2*m,:)=-grrampl(m+1:2*m,:);
+    grtampl=-(n1.*sin1)./k1+(n1.*sin2)./k2;
+    grtampl(m+1:2*m,:)=-grtampl(m+1:2*m,:);
+    //
+    //CRITERE D'ERREUR EN AMPLITUDE ET SON GRADIENT
+    a=ampl-ad;a1=a**(2*p);a1=a1.*wa;la=sum(a1);
+    a1=(a1./a)*2*p;a1=a1.*ampl;
+    grad=[grrampl*a1' grtampl*a1'];
+endfunction
diff --git a/modules/signal_processing/macros/intdec.bin b/modules/signal_processing/macros/intdec.bin
new file mode 100755
index 000000000..ee56ed02f
--- /dev/null
+++ b/modules/signal_processing/macros/intdec.bin
diff --git a/modules/signal_processing/macros/intdec.sci b/modules/signal_processing/macros/intdec.sci
new file mode 100755
index 000000000..dde6f2386
--- /dev/null
+++ b/modules/signal_processing/macros/intdec.sci
@@ -0,0 +1,112 @@
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) INRIA - 1990 - C. Bunks
+//
+// 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.1-en.txt
+
+function [y]=intdec(x,lom)
+    //y=intdec(x,lom)
+    //Changes the sampling rate of a 1D or 2D signal by the rates in lom
+    //  x      :Input sampled signal
+    //  lom    :For a 1D signal this is a scalar which gives the rate change
+    //         :For a 2D signal this is a 2-Vector of sampling rate
+    //         :changes lom=(col rate change,row rate change)
+    //  y      :Output sampled signal
+    //!
+
+    //Get dimensions of vectors
+
+    xsize=size(x);
+    xmin=min(x);xmax=max(x);
+    if xsize(1)==1 then, lom=[1 lom]; end,
+    if xsize(2)==1 then, lom=[lom 1]; end,
+
+    //compute sampling rate change as ratio of two integers
+
+    [l,m]=rat(lom);
+
+    //Assuming that the signal length is N (=xsize)
+    //the interpolated signal is N*l and the decimated
+    //signal is N*l/m.  The resulting output will have
+    //length int(N*l/m).
+
+    xlsize=xsize.*l;
+    xmsize=int(xlsize./m);
+
+    //Since the location of %pi in the frequency domain
+    //falls on a sample point for N even and between two
+    //sample points for N odd care must be taken to differentiate
+    //between the two cases in the following manipulations.
+
+    leven=2*(int(xsize/2)-xsize/2)+ones(xsize);
+    meven=2*(int(xmsize/2)-xmsize/2)+ones(xmsize);
+
+    //The position of %pi for the Fourier transform of the
+    //original signal is different for odd and even length signals.
+    //For an even length signal %pi is at the (N/2)+1 sample and
+    //for an odd length signal %pi is between the (N+1)/2 and the
+    //(N+1)/2 + 1 samples.
+
+    fp=int(xsize/2)+ones(xsize);
+    fpc=xsize-fp+leven;
+
+    fm=int(xmsize/2)+ones(xmsize);
+    fmc=fm-ones(fm)-meven;
+
+    //If the input is a constant then don't do the work
+
+    if xmax==xmin then,
+        y=xmax*ones(xmsize(1),xmsize(2));
+    else,
+
+        //For interpolation we, theoretically, would upsample x by putting l-1
+        //zeroes between each sample and then low pass filter at w=%pi.
+        //However, this corresponds to cutting the Fourier transform of the
+        //original signal into two parts at w=%pi and inserting l times the
+        //length of the signal in zeroes.
+        //
+        //For decimation we, theoretically, would low pass filter at the
+        //Nyquist frequency and then remove m-1 samples out of m of the
+        //time domain signal.  However, this corresponds to saving the
+        //N/m points of the signal at the two ends of the Fourier transform
+        //and then inverse transforming.
+
+        //Fourier transform the signal
+
+        xf=fft(x,-1);
+
+        //Re-assemble the correct portions of the frequency domain signal
+
+        if fm(1)<fp(1) then,//reduce row length (decimation)
+            xlf=[xf(1:fm(1),:);xf(xsize(1)-fmc(1)+1:xsize(1),:)];
+        else,
+            if xmsize(1)==xsize(1) then,//no changes in row length
+                xlf=xf;
+            else,//increase row length (interpolation)
+                xlf=[xf(1:fp(1),:);...
+                0*ones(xmsize(1)-fpc(1)-fp(1),xsize(2));...
+                xf(xsize(1)-fpc(1)+1:xsize(1),:)];
+            end,
+        end,
+        if fm(2)<fp(2) then,//decrease column length (decimation)
+            xlf=[xlf(:,1:fm(2)),xlf(:,xsize(2)-fmc(2)+1:xsize(2))];
+        else,
+            if xmsize(2)==xsize(2) then,//no changes in column length
+                xlf=xlf;
+            else,//increase column length (interpolation)
+                xlf=[xlf(:,1:fp(2)),...
+                0*ones(xmsize(1),xmsize(2)-fpc(2)-fp(2)),...
+                xlf(:,xsize(2)-fpc(2)+1:xsize(2))];
+            end,
+        end,
+
+        //recuperate new signal and rescale
+
+        y=real(fft(xlf,1));
+        y=prod(lom)*y;
+
+    end
+endfunction
diff --git a/modules/signal_processing/macros/kalm.bin b/modules/signal_processing/macros/kalm.bin
new file mode 100755
index 000000000..1f936e014
--- /dev/null
+++ b/modules/signal_processing/macros/kalm.bin
diff --git a/modules/signal_processing/macros/kalm.sci b/modules/signal_processing/macros/kalm.sci
new file mode 100755
index 000000000..baa54ec78
--- /dev/null
+++ b/modules/signal_processing/macros/kalm.sci
@@ -0,0 +1,35 @@
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) INRIA - 1988 - C. Bunks
+//
+// 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.1-en.txt
+
+function [x1,p1,x,p]=kalm(y,x0,p0,f,g,h,q,r)
+    //[x1,p1,x,p]=kalm(y,x0,p0,f,g,h,q,r)
+    //macro which gives the Kalman update and error variance
+    //Input to the macro is:
+    // f,g,h  :current system matrices
+    // q,r    :covariance matrices of dynamics and observation noise
+    // x0,p0  :state estimate and error variance at t=0 based
+    //        :on data up to t=-1
+    // y      :current observation
+    //
+    //Output from the macro is:
+    // x1,p1  :updated estimate and error covariance at t=1
+    //        :based on data up to t=0
+    // x,p    :updated estimate and error covariance at t=0
+    //        :based on data up to t=0
+    //!
+
+    if argn(2) <> 8 then
+        error(sprintf(gettext("%s: Wrong number of input argument(s): %d expected.\n"), "kalm", 8));
+    end
+    k=p0*h'*(h*p0*h'+r)**(-1);
+    p=(eye(p0)-k*h)*p0;
+    p1=f*p*f'+g*q*g';
+    x=x0+k*(y-h*x0);
+    x1=f*x;
+endfunction
diff --git a/modules/signal_processing/macros/lattn.bin b/modules/signal_processing/macros/lattn.bin
new file mode 100755
index 000000000..84266e8b5
--- /dev/null
+++ b/modules/signal_processing/macros/lattn.bin
diff --git a/modules/signal_processing/macros/lattn.sci b/modules/signal_processing/macros/lattn.sci
new file mode 100755
index 000000000..fb262e482
--- /dev/null
+++ b/modules/signal_processing/macros/lattn.sci
@@ -0,0 +1,116 @@
+// 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.1-en.txt
+
+function [la,lb]=lattn(n,p,mat_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(mat_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(mat_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=mat_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=[mat_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,
+        r1=flipdim(cv((p+1)*d+1:(p+2+j)*d,:), 1, d);
+        r2=flipdim(cv(p*d+1:(p+1+j)*d,:), 1, d);
+        r3=flipdim(cv((p-1-j)*d+1:p*d,:), 1, d);
+        r4=flipdim(cv((p-j)*d+1:(p+1)*d,:), 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/modules/signal_processing/macros/lattp.bin b/modules/signal_processing/macros/lattp.bin
new file mode 100755
index 000000000..e404ccb97
--- /dev/null
+++ b/modules/signal_processing/macros/lattp.bin
diff --git a/modules/signal_processing/macros/lattp.sci b/modules/signal_processing/macros/lattp.sci
new file mode 100755
index 000000000..790ba3865
--- /dev/null
+++ b/modules/signal_processing/macros/lattp.sci
@@ -0,0 +1,33 @@
+// 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.1-en.txt
+
+function [la,lb]=lattp(n,p,mat_cov)
+    // See lattn for more information
+    [l,d]=size(mat_cov);
+    id=eye(d);
+    [a,b]=lattn(n,0,mat_cov);
+    a=a(n);
+    b=b(n);
+    z=poly(0,"z");
+    la=list();
+    lb=list();
+    for j=0:p-1,
+        r1=flipdim(mat_cov((j+1)*d+1:(j+n+2)*d,:), 1, d);
+        r2=flipdim(mat_cov(j*d+1:(j+n+1)*d,:), 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
diff --git a/modules/signal_processing/macros/lev.bin b/modules/signal_processing/macros/lev.bin
new file mode 100755
index 000000000..5060c04b0
--- /dev/null
+++ b/modules/signal_processing/macros/lev.bin
diff --git a/modules/signal_processing/macros/lev.sci b/modules/signal_processing/macros/lev.sci
new file mode 100755
index 000000000..9a5d40416
--- /dev/null
+++ b/modules/signal_processing/macros/lev.sci
@@ -0,0 +1,48 @@
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) INRIA - 1988 - C. Bunks
+// Copyright (C) INRIA - 1991 - C. Bunks
+//
+// 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.1-en.txt
+
+function [ar,sigma2,rc]=lev(r)
+    //[ar,sigma2,rc]=lev(r)
+    //Resolve the Yule-Walker equations:
+    //
+    //       |r(0)   r(1)   ... r(N-1)|| a(1) | |sigma2|
+    //       |r(1)   r(0)   ... r(n-1)|| a(2) | |  0   |
+    //       |  :      :    ...    :  ||   :  |=|  0   |
+    //       |  :      :    ...    :  ||   :  | |  0   |
+    //       |r(N-1) r(N-2) ...  r(0) ||a(N-1)| |  0   |
+    //
+    //using Levinson's algorithm.
+    //  r      :Correlation coefficients
+    //  ar     :Auto-Regressive model parameters
+    //  sigma2 :Scale constant
+    //  rc     :Reflection coefficients
+    //!
+
+    //get the size of the correlation vector
+
+    rsize=max(size(r));
+    r=matrix(r,1,rsize);
+
+    //initialize levinson's algorithm
+
+    ar=-r(2)/r(1);
+    rc(1)=ar;
+    sigma2=(1-ar*conj(ar))*r(1);
+
+    //iterative solution to yule-walker equations
+
+    for k=2:rsize-1,
+        ak1(k)=-(r(k+1)+ar(1:k-1)'*r(k:-1:2)')/sigma2;
+        rc(k)=ak1(k);
+        ak1(1:k-1)=ar(1:k-1)+ak1(k)*conj(ar(k-1:-1:1));
+        sigma2=(1-ak1(k)*conj(ak1(k)))*sigma2;
+        ar=ak1;
+    end,
+endfunction
diff --git a/modules/signal_processing/macros/levin.bin b/modules/signal_processing/macros/levin.bin
new file mode 100755
index 000000000..e5f048391
--- /dev/null
+++ b/modules/signal_processing/macros/levin.bin
diff --git a/modules/signal_processing/macros/levin.sci b/modules/signal_processing/macros/levin.sci
new file mode 100755
index 000000000..f2f2eabb2
--- /dev/null
+++ b/modules/signal_processing/macros/levin.sci
@@ -0,0 +1,89 @@
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) INRIA - G. Le Vey
+//
+// 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.1-en.txt
+
+function [la, sig, lb] = levin(n, Cov)
+    // [la, sig, lb] = levin(n, Cov)
+    // macro which solves recursively on n
+    // the following Toeplitz system (normal equations)
+    //
+    //
+    //		|R1   R2   . . . Rn  |
+    //		|R0   R1   . . . Rn-1|
+    //		|R-1  R0   . . . Rn-2|
+    //		| .    .   . . .  .  |
+    // |I -A1 . -An|| .    .   . . .  .  |
+    //		| .    .   . . .  .  |
+    //		| .    .   . . .  .  |
+    //		|R2-n R3-n . . . R1  |
+    //		|R1-n R2-n . . . R0  |
+    //
+    //  where {Rk;k=1,nlag} is the sequence of nlag empirical covariances
+    //
+    //  n   : maximum order of the filter
+    //  Cov : matrix containing the Rk (d*d matrices for a
+    //      : d-dimensional process). It must be given the
+    //      : following way:
+    //
+    //		|  R0  |
+    //		|  R1  |
+    //		|  R2  |
+    //		|  .   |
+    //		|  .   |
+    //		| Rnlag|
+    //
+    //  la  : list-type variable, giving the successively calculated
+    //      : polynomials (degree 1 to degree n), with coefficients Ak
+    //  sig : list-type variable, giving the successive
+    //      : mean-square errors.
+    //!
+
+    [lhs, rhs] = argn(0);
+    if rhs <> 2 then
+        error(msprintf(gettext("%s: Wrong number of input argument(s): %d expected.\n"), "levin", 2));
+    end
+    [l, d] = size(Cov);
+    if d > l then
+        error(msprintf(gettext("%s: Wrong size for input argument #%d: A tall matrix expected.\n"), "levin", 2));
+    end
+    //
+    //   Initializations
+    //
+    a = eye(d, d);
+    b = a;
+    z = poly(0, "z");
+    la = list();
+    lb = list();
+    sig = list();
+    p = n+1;
+    cv = Cov;
+    for j=1:p
+        cv = [Cov(j*d+1:(j+1)*d, :)'; cv];
+    end
+    for j=0:n-1
+        //
+        //   Levinson algorithm
+        //
+        r1 = flipdim(cv((p+1)*d+1:(p+2+j)*d, :), 1, d);
+        r2 = flipdim(cv(p*d+1:(p+1+j)*d, :),     1, d);
+        r3 = flipdim(cv((p-1-j)*d+1:p*d, :),     1, d);
+        r4 = flipdim(cv((p-j)*d+1:(p+1)*d, :),   1, d);
+        c1 = coeff(a);
+        c2 = coeff(b);
+        sig1 = c1*r4;
+        gam1 = c2*r2;
+        k1 = (c1*r1)*inv(gam1);
+        k2 = (c2*r3)*inv(sig1);
+        a1 = a-k1*z*b;
+        b = -k2*a+z*b;
+        a = a1;
+        la(j+1) = a;
+        lb(j+1) = b;
+        sig(j+1) = sig1;
+    end
+endfunction
diff --git a/modules/signal_processing/macros/lib b/modules/signal_processing/macros/lib
new file mode 100755
index 000000000..cf370ed9c
--- /dev/null
+++ b/modules/signal_processing/macros/lib
diff --git a/modules/signal_processing/macros/lindquist.bin b/modules/signal_processing/macros/lindquist.bin
new file mode 100755
index 000000000..a38467568
--- /dev/null
+++ b/modules/signal_processing/macros/lindquist.bin
diff --git a/modules/signal_processing/macros/lindquist.sci b/modules/signal_processing/macros/lindquist.sci
new file mode 100755
index 000000000..01d506efb
--- /dev/null
+++ b/modules/signal_processing/macros/lindquist.sci
@@ -0,0 +1,45 @@
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) INRIA - 1989 - G. Le Vey
+//
+// 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.1-en.txt
+
+function [P,R,T]=lindquist(n,H,F,G,R0)
+    //[Pn,Rn,Tn]=lindquist(n,H,F,G,R0)
+    //macro which computes iteratively the minimal solution of the algebraic
+    //Riccati equation and gives the matrices Rn and Tt of the filter model,
+    //by the lindquist algorithm.
+    //   n     : number of iterations.
+    //   H,F,G : estimated triple from the covariance sequence of y.
+    //   R0    : E(yk*yk')
+    //   Pn    : solution of the Riccati equation after n iterations.
+    //   Rn,Tn : gain matrices of the filter.
+    //!
+    [d,m]=size(H);
+    //initialization
+    Gn=G;
+    Rn=R0;
+    Pn=zeros(m,m)
+    Kn=0*ones(m,d);
+
+    //recursion
+    for j=1:n,
+        //  Pn=Pn+Gn/Rn*Gn'
+        //  Kn=Pn*H'
+        Kn=Kn+Gn/Rn*Gn'*H';
+        r1=R0-H*Kn;
+        Rn=Rn-Gn'*H'/r1*H*Gn;
+        Gn=(F-(G-F*Kn)/r1*H)*Gn;
+    end
+
+    //gain matrices of the filter.
+    //P=Pn
+    //R=R0-H*P*H'
+    //T=(G-F*P*H')/R
+    R=R0-H*Kn
+    T=(G-F*Kn)/R
+    P=Kn
+endfunction
diff --git a/modules/signal_processing/macros/mese.bin b/modules/signal_processing/macros/mese.bin
new file mode 100755
index 000000000..b8615aaa0
--- /dev/null
+++ b/modules/signal_processing/macros/mese.bin
diff --git a/modules/signal_processing/macros/mese.sci b/modules/signal_processing/macros/mese.sci
new file mode 100755
index 000000000..4d7405abe
--- /dev/null
+++ b/modules/signal_processing/macros/mese.sci
@@ -0,0 +1,43 @@
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) INRIA - 1988 - C. Bunks
+// Copyright (C) INRIA - 1991 - C. Bunks
+//
+// 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.1-en.txt
+
+function [sm,fr]=mese(x,npts);
+    //<sm,fr]=mese(x [,npts]);
+    //Calculate the maximum entropy spectral estimate of x
+    //  x    :Input sampled data sequence
+    //  npts :Optional parameter giving number of points of fr and sm
+    //        (default is 256)
+    //  sm   :Samples of spectral estimate on the frequency grid fr
+    //  fr   :npts equally spaced frequency samples in [0,.5)
+    //!
+
+    //default evaluation
+
+    [rhs,lhs]=argn(0);
+    if lhs==1 then,
+        npts=256;
+    end,
+
+    //estimate autocorrelation function of x
+
+    Nx=length(x);
+    r=convol(x,x(Nx:-1:1))
+    r=r(Nx:-1:1)/Nx;
+
+    //get solution to the Yule-Walker equations
+
+    [ar,sigma2,rc]=lev(r);
+
+    //compute spectrum
+
+    ak=[1;ar];
+    [sf,fr]=frmag(ak,npts);
+    sm=sigma2*ones(sf)./(sf.*conj(sf));
+endfunction
diff --git a/modules/signal_processing/macros/mrfit.bin b/modules/signal_processing/macros/mrfit.bin
new file mode 100755
index 000000000..5d8b8a2b5
--- /dev/null
+++ b/modules/signal_processing/macros/mrfit.bin
diff --git a/modules/signal_processing/macros/mrfit.sci b/modules/signal_processing/macros/mrfit.sci
new file mode 100755
index 000000000..544a218e0
--- /dev/null
+++ b/modules/signal_processing/macros/mrfit.sci
@@ -0,0 +1,80 @@
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) INRIA
+//
+// 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.1-en.txt
+
+function [num,den]=mrfit(w,mod,r)
+    //Calling sequence:
+    //[num,den]=mrfit(w,mod,r)
+    //sys = mrfit(w,mod,r)
+    //
+    //w: vector of frequencies
+    //mod: vector of frequency response magnitude at these frequencies.
+    //weight: vector of weights given to each point
+    //
+    //Fits frequency response magnitude data points by a bi-stable transfer
+    //function
+    //G(s) = num(s)/den(s) of order r.
+    //
+    //  abs(freq(num,den,%i*w)) should be close to mod
+    //
+
+    function ww=mrfitdiff(ww)
+        //Utility fct
+        p=size(ww(:),"*");
+        ww=ww(2:p)-ww(1:p-1);
+    endfunction
+
+    w=w(:);mod=mod(:);
+    [LHS,RHS]=argn(0);
+    if w(1)==0 then w(1)=%eps;end
+
+    if r==0 then num=sum(mod)/size(mod,"*");den=1;return;end
+
+    sl0=round(log10(mod(2)/mod(1))/log10(w(2)/w(1)));w(1)=w(2)/10;
+    if sl0~=0 then mod(1)=mod(2)/10^sl0;end
+
+    n=length(w);
+    slinf=round(log10(mod($)/mod($-1))/log10(w($)/w($-1)));
+    w($)=w($-1)*10;
+    if slinf~=0 then mod($)=mod($-1)*10^slinf;end
+    logw=log10(w);logmod=log10(mod);delw=mrfitdiff(logw);delmod=mrfitdiff(logmod);
+
+    weight=ones(length(w),1);
+
+    junk=find(abs(mrfitdiff(delmod./delw)) > .6);
+    ind=1+junk;
+    if junk==[] then ind=[];end;
+    weight(ind)=10*ones(length(ind),1);
+
+    lwnew=[]; modnew=[];
+    for i=1:length(w)-1,
+        nadd=floor(delw(i)/.3)-1;
+        if nadd>0 then
+            slope=delmod(i)/delw(i);
+            lwnew=[lwnew logw(i)+.3*(1:nadd)];
+            modnew=[modnew logmod(i)+.3*slope*(1:nadd)];
+            weight=[weight ; ones(nadd,1)];
+        end
+    end
+    log10is=log(10);
+    w=[w ; exp(lwnew'*log10is)];
+    mod=[mod ; exp(modnew'*log10is)];
+
+    [w,ind] = gsort(-w);
+    w=-w; mod=mod(ind); weight=weight(ind);
+    ind=find(mrfitdiff(w)>0); ind=[ind(:);length(w)];
+    w=w(ind); mod=mod(ind); weight=weight(ind);
+
+    fresp=cepstrum(w,mod);
+
+    [num,den]=frfit(w,fresp,r,weight);
+    if LHS==1 then
+        num=syslin("c",num/den);
+    end
+
+endfunction
diff --git a/modules/signal_processing/macros/names b/modules/signal_processing/macros/names
new file mode 100755
index 000000000..baadbfdda
--- /dev/null
+++ b/modules/signal_processing/macros/names
@@ -0,0 +1,71 @@
+%k
+%sn
+analpf
+bilt
+buttmag
+casc
+cepstrum
+cheb1mag
+cheb2mag
+conv
+convol
+convol2d
+cspect
+czt
+detrend
+ell1mag
+eqfir
+eqiir
+faurre
+ffilt
+fft2
+fftshift
+filt_sinc
+filter
+find_freq
+findm
+frfit
+frmag
+fsfirlin
+group
+hank
+hilb
+hilbert
+idct
+idst
+ifft
+ifftshift
+iir
+iirgroup
+iirlp
+iirmod
+intdec
+kalm
+lattn
+lattp
+lev
+levin
+lindquist
+mese
+mrfit
+phc
+pspect
+remezb
+sincd
+srfaur
+srkf
+sskf
+system
+trans
+wfir
+wfir_gui
+wiener
+wigner
+window
+xcorr
+xcov
+yulewalk
+zpbutt
+zpch1
+zpch2
+zpell
diff --git a/modules/signal_processing/macros/phc.bin b/modules/signal_processing/macros/phc.bin
new file mode 100755
index 000000000..ef744324f
--- /dev/null
+++ b/modules/signal_processing/macros/phc.bin
diff --git a/modules/signal_processing/macros/phc.sci b/modules/signal_processing/macros/phc.sci
new file mode 100755
index 000000000..fca2c0f9d
--- /dev/null
+++ b/modules/signal_processing/macros/phc.sci
@@ -0,0 +1,32 @@
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) INRIA - 1989 - G. Le Vey
+//
+// 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.1-en.txt
+
+function [h,f,g]=phc(hk,d,r)
+    //[h,f,g]=phc(hk,d,r)
+    //macro which computes the matrices h, f, g by the principal hankel
+    //component approximation method, from the hankel matrix built
+    //from the covariance sequence of a stochastic process.
+    //   hk  : hankel matrix
+    //   d   : dimension of the observation
+    //   r   : desired dimension of the state vector
+    //       : for the approximated model
+    //   h,f,g  : relevant matrices of the state-space model
+    //!
+
+
+    [p,q]=size(hk);
+    [u,s,v]=svd(hk);
+    s=diag(sqrt(diag(s(1:r,1:r))));
+    obs=u(:,1:r)*s;
+    g=s*v(:,1:r)';
+    g=g(:,1:d);
+    ofl=obs(d+1:p,:);opr=obs(1:p-d,:);
+    f=opr\ofl;
+    h=obs(1:d,:);
+endfunction
diff --git a/modules/signal_processing/macros/pspect.bin b/modules/signal_processing/macros/pspect.bin
new file mode 100755
index 000000000..a1b543daf
--- /dev/null
+++ b/modules/signal_processing/macros/pspect.bin
diff --git a/modules/signal_processing/macros/pspect.sci b/modules/signal_processing/macros/pspect.sci
new file mode 100755
index 000000000..37d357ec5
--- /dev/null
+++ b/modules/signal_processing/macros/pspect.sci
@@ -0,0 +1,126 @@
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) INRIA - 1988 - C. Bunks
+//
+// 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.1-en.txt
+
+function [sm,cwp]=pspect(sec_step,sec_leng,wtype,x,y,wpar)
+    //[sm [,cwp] ]=pspect(sec_step,sec_leng,wtype,x [,y] [,wpar])
+    //Cross-spectral estimate between x and y if both are given
+    //and auto-spectral estimate of x otherwise.
+    //Spectral estimate obtained using the modified periodogram method
+    // x        :
+    //            - if vector, the first signal data
+    //            - amount of input data if scalar. In this case the first signal
+    //              is read by batch by a   user defined function xbatch=getx(n,offset).)
+    // y        : optional, if present cross correlation is made else auto correlation
+    //            - if vector, the second  signal data .It must have the same number of element than x
+    //            - if scalar, the second signal
+    //              is read by batch by a  user defined function ybatch=gety(n,offset).)
+    // sec_step : offset of each data window.
+    //            if sec_step==sec_leng/2 50% overlap is made
+    // sec_leng : Number of points of the window
+    // wtype    : window type ('re','tr','hm','hn','kr','ch')
+    // wpar     : optional window parameters
+    //          :         for wtype='kr', wpar>0
+    //          :         for wtype='ch', 0<wpar(1)<.5, wpar(2)>0
+    // sm       : power spectral estimate in the interval [0,1]
+    // cwp      : unspecified Chebyshev window parameter
+    //!
+
+    [lhs,rhs]=argn(0);
+
+
+    if sec_step>=sec_leng then
+        error(msprintf(gettext("%s: Incompatible input arguments #%d and #%d: Argument #%d expected to be less than argument #%d.\n"),"pspect",1,2,1,2));
+    end
+    if and(wtype<>["re","tr","hm","hn","kr","ch"]) then
+        error(msprintf(gettext("%s: Wrong value for input argument #%d: Must be in the set {%s}.\n"),..
+        "pspect",3,"''re'',''tr'',''hm'',''hn'',''kr'',''ch''"));
+    end
+
+    //Analyze calling sequence and construct window
+    if rhs==4 then,
+        //pspect(sec_step,sec_leng,wtype,x)
+        w=window(wtype,sec_leng);
+        crossFlag=%f; //autocorrelation
+    elseif rhs==5 then,
+        //pspect(sec_step,sec_leng,wtype,x,y) or pspect(sec_step,sec_leng,wtype,x,wpar)
+        if wtype=="kr" then,//pspect(sec_step,sec_leng,'kr',x,wpar)
+            wpar=y;
+            w=window(wtype,sec_leng,wpar);cwp=[];
+            crossFlag=%f; //autocorrelation
+        elseif wtype=="ch" then,//pspect(sec_step,sec_leng,'ch',x,wpar)
+            wpar=y;
+            [w,cwp]=window(wtype,sec_leng,wpar);
+            crossFlag=%f; //autocorrelation
+        else,//pspect(sec_step,sec_leng,wtype,x,y)
+            w=window(wtype,sec_leng);cwp=[];
+            crossFlag=%t; //cross-correlation
+        end,
+    else //pspect(sec_step,sec_leng,wtype,x,y,wpar)
+        [w,cwp]=window(wtype,sec_leng,wpar);
+        crossFlag=%t; //cross-correlation
+    end,
+    wpower=w*w';//the window energy
+
+    //Make x and y row vectors
+    x=matrix(x,1,-1);
+    if crossFlag then
+        y=matrix(y,1,-1);
+        if size(x,"*")<>size(y,"*") then
+            error(msprintf(gettext("%s: Arguments #%d and #%d must have the same sizes.\n"),"pspect",4,5));
+        end
+    end
+
+    //Average periodograms
+    sm=0*w;
+    if size(x,"*")==1 then  //Batch processing of x and y data
+
+        //get number of sections to be used
+        nsecs=int((x-sec_leng+sec_step)/sec_step);//x  contains the number of points of the signals
+        ovrlp=sec_leng-sec_step;
+        xd=[0*ones(1,sec_step) getx(ovrlp,1)];
+        if crossFlag then
+            yd=[0*ones(1,sec_step) gety(ovrlp,1)];
+        end
+        for k=1:nsecs
+            xd(1:ovrlp)=xd(sec_step+1:sec_leng);
+            xd(ovrlp+1:sec_leng)=getx(sec_step,sec_leng+(k-1)*sec_step+1);
+            xw=w.*(xd-(sum(xd)/sec_leng));
+            fx=fft(xw,-1);
+            if crossFlag then
+                yd(1:ovrlp)=yd(sec_step+1:sec_leng);
+                yd(ovrlp+1:sec_leng)=gety(sec_step,sec_leng+(k-1)*sec_step+1);
+                yw=w.*(yd-(sum(yd)/sec_leng));
+                sm=sm+fx.*conj(fft(yw,-1));
+            else
+                sm=sm+real(fx.*conj(fx));
+            end
+        end
+    else // Signal data given in x and y if cross-correlation
+
+        //get number of sections to be used
+        nsecs=int((size(x,"*")-sec_leng+sec_step)/sec_step);
+        ind=(1:sec_leng);
+        for k=1:nsecs
+            indd=ind+(k-1)*sec_step*ones(1:sec_leng);
+            xd=x(indd);
+            xe=w.*(xd-(sum(xd)/sec_leng));
+            fx=fft(xe,-1);
+            if crossFlag then
+                yd=y(indd);
+                ye=w.*(yd-(sum(yd)/sec_leng));
+                sm=sm+fx.*conj(fft(ye,-1));
+            else
+                sm=sm+real(fx.*conj(fx));
+            end
+        end
+    end
+
+    //Normalization
+    sm=sm/(nsecs*wpower);
+endfunction
diff --git a/modules/signal_processing/macros/remezb.bin b/modules/signal_processing/macros/remezb.bin
new file mode 100755
index 000000000..3722c4220
--- /dev/null
+++ b/modules/signal_processing/macros/remezb.bin
diff --git a/modules/signal_processing/macros/remezb.sci b/modules/signal_processing/macros/remezb.sci
new file mode 100755
index 000000000..c81072a30
--- /dev/null
+++ b/modules/signal_processing/macros/remezb.sci
@@ -0,0 +1,61 @@
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) INRIA - 1988 - C. Bunks
+//
+// 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.1-en.txt
+function an=remezb(nc,fg,ds,wt)
+    //an=remezb(nc,fg,ds,wt)
+    //Minimax approximation of a frequency domain
+    //magnitude response.  The approximation takes
+    //the form
+    //
+    //          h = sum[a(n)cos(wn)]
+    //
+    //for n=0,1,...,nc.  An FIR, linear-phase filter
+    //can be obtained from the output of the macro
+    //by using the following Scilab commands
+    //
+    //          hn(1:nc-1)=an(nc:-1:2)/2;
+    //          hn(nc)=an(1);
+    //          hn(nc+1:2*nc-1)=an(2:nc)/2;
+    //
+    //  nc :Number of cosine functions
+    //  fg :Grid of frequency points in [0,.5)
+    //  ds :Desired magnitude on grid fg
+    //  wt :Weighting function on error on grid fg
+    //  an :Cosine filter coefficients
+    //!
+
+
+    //get frequency grid size
+
+    ngrid=max(size(fg));
+
+    //check for errors in input
+
+    if min(fg)<0|max(fg)>0.5 then
+        error(msprintf(gettext("%s: Wrong value for input argument #%d: Must be in the interval [%s, %s].\n"),"remezb",2,"0","0.5"));
+    end,
+    dsize=max(size(ds));
+    wsize=max(size(wt));
+    if dsize<>ngrid then
+        error(msprintf(gettext("%s: Incompatible input arguments #%d and #%d: Same sizes expected.\n"),"remezb",2,3));
+    end,
+    if wsize<>ngrid then
+        error(msprintf(gettext("%s: Incompatible input arguments #%d and #%d: Same sizes expected.\n"),"remezb",2,4));
+    end,
+
+    //set up the initial guess for the extremal frequencies
+
+    iext=round(1:ngrid/nc:ngrid);
+    iext(nc+1)=ngrid;
+    iext(nc+2)=ngrid;
+
+    //call remez.f
+
+    an=remez(iext,ds,fg,wt);
+    an=an(1:nc)
+endfunction
diff --git a/modules/signal_processing/macros/sincd.bin b/modules/signal_processing/macros/sincd.bin
new file mode 100755
index 000000000..5c349d174
--- /dev/null
+++ b/modules/signal_processing/macros/sincd.bin
diff --git a/modules/signal_processing/macros/sincd.sci b/modules/signal_processing/macros/sincd.sci
new file mode 100755
index 000000000..59fb1e45e
--- /dev/null
+++ b/modules/signal_processing/macros/sincd.sci
@@ -0,0 +1,42 @@
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) INRIA - 1989 - G. Le Vey
+//
+// 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.1-en.txt
+
+function [s]=sincd(n,flag)
+    //<s>=sincd(n,flag)
+    //macro which calculates the function Sin(N*x)/Sin(x)
+    //  n    : value of N in the preceding function
+    //  flag : flag=1 : the function is centered around the origin
+    //       : flag=2 : the function is delayed by %pi/(2*n)
+    //  s    : vector of values of the function on a dense
+    //       : grid of frequencies
+    //!
+
+    npt=4*n;
+    pas=%pi/npt;
+    om=0:pas:%pi;
+    eps=(-1)**(n-1);
+    select flag
+    case 1,
+        s1=sin(n*om);s2=sin(om);
+        s1(1)=n;s2(1)=1;s1(npt+1)=n*eps;s2(npt+1)=1;
+        s=s1./s2;
+        s=[s(npt+1:-1:2) s];
+        s=s/n;
+    case 2,
+        om=om-ones(om)*%pi/2/n;
+        s1=sin(n*om);
+        s2=sin(om);
+        s1(3)=n;s2(3)=1;
+        s=s1./s2;
+        s=[eps*s s(2:npt+1)];
+        s=s/n;
+    else
+        error(msprintf(gettext("%s: Wrong value for input argument #%d: Must be in the set {%s}.\n"),"sincd",2,"1,2"));
+    end;
+endfunction
diff --git a/modules/signal_processing/macros/srfaur.bin b/modules/signal_processing/macros/srfaur.bin
new file mode 100755
index 000000000..69c39a6e2
--- /dev/null
+++ b/modules/signal_processing/macros/srfaur.bin
diff --git a/modules/signal_processing/macros/srfaur.sci b/modules/signal_processing/macros/srfaur.sci
new file mode 100755
index 000000000..8855f1283
--- /dev/null
+++ b/modules/signal_processing/macros/srfaur.sci
@@ -0,0 +1,51 @@
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) INRIA
+//
+// 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.1-en.txt
+function [p,s,t,l,rt,tt]=srfaur(h,f,g,r0,n,p,s,t,l)
+    //square-root algorithm for the algebraic Riccati equation.
+    //
+    //   h,f,g  : convenient matrices of the state-space model.
+    //   r0     : E(yk*yk').
+    //   n      : number of iterations.
+    //
+    //   p      : estimate of the solution after n iterations.
+    //   s,t,l  : intermediate matrices for
+    //          : successive iterations;
+    //   rt,tt  : gain matrices of the filter model after n iterations.
+    //
+    //          : p,s,t,l may be given as input if more than one recursion
+    //          : is desired (evaluation of intermediate values of p, e.g.).
+    //!
+
+    [lhs,rhs]=argn(0);
+    [d,m]=size(h);
+    if rhs==5,
+        //initializations
+        r0=.5*(r0+r0');
+        s=sqrtm(r0);
+        t=(g/s)';s=s';
+        l=-%i*t;
+        p=l'*l;
+    else,
+        if rhs<>9 then
+            error(msprintf(gettext("%s: Wrong number of input argument(s): %d expected.\n"),"srfaur",9));
+        end;
+    end;
+    //recursion
+    for j=1:n,
+        a=[s t;l*h' l*f'];
+        [q,r]=qr(a);
+        s=r(1:d,1:d);
+        t=r(1:d,d+1:d+m);
+        l=r(d+1:2*d,d+1:d+m);
+        p=p+l'*l;
+    end;
+    //gain matrices of the filter.
+    rt=r0-h*p*h';
+    tt=(g-f*p*h')*inv(rt);
+endfunction
diff --git a/modules/signal_processing/macros/srkf.bin b/modules/signal_processing/macros/srkf.bin
new file mode 100755
index 000000000..581f8113b
--- /dev/null
+++ b/modules/signal_processing/macros/srkf.bin
diff --git a/modules/signal_processing/macros/srkf.sci b/modules/signal_processing/macros/srkf.sci
new file mode 100755
index 000000000..7f6792c78
--- /dev/null
+++ b/modules/signal_processing/macros/srkf.sci
@@ -0,0 +1,39 @@
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) INRIA - 1988 - C. Bunks
+//
+// 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.1-en.txt
+
+function [x1,p1]=srkf(y,x0,p0,f,h,q,r)
+    //square root kalman filter algorithm
+    //Input to the macro is:
+    // f,h    :current system matrices
+    // q,r    :covariance matrices of dynamics
+    //        :and observation noise
+    // x0,p0  :state estimate and error variance
+    //        :at t=0 based on data up to t=-1
+    // y      :current observation
+    //
+    //Output from the macro is:
+    // x1,p1  :updated estimate and error covariance
+    //        :at t=1 based on data up to t=0
+    //!
+
+    n=max(size(x0));
+    p=max(size(y));
+
+    j=[chol(r)',0*r];
+    g=[0*q,chol(q)'];
+
+    mat=[h*p0,j;f*p0,g];
+    [q,tmat]=qr(mat')';
+    p1=tmat(p+1:p+n,p+1:p+n);
+    k=tmat(p+1:p+n,1:p);
+    re=tmat(1:p,1:p);
+
+    epsilon=y-h*x0;
+    x1=f*x0+k*(re**(-1))*epsilon;
+endfunction
diff --git a/modules/signal_processing/macros/sskf.bin b/modules/signal_processing/macros/sskf.bin
new file mode 100755
index 000000000..2b461ffa8
--- /dev/null
+++ b/modules/signal_processing/macros/sskf.bin
diff --git a/modules/signal_processing/macros/sskf.sci b/modules/signal_processing/macros/sskf.sci
new file mode 100755
index 000000000..2937e093c
--- /dev/null
+++ b/modules/signal_processing/macros/sskf.sci
@@ -0,0 +1,30 @@
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) INRIA
+//
+// 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.1-en.txt
+
+function [xe, x]=sskf(y,f,h,q,r,x0)
+    //<xe>=sskf(y,f,h,q,r,x0)
+    //steady-state kalman filter
+    // y   :data in form [y0,y1,...,yn], yk a column vector
+    // f   :system matrix dim(NxN)
+    // h   :observations matrix dim(NxM)
+    // q   :dynamics noise matrix dim(NxN)
+    // r   :observations noise matrix dim(MxM)
+    //
+    // xe  :estimated state
+    //!
+
+    //get steady-state Kalman gain
+
+    x=ricc(f',h'/r*h,q,"disc") // steady state err cov
+    k=x*h'/(h*x*h'+r)
+
+    // estimate state
+    kfd=(eye(f)-k*h)*f;
+    [xe]=ltitr(kfd,k,y,x0);
+endfunction
diff --git a/modules/signal_processing/macros/system.bin b/modules/signal_processing/macros/system.bin
new file mode 100755
index 000000000..526e96825
--- /dev/null
+++ b/modules/signal_processing/macros/system.bin
diff --git a/modules/signal_processing/macros/system.sci b/modules/signal_processing/macros/system.sci
new file mode 100755
index 000000000..ae80bc280
--- /dev/null
+++ b/modules/signal_processing/macros/system.sci
@@ -0,0 +1,43 @@
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 1988 - INRIA - C. Bunks
+// Copyright (C) DIGITEO - 2011 - Allan CORNET
+//
+// 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.1-en.txt
+//
+function [x1, y] = system(x0, f, g, h, q, r)
+    //<x1,y>=system(x0,f,g,h,q,r)
+    //define system macro which generates the next
+    //observation given the old state
+    //  x0 :Input state vector
+    //  f  :System matrix
+    //  g  :Input matrix
+    //  h  :Output matrix
+    //  q  :Input noise covariance matrix
+    //  r  :Output noise covariance matrix
+    //  x1 :Output state vector
+    //  y  :Output observation
+    //System recursively calculates
+    //
+    //     x1=f*x0+g*u
+    //      y=h*x0+v
+    //
+    //where u is distributed N(0,q)
+    //and v is distribute N(0,r).
+
+    [lhs, rhs] = argn(0);
+    if rhs == 0 then
+        error(999, msprintf(_("%s: Wrong number of input argument(s).\n"), "system"));
+    end
+
+    rand("normal");
+    q2 = chol(q);
+    r2 = chol(r);
+    u = q2' * rand(ones(x0));
+    v = r2' * rand(ones(x0));
+    x1 = f * x0 + g * u;
+    y = h * x0 + v;
+endfunction
diff --git a/modules/signal_processing/macros/trans.bin b/modules/signal_processing/macros/trans.bin
new file mode 100755
index 000000000..22de8f5f3
--- /dev/null
+++ b/modules/signal_processing/macros/trans.bin
diff --git a/modules/signal_processing/macros/trans.sci b/modules/signal_processing/macros/trans.sci
new file mode 100755
index 000000000..148915454
--- /dev/null
+++ b/modules/signal_processing/macros/trans.sci
@@ -0,0 +1,105 @@
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) INRIA - 1988 - C. Bunks
+//
+// 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.1-en.txt
+
+function [hzt,zt,gt]=trans(pd,zd,gd,tr_type,frq)
+    //hzt=trans(pd,zd,gd,tr_type,frq)
+    //macro for transforming standardized low-pass filter into
+    //one of the following filters:
+    //     low-pass, high-pass, band-pass, stop-band.
+    // hz      :input polynomial
+    // tr_type :type of transformation
+    // frq     :frequency values
+    // hzt     :output polynomial
+    //!
+    if and(argn(1)<>[1 3]) then
+        error(msprintf(gettext("%s: Wrong number of output arguments: %d or %d expected.\n"),"trans",1,3))
+    end
+    select argn(2)
+    case 3 then //trans(hz,tr_type,frq): filter given by a siso tranfer function
+        hz=pd
+        if typeof(hz)<>"rational" then
+            error(msprintf(_("%s: Wrong type for input argument #%d: rational fraction array expected.\n"),"trans",1))
+        end
+        if size(hz,"*")<>1 then
+            error(msprintf(_("%s: Wrong size for input argument #%d: A single input, single output system expected.\n"),"trans",1))
+        end
+        tr_type=zd;tr_pos=2;
+        frq=gd;frq_pos=3;
+        pd=roots(hz.den)
+        zd=roots(hz.num)
+        gd=coeff(hz.num,degree(hz.num))/coeff(hz.den,degree(hz.den))
+    case 5 then //trans(pd,zd,gd,tr_type,frq):  filter given by zeros,poles and gain
+        if type(pd)<>1 then
+            error(msprintf(_("%s: Wrong type for input argument #%d: Array of floating point numbers expected.\n"),"trans",1))
+        end
+        if type(zd)<>1 then
+            error(msprintf(_("%s: Wrong type for input argument #%d: Array of floating point numbers expected.\n"),"trans",2))
+        end
+        if type(gd)<>1|size(gd,"*")<>1 then
+            error(msprintf(_("%s: Wrong size for input argument #%d: A scalar expected.\n"),"trans",3))
+        end
+        tr_pos=4;
+        frq_pos=5;
+    else
+        error(msprintf(_("%s: Wrong number of input arguments: %d or %d expected.\n"),"trans",3,5))
+    end
+    if and(tr_type<>["lp","hp","bp","sb"]) then
+        error(msprintf(_("%s: Wrong value for input argument #%d: Must be in the set {%s}.\n"),"trans",tr_pos,"''lp'',''hp'',''bp'',''sb''"))
+    end
+    if type(frq)<>1 then
+        error(msprintf(_("%s: Wrong type for input argument #%d: Array of floating point numbers expected.\n"),"trans",frq_pos))
+    end
+    if or(tr_type==["lp","hp"]) then
+        if and(size(frq,"*")<>[1,2]) then
+            error(msprintf(_("%s: Wrong size for input argument #%d: A scalar expected.\n"),"trans",frq_pos))
+        end
+        fu=frq(1)
+    else
+        if size(frq,"*")<>2 then
+            error(msprintf(_("%s: Wrong size for input argument #%d: A %d elements array expected.\n"),"trans",frq_pos,2))
+        end
+        fu=frq(1);fl=frq(2);
+    end
+
+    z=poly(0,"z");fc=.25;
+    //make filter type using all-pass change of variables
+    select tr_type
+    case "lp" then
+        alpha=sin(%pi*(fc-fu))/sin(%pi*(fc+fu));
+        num=z-alpha;
+        den=1-alpha*z;
+    case "hp" then
+        alpha=-cos(%pi*(fc-fu))/cos(%pi*(fc+fu));
+        num=-(1+alpha*z);
+        den=z+alpha;
+    case "bp" then
+        k=tan(%pi*fc)/tan(%pi*(fu-fl));
+        alpha=cos(%pi*(fu+fl))/cos(%pi*(fu-fl));
+        num=-((k+1)-2*alpha*k*z+(k-1)*z^2);
+        den=(k+1)*z^2-2*alpha*k*z+(k-1);
+    case "sb" then
+        k=tan(%pi*fc)*tan(%pi*(fu-fl));
+        alpha=cos(%pi*(fu+fl))/cos(%pi*(fu-fl));
+        num=(k+1)-2*alpha*z+(1-k)*z^2;
+        den=(k+1)*z^2-2*alpha*z+(1-k);
+    end
+    if tr_type == "lp" & fu == 0.25 then //no change of variables is necessary if user wants a low pass filter with cut off frequency of 0.25Hz.
+        pt=pd;
+        zt=zd;
+        gt=gd;
+    else
+        [pt,zt,gt]=bilt(pd,zd,gd,num,den);
+    end
+    if argn(1)==1 then
+        hzt=rlist(gt*real(poly(zt,"z")),real(poly(pt,"z")),"d");
+    else
+        hzt=pt(:)
+        zt=zt(:)
+    end
+endfunction
diff --git a/modules/signal_processing/macros/wfir.bin b/modules/signal_processing/macros/wfir.bin
new file mode 100755
index 000000000..20ddc0fed
--- /dev/null
+++ b/modules/signal_processing/macros/wfir.bin
diff --git a/modules/signal_processing/macros/wfir.sci b/modules/signal_processing/macros/wfir.sci
new file mode 100755
index 000000000..180d4fb9c
--- /dev/null
+++ b/modules/signal_processing/macros/wfir.sci
@@ -0,0 +1,105 @@
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) INRIA - 1988 - C. Bunks
+//
+// 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.1-en.txt
+
+function [wft,wfm,fr]=wfir(ftype,forder,cfreq,wtype,fpar)
+    //[wft,wfm,fr]=wfir(ftype,forder,cfreq,wtype,fpar)
+    //Macro which makes linear-phase, FIR low-pass, band-pass,
+    //high-pass, and stop-band filters
+    //using the windowing technique.
+    //Works interactively if called with no arguments.
+    //  ftype  :Filter type ('lp','hp','bp','sb')
+    //  forder :Filter order (pos integer)(odd for ftype='hp' or 'sb')
+    //  cfreq  :2-vector of cutoff frequencies (0<cfreq(1),cfreq(2)<.5)
+    //         :only cfreq(1) is used when ftype='lp' or 'hp'
+    //  wtype  :Window type ('re','tr','hm','hn','kr','ch')
+    //  fpar   :2-vector of window parameters
+    //         :     Kaiser window: fpar(1)>0 fpar(2)=0
+    //         :     Chebyshev window: fpar(1)>0 fpar(2)<0 or
+    //         :                       fpar(1)<0 0<fpar(2)<.5
+    //  wft    :Time domain filter coefficients
+    //  wfm    :Frequency domain filter response on the grid fr
+    //  fr     :Frequency grid
+    //!
+    wft=[];wfm=[];fr=[]
+
+    FT=["lp","hp","bp","sb"]
+
+    [lhs,rhs]=argn(0);
+
+    if rhs<=0 then,
+        //if macro called with no arguments query user for values
+        [ok,values,exprs]=wfir_gui()
+        if ~ok then return,end
+        ftype=values.ftype
+        forder=values.forder
+        cfreq=[values.low,values.high]/values.freq_ech;
+        wtype=values.wtype
+        fpar=values.fpar
+
+        fl=values.low/values.freq_ech;
+        fh=values.high/values.freq_ech;
+    else
+        //check arguments of macro call
+        if and(ftype<>FT) then
+            error(msprintf(_("%s: Wrong value for input argument #%d: Must be in the set {%s}.\n"),"wfir",1,strcat(FT,",")))
+        end
+        if type(forder)<>1|size(forder,"*")<>1|~isreal(forder)|size(forder,"*")<>1|int(forder)<>forder|forder<1 then
+            error(msprintf(_("%s: Wrong type for input argument #%d: A positive integer expected.\n"),"wfir",2))
+        end
+        if or(ftype==["hp" "sb"]) then
+            if 2*int(forder/2)==forder then
+                error(msprintf(_("%s:  Wrong value for input argument #%d: Must be odd.\n"),"wfir",2))
+            end
+        end
+
+        if type(cfreq)<>1|~isreal(cfreq) then
+            error(msprintf(_("%s: Wrong type for input argument #%d: A %d-by-%d real vector expected.\n"),"wfir",3,1,2))
+        end
+        if or(ftype==["hp" "lp"]) then
+            if size(cfreq,"*")==0| size(cfreq,"*")>2 then
+                error(msprintf(_("%s: Wrong size for input argument #%d: A %d-by-%d real vector expected.\n"),"wfir",3,1,2))
+            end
+            if cfreq(1)<0|cfreq(1)>0.5 then
+                error(msprintf(_("%s: Wrong value for input argument #%d: Must be in the interval %s.\n"),"wfir",3,"]0,0.5["))
+            end
+            fl=cfreq(1);fh=[]
+        else
+            if size(cfreq,"*")<>2 then
+                error(msprintf(_("%s: Wrong size for input argument #%d: A %d-by-%d real vector expected.\n"),"wfir",3,1,2))
+            end
+            if or(cfreq<0|cfreq>0.5) then
+                error(msprintf(_("%s: Wrong value for input argument #%d: Must be in the interval %s.\n"),"wfir",3,"]1,0.5["))
+            end
+            if cfreq(1)>=cfreq(2) then
+                error(msprintf(_("%s: Elements of %dth argument must be in increasing order.\n"),"wfir",3))
+            end
+            fl=cfreq(1);
+            fh=cfreq(2);
+        end
+    end
+
+
+    //Calculate window coefficients
+
+    [win_l,cwp]=window(wtype,forder,fpar);
+    [dummy,forder]=size(win_l);
+
+    //Get forder samples of the appropriate filter type
+
+    hfilt=ffilt(ftype,forder,fl,fh);
+
+    //Multiply window with sinc function
+
+    wft=win_l.*hfilt;
+
+    //Calculate frequency response of the windowed filter
+
+    [wfm,fr]=frmag(wft,256);
+endfunction
+
diff --git a/modules/signal_processing/macros/wfir_gui.bin b/modules/signal_processing/macros/wfir_gui.bin
new file mode 100755
index 000000000..5b6f7fa50
--- /dev/null
+++ b/modules/signal_processing/macros/wfir_gui.bin
diff --git a/modules/signal_processing/macros/wfir_gui.sci b/modules/signal_processing/macros/wfir_gui.sci
new file mode 100755
index 000000000..25e8624b4
--- /dev/null
+++ b/modules/signal_processing/macros/wfir_gui.sci
@@ -0,0 +1,843 @@
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2011-2014 - INRIA - Serge Steer
+//
+// 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.1-en.txt
+
+function [ok,values,exprs]=wfir_gui(exprs)
+    FT=["lp","hp","bp","sb"]
+    WT=["re","tr","hn","hm","kr","ch","ch"]
+    //  errcatch(-1,'continue')
+    if argn(2)<1 then
+        exprs=["""lp""";
+        """re"""
+        "48";
+        "1";
+        "0.05"
+        "0.45";
+        "0.5"]
+    end
+    ftype=exprs(1)
+    wtype=exprs(2)
+    if execstr("forder="+exprs(3)+";freq_ech="+exprs(4)+";low="+exprs(5)+";high="+exprs(6)+";fp="+exprs(7),"errcatch")<>0 then
+        values=[]
+        exprs=[]
+        ok=%f
+        return
+    end
+    global ret
+    ret=0
+
+    margin_x     = 5;      // Horizontal margin between each elements
+    margin_y     = 5;      // Vertical margin between each elements
+    button_w     = 70;
+    button_h     = 30;
+    frame_h     = 300;
+    frame_w     = 440;
+    label_h      = 20;
+    list_h       = 120;
+    list_w       = 100;
+    list_incr    = list_w;
+    axes_w       = 2*margin_x+frame_w;
+    axes_h       = 2*margin_y+frame_h;// Frame height
+    defaultfont  = "arial"; // Default Font
+
+    fig_id=max(winsid())+1
+    fig = scf(fig_id)
+
+    // Remove Scilab graphics menus & toolbar
+    //  drawlater (bug)
+    delmenu(fig.figure_id, gettext("&File"));
+    delmenu(fig.figure_id, gettext("&Tools"));
+    delmenu(fig.figure_id, gettext("&Edit"));
+    delmenu(fig.figure_id, gettext("&?"));
+    toolbar(fig.figure_id, "off");
+    fig.axes_size       = [axes_w axes_h];
+
+
+    fig.background      = addcolor([0.8 0.8 0.8]);
+    fig.figure_name     = _("WFIR settings");
+    ax=fig.children;
+    ax.background= fig.background ;
+
+    gui=uicontrol( ...
+    "parent"              , fig,...
+    "style"               , "frame",...
+    "units"               , "pixels",...
+    "position"            , [0 0 axes_w axes_h],...
+    "background"          , [1 1 1]*0.8, ...
+    "visible"             , "on");
+
+    yo=axes_h-button_h;
+
+
+    yft=frame_h-margin_y-label_h
+    xft=margin_x;
+    //Filter type selection -----------------------------------------------------------------
+    listf_h       = 80;
+    listf_w       = 150;
+    label_w       = 75;
+    uicontrol( ...
+    "parent"              , gui,...
+    "style"               , "text",...
+    "string"              , _("Filter type"),...
+    "units"               , "pixels",...
+    "position"            , [xft yft label_w label_h],....
+    "fontname"            , defaultfont,...
+    "fontunits"           , "points",...
+    "fontsize"            , 12,...
+    "fontweight"          , "bold", ...
+    "horizontalalignment" , "left", ...
+    "background"          , [1 1 1]*0.8, ...
+    "visible"             , "on");
+    Ftv=uicontrol( ...
+    "parent"              , gui,...
+    "style"               , "edit",...
+    "string"              , ftype,...
+    "units"               , "pixels",...
+    "position"            , [xft+label_w+margin_x yft listf_w-label_w-margin_x label_h],....
+    "fontname"            , defaultfont,...
+    "fontunits"           , "points",...
+    "fontsize"            , 12,...
+    "fontweight"          , "bold", ...
+    "horizontalalignment" , "left", ...
+    "background"          , [1 1 1]*0.8, ...
+    "callback"            , "wfirSetFilterType()", ...
+    "visible"             , "on");
+    Filtertype=uicontrol( ...
+    "parent"              , gui,...
+    "style"               , "listbox",...
+    "string"              , [_("Low pass");_("High pass");_("Band pass");_("Stop Band")],...
+    "value"               , find(FT==evstr(ftype)),...
+    "units"               , "pixels",...
+    "position"            , [xft,yft-listf_h,listf_w,listf_h],...
+    "fontname"            , defaultfont,...
+    "fontunits"           , "points",...
+    "fontsize"            , 12,...
+    "horizontalalignment" , "left", ...
+    "BackgroundColor"     , [1 1 1], ...
+    "ForegroundColor"     , [1 1 1]*0, ...
+    "visible"             , "on", ...
+    "callback"            , "wfirSetFilterType()");
+    set(Ftv,"userdata",list([],Filtertype))
+
+    //Window type selection -----------------------------------------------------------------
+    ywt=yft-listf_h-4*margin_y
+    xwt=xft
+    listw_h       = 140;
+    listw_w       = listf_w;
+    uicontrol( ...
+    "parent"              , gui,...
+    "style"               , "text",...
+    "string"              , _("Window type"),...
+    "units"               , "pixels",...
+    "position"            , [xwt ywt label_w label_h],....
+    "fontname"            , defaultfont,...
+    "fontunits"           , "points",...
+    "fontsize"            , 12,...
+    "fontweight"          , "bold", ...
+    "horizontalalignment" , "left", ...
+    "background"          , [1 1 1]*0.8, ...
+    "visible"             , "on");
+    Wtv=uicontrol( ...
+    "parent"              , gui,...
+    "style"               , "edit",...
+    "string"              , wtype,...
+    "units"               , "pixels",...
+    "position"            , [xwt+label_w+margin_x ywt listw_w-label_w-margin_x label_h],....
+    "fontname"            , defaultfont,...
+    "fontunits"           , "points",...
+    "fontsize"            , 12,...
+    "fontweight"          , "bold", ...
+    "horizontalalignment" , "left", ...
+    "background"          , [1 1 1]*0.8, ...
+    "callback"            , "wfirSetWindowType()",...
+    "visible"             , "on");
+
+
+
+    W=[_("Rectangular");_("Triangular");_("Hanning");
+    _("Hamming");_("Kaiser");_("Chebychev main lobe");_("Chebychev side lobe")];
+    index=find(WT==evstr(wtype))
+    if index==[] then index=1,end//for debug
+    Windowtype=uicontrol( ...
+    "parent"              , gui,...
+    "style"               , "listbox",...
+    "string"              , W,...
+    "value"               , find(WT==evstr(wtype)),...
+    "units"               , "pixels",...
+    "position"            , [xwt,ywt-listw_h,listw_w,listw_h],...
+    "fontname"            , defaultfont,...
+    "fontunits"           , "points",...
+    "fontsize"            , 12,...
+    "horizontalalignment" , "left", ...
+    "BackgroundColor"     , [1 1 1], ...
+    "ForegroundColor"     , [1 1 1]*0, ...
+    "visible"             , "on", ...
+    "callback"            , "wfirSetWindowType()");
+    set(Wtv,"userdata",list([],Windowtype))
+
+    //Sampling Frequency -----------------------------------------------------------------
+    ysf=frame_h-margin_y-label_h
+    xsf=xft+listf_w+2*margin_x
+    labelsf_w=180
+    editsf_w=80
+    sfreq=1
+    uicontrol( ...
+    "parent"              , gui,...
+    "style"               , "text",...
+    "string"              ,  _("Sampling Frequency (Hz)"),...
+    "units"               , "pixels",...
+    "position"            , [xsf ysf labelsf_w label_h],....
+    "fontname"            , defaultfont,...
+    "fontunits"           , "points",...
+    "fontsize"            , 12,...
+    "fontweight"          , "bold", ...
+    "horizontalalignment" , "left", ...
+    "background"          , [1 1 1]*0.8, ...
+    "visible"             , "on");
+    Sfreq =uicontrol( ...
+    "parent"              , gui,...
+    "style"               , "edit",...
+    "string"               , exprs(4),...
+    "units"               , "pixels",...
+    "position"            , [xsf+margin_x+labelsf_w,ysf,editsf_w,label_h],...
+    "fontname"            , defaultfont,...
+    "fontunits"           , "points",...
+    "fontsize"            , 12,...
+    "horizontalalignment" , "left", ...
+    "BackgroundColor"     , [1 1 1], ...
+    "ForegroundColor"     , [1 1 1]*0, ...
+    "callback"            , "wfirSetSamplingFrequency()", ...
+    "userdata"            ,sfreq,...
+    "visible"             , "on")
+    //Filter order -----------------------------------------------------------------
+    yfo=ysf-margin_y-label_h
+    xfo=xsf
+    labelfo_w=labelsf_w;
+    editfo_w= editsf_w;
+    scroll_w=160
+    uicontrol( ...
+    "parent"              , gui,...
+    "style"               , "text",...
+    "string"              ,  _("Filter Order"),...
+    "units"               , "pixels",...
+    "position"            , [xfo yfo labelfo_w label_h],....
+    "fontname"            , defaultfont,...
+    "fontunits"           , "points",...
+    "fontsize"            , 12,...
+    "fontweight"          , "bold", ...
+    "horizontalalignment" , "left", ...
+    "background"          , [1 1 1]*0.8, ...
+    "visible"             , "on");
+    Forderv=uicontrol( ...
+    "parent"              , gui,...
+    "style"               , "edit",...
+    "string"              , exprs(3),...
+    "units"               , "pixels",...
+    "position"            , [xfo,yfo-label_h-margin_y,editfo_w,label_h],...
+    "fontname"            , defaultfont,...
+    "fontunits"           , "points",...
+    "fontsize"            , 12,...
+    "horizontalalignment" , "left", ...
+    "BackgroundColor"     , [1 1 1], ...
+    "ForegroundColor"     , [1 1 1]*0, ...
+    "callback"            , "wfirValue2Sliderpos()", ...
+    "visible"             , "on")
+
+    Forders=uicontrol( ...
+    "parent"              , gui,...
+    "style"               , "slider",...
+    "min"                 , 1,...
+    "max"                 , 1000,...
+    "sliderstep"          , [1 100],...
+    "value"               , forder,...
+    "units"               , "pixels",...
+    "position"            , [xfo+editfo_w+margin_x,yfo-label_h-margin_y,scroll_w,label_h],...
+    "fontname"            , defaultfont,...
+    "fontunits"           , "points",...
+    "fontsize"            , 12,...
+    "horizontalalignment" , "left", ...
+    "BackgroundColor"     , [1 1 1], ...
+    "ForegroundColor"     , [1 1 1]*0, ...
+    "userdata"            , list(Forderv,1), ...
+    "callback"            , "wfirSliderpos2ValueI()",...
+    "visible"             , "on")
+    set(Forderv,"userdata",Forders);
+
+    //Low cutoff frequency -----------------------------------------------------------------
+    ylcf=yfo-2*(margin_y+label_h)
+    xlcf=xsf
+    labellcf_w=labelsf_w;
+    editlcf_w= editsf_w;
+    if ftype=="""hp""" then visible="off";else visible="on";end
+    Lcfl=uicontrol( ...
+    "parent"              , gui,...
+    "style"               , "text",...
+    "visible"             , visible,...
+    "string"              ,  _("Low cutoff frequency (Hz)"),...
+    "units"               , "pixels",...
+    "position"            , [xlcf ylcf labellcf_w label_h],....
+    "fontname"            , defaultfont,...
+    "fontunits"           , "points",...
+    "fontsize"            , 12,...
+    "fontweight"          , "bold", ...
+    "horizontalalignment" , "left", ...
+    "background"          , [1 1 1]*0.8);
+    fact=1000;
+    Lcfv=uicontrol( ...
+    "parent"              , gui,...
+    "style"               , "edit",...
+    "visible"             , visible,...
+    "units"               , "pixels",...
+    "string"              , exprs(5),...
+    "position"            , [xlcf,ylcf-label_h-margin_y,editlcf_w,label_h],...
+    "fontname"            , defaultfont,...
+    "fontunits"           , "points",...
+    "fontsize"            , 12,...
+    "horizontalalignment" , "left", ...
+    "BackgroundColor"     , [1 1 1], ...
+    "ForegroundColor"     , [1 1 1]*0, ...
+    "callback"            , "wfirValue2Sliderpos()")
+    Lcfs=uicontrol( ...
+    "parent"              , gui,...
+    "style"               , "slider",...
+    "visible"             , visible,...
+    "min"                 , 0,...
+    "max"                 , 500,...
+    "sliderstep"          , [1 10],...
+    "value"               , low*fact,...
+    "units"               , "pixels",...
+    "position"            , [xlcf+editlcf_w+margin_x,ylcf-label_h-margin_y,scroll_w,label_h],...
+    "fontname"            , defaultfont,...
+    "fontunits"           , "points",...
+    "fontsize"            , 12,...
+    "horizontalalignment" , "left", ...
+    "BackgroundColor"     , [1 1 1], ...
+    "ForegroundColor"     , [1 1 1]*0, ...
+    "userdata"            , list(Lcfv,fact), ...
+    "callback"            , "wfirSliderpos2Value()")
+    set(Lcfv,"userdata",Lcfs);
+
+    //High cutoff frequency -----------------------------------------------------------------
+    yhcf=ylcf-2*(margin_y+label_h)
+    xhcf=xsf
+    labelhcf_w=labelsf_w;
+    edithcf_w= editsf_w;
+    if ftype=="""lp""" then visible="off";else visible="on";end
+    Hcfl=uicontrol( ...
+    "parent"              , gui,...
+    "style"               , "text",...
+    "visible"             , visible,...
+    "string"              ,  _("High cutoff frequency (Hz)"),...
+    "units"               , "pixels",...
+    "position"            , [xhcf yhcf labelhcf_w label_h],....
+    "fontname"            , defaultfont,...
+    "fontunits"           , "points",...
+    "fontsize"            , 12,...
+    "fontweight"          , "bold", ...
+    "horizontalalignment" , "left", ...
+    "background"          , [1 1 1]*0.8);
+    fact=1000;
+    Hcfv=uicontrol( ...
+    "parent"              , gui,...
+    "style"               , "edit",...
+    "visible"             , visible,...
+    "units"               , "pixels",...
+    "string"              , exprs(6),...
+    "position"            , [xhcf,yhcf-label_h-margin_y,edithcf_w,label_h],...
+    "fontname"            , defaultfont,...
+    "fontunits"           , "points",...
+    "fontsize"            , 12,...
+    "horizontalalignment" , "left", ...
+    "BackgroundColor"     , [1 1 1], ...
+    "ForegroundColor"     , [1 1 1]*0, ...
+    "callback"            , "wfirValue2Sliderpos()")
+    Hcfs=uicontrol( ...
+    "parent"              , gui,...
+    "style"               , "slider",...
+    "visible"             , visible,...
+    "min"                 , 0,...
+    "max"                 , 500,...
+    "sliderstep"          , [1 10],...
+    "value"               , high*fact,...
+    "units"               , "pixels",...
+    "position"            , [xhcf+edithcf_w+margin_x,yhcf-label_h-margin_y,scroll_w,label_h],...
+    "fontname"            , defaultfont,...
+    "fontunits"           , "points",...
+    "fontsize"            , 12,...
+    "horizontalalignment" , "left", ...
+    "BackgroundColor"     , [1 1 1], ...
+    "ForegroundColor"     , [1 1 1]*0, ...
+    "userdata"            , list(Hcfv,fact), ...
+    "callback"            , "wfirSliderpos2Value()")
+    set(Hcfv,"userdata",Hcfs);
+
+
+    //Chebychev or Kaiser window parameter -----------------------------------------------------------------
+    yfp=yhcf-2*(margin_y+label_h)
+    xfp=xsf
+    labelfp_w=labelsf_w;
+    editfp_w= editsf_w;
+    if and(wtype<>["""kr""" """ch"""]) then visible="off";else visible="on";end
+    Fpl=uicontrol( ...
+    "parent"              , gui,...
+    "style"               , "text",...
+    "visible"             , visible,...
+    "string"              ,  _("Parameter"),...
+    "units"               , "pixels",...
+    "position"            , [xfp yfp labelfp_w label_h],....
+    "fontname"            , defaultfont,...
+    "fontunits"           , "points",...
+    "fontsize"            , 12,...
+    "fontweight"          , "bold", ...
+    "horizontalalignment" , "left", ...
+    "background"          , [1 1 1]*0.8);
+    fact=1000;
+
+    Fpv=uicontrol( ...
+    "parent"              , gui,...
+    "style"               , "edit",...
+    "visible"             , visible,...
+    "units"               , "pixels",...
+    "string"              , exprs(7),...
+    "position"            , [xfp,yfp-label_h-margin_y,editfp_w,label_h],...
+    "fontname"            , defaultfont,...
+    "fontunits"           , "points",...
+    "fontsize"            , 12,...
+    "horizontalalignment" , "left", ...
+    "BackgroundColor"     , [1 1 1], ...
+    "ForegroundColor"     , [1 1 1]*0, ...
+    "callback"            , "wfirValue2Sliderpos()")
+    Fps=uicontrol( ...
+    "parent"              , gui,...
+    "style"               , "slider",...
+    "visible"             , visible,...
+    "min"                 , 0,...
+    "max"                 , 500,...
+    "sliderstep"          , [1 10],...
+    "value"               , fp*fact,...
+    "units"               , "pixels",...
+    "position"            , [xfp+editfp_w+margin_x,yfp-label_h-margin_y,scroll_w,label_h],...
+    "fontname"            , defaultfont,...
+    "fontunits"           , "points",...
+    "fontsize"            , 12,...
+    "horizontalalignment" , "left", ...
+    "BackgroundColor"     , [1 1 1], ...
+    "ForegroundColor"     , [1 1 1]*0, ...
+    "userdata"            , list(Fpv,fact), ...
+    "callback"            , "wfirSliderpos2Value()")
+
+    set(Fpv,"userdata",Fps);
+
+
+    set(Filtertype,"userdata",list([Lcfl Lcfv Lcfs;Hcfl Hcfv Hcfs],Ftv))
+    set(Windowtype,"userdata",list([Fpl Fpv Fps], Wtv))
+
+    //Display -----------------------------------------------------------------
+    yd=yfp-2*(margin_y+label_h)
+    xd=xsf
+    labeld_w=40
+    uicontrol( ...
+    "parent"              , gui,...
+    "style"               , "text",...
+    "string"              ,  _("View"),...
+    "units"               , "pixels",...
+    "position"            , [xd yd labeld_w label_h],...
+    "fontname"            , defaultfont,...
+    "fontunits"           , "points",...
+    "fontsize"            , 12,...
+    "fontweight"          , "bold", ...
+    "horizontalalignment" , "left", ...
+    "background"          , [1 1 1]*0.8, ...
+    "visible"             , "on");
+    Fview=uicontrol( ...
+    "parent"              , gui,...
+    "style"               , "checkbox",...
+    "units"               , "pixels",...
+    "min"                 , 0, ...
+    "max"                 , 1, ...
+    "value"               , 0, ...
+    "position"            , [xd+margin_x+labeld_w yd label_h label_h],...
+    "background"          , [1 1 1]*0.8, ...
+    "callback"            , "wfirSetFilterView()",...
+    "visible"             , "on");
+
+    ybuttons=margin_y
+    xbuttons=xhcf
+    Ok = uicontrol( ...
+    "parent"              , gui,...
+    "relief"              , "groove",...
+    "style"               , "pushbutton",...
+    "string"              , _("Ok"),...
+    "units"               , "pixels",...
+    "position"            , [xbuttons,ybuttons,button_w button_h],..
+    "fontname"            , "arial",...
+    "fontunits"           , "points",...
+    "fontsize"            , 12,...
+    "fontweight"          , "bold", ...
+    "horizontalalignment" , "center", ...
+    "visible"             , "on", ...
+    "callback"            , "global ret;ret=1;");
+    xbuttons=xbuttons+margin_x+button_w
+    Cancel = uicontrol( ...
+    "parent"              , gui,...
+    "relief"              , "groove",...
+    "style"               , "pushbutton",...
+    "string"              , _("Cancel"),...
+    "units"               , "pixels",...
+    "position"            , [xbuttons,ybuttons,button_w button_h],..
+    "fontname"            , "arial",...
+    "fontunits"           , "points",...
+    "fontsize"            , 12,...
+    "fontweight"          , "bold", ...
+    "horizontalalignment" , "center", ...
+    "visible"             , "on", ...
+    "callback"            , "global ret;ret=2;");
+    xbuttons=xbuttons+margin_x+button_w
+    Help = uicontrol( ...
+    "parent"              , gui,...
+    "relief"              , "groove",...
+    "style"               , "pushbutton",...
+    "string"              , _("Help"),...
+    "units"               , "pixels",...
+    "position"            , [xbuttons,ybuttons,button_w button_h],..
+    "fontname"            , "arial",...
+    "fontunits"           , "points",...
+    "fontsize"            , 12,...
+    "fontweight"          , "bold", ...
+    "horizontalalignment" , "center", ...
+    "visible"             , "on", ...
+    "callback"            , "help wfir");
+    // next used by wfirGetFilterParameters
+    set(gui,"userdata",[Fview,Filtertype,Windowtype,Forderv,Forders,Lcfv,Lcfs,Hcfv,Hcfs,Fpv,Fps,Sfreq,Ftv,Wtv])
+    realtimeinit(0.1);
+    t=0;
+    while ret==0&or(winsid()==fig_id) then
+        //slow down the waiting loop
+        t=t+1;
+        realtime(t);
+    end
+    if ret==1&or(winsid()==fig_id) then
+        ok=%t
+        [ftype,forder,low,high,wtype,fpar,freq_ech]=wfirGetFilterParameters(gui.userdata)
+        values=tlist(["wfir","ftype","forder","low","high","wtype","fpar","freq_ech"],ftype,forder,low,high,wtype,fpar,freq_ech)
+        exprs= wfirGetFilterExprs(gui.userdata)
+        delete(fig)
+    else
+        //user had canceled or closed the gui window
+        ok=%f
+        values=[]
+        exprs=[]
+    end
+    clearglobal ret idle;
+
+endfunction
+
+function wfirSliderpos2Value(hs)
+    //hs handle on the slider
+    if argn(2)<1 then hs=gcbo,end
+    ud=hs.userdata
+    hv=ud(1);//handle on edition zone
+    fact=ud(2)// scale factor
+    hv.string=string(hs.value/fact)
+    wfirUpdate(hs)
+endfunction
+function wfirSliderpos2ValueI(hs)
+    //hs handle on the slider
+    if argn(2)<1 then hs=gcbo,end
+    hs.value=round(hs.value)
+    ud=hs.userdata
+    hv=ud(1);//handle on edition zone
+    fact=ud(2)// scale factor
+    hv.string=string(hs.value/fact)
+    wfirUpdate(hs)
+endfunction
+
+function wfirValue2Sliderpos(hv)
+    // hv handle on edition zone
+    if argn(2)<1 then hv=gcbo,end
+    if execstr("v="+hv.string,"errcatch")==0 then
+        hs=hv.userdata //handle on the slider
+        ud=hs.userdata
+        fact=ud(2)// scale factor
+        hs.value=fact*v
+        hv.ForegroundColor=[0 0 0];
+        wfirUpdate(hv)
+    else
+        hv.ForegroundColor=[1 0 0];
+    end
+endfunction
+
+function wfirSetFilterType()
+    FT=["lp","hp","bp","sb"]
+    ud=gcbo.userdata
+
+    if gcbo.style=="edit" then
+        he=gcbo //handle on the edit area
+        hl=ud(2) //handle on the list area
+
+        ok=execstr("k=find("+he.string+"==FT)","errcatch")==0
+        if ~ok|k==[] then
+            he.ForegroundColor=[1 0 0];
+            return
+        else
+            he.ForegroundColor=[0 0 0];
+            hl.value=k
+        end
+    else
+        hl=gcbo //handle on the list area
+        he=ud(2) //handle on the edit area
+        he.string=sci2exp(FT(hl.value),0)
+    end
+    ud=hl.userdata
+    Handles=ud(1)
+    select hl.value
+    case 1 then //low pass
+        Handles(1,:).visible="on";
+        Handles(2,:).visible="off";
+    case 2 then //high pass
+        Handles(1,:).visible="off";
+        Handles(2,:).visible="on";
+    case 3 then //band pass
+        Handles(1,:).visible="on";
+        Handles(2,:).visible="on";
+    case 4 then //stop band
+        Handles(1,:).visible="on";
+        Handles(2,:).visible="on";
+    end
+    wfirUpdate(gcbo)
+endfunction
+
+function wfirSetWindowType()
+    WT=["re","tr","hn","hm","kr","ch","ch"]
+    ud=gcbo.userdata
+    if gcbo.style=="edit" then
+        he=gcbo //handle on the edit area
+        hl=ud(2) //handle on the list area
+        ok=execstr("k=find("+he.string+"==WT)","errcatch")==0
+
+        if ~ok|k==[] then
+            he.ForegroundColor=[1 0 0];
+            return
+        else
+            he.ForegroundColor=[0 0 0];
+            hl.value=k
+        end
+    else
+        hl=gcbo //handle on the list area
+        he=ud(2) //handle on the edit area
+        he.string=sci2exp(WT(hl.value),0)
+    end
+
+    ud=hl.userdata
+    Handles=ud(1)
+    k=hl.value
+    if and(k<>[5 6 7]) then
+        Handles.visible="off"
+    else
+        label=Handles(1)
+        value=Handles(2)
+        scroll=Handles(3)
+
+        if k==5 then //Kaiser
+            label.string="Beta"
+            label.visible="on";
+            value.visible="on";
+
+            init=2.5;fact=100;
+            value.string=string(init)
+            scroll.min=0
+            scroll.max=10*fact
+            scroll.value=init*fact
+            d=scroll.userdata;d(2)=fact
+            set(scroll,"userdata",d)
+            scroll.visible="on";
+        elseif k==6 then //Chebychev, main lobe constraint df in [0 0.5]
+            label.string=_("Window main lobe width")
+            label.visible="on";
+            value.visible="on";
+            init=0.2;fact=1000;
+            value.string=string(init)
+            scroll.min=1d-10*fact
+            scroll.max=0.5*fact
+            scroll.value=init*fact
+            d=scroll.userdata;d(2)=fact
+            set(scroll,"userdata",d)
+            scroll.visible="on";
+
+        elseif k==7 then //Chebychev, side lobe constraint dp>0
+            label.string=_("maximum side-lobe height")
+            label.visible="on";
+            value.visible="on";
+            init=0.001;fact=1000;
+            value.string=string(init)
+            scroll.min=1d-10*fact
+            scroll.max=1*fact
+            scroll.value=init*fact
+            d=scroll.userdata;d(2)=fact
+            set(scroll,"userdata",d)
+            scroll.visible="on";
+        end
+    end
+    wfirUpdate(gcbo)
+endfunction
+
+function wfirSetFilterView()
+    h=gcbo //handle on the view check box
+    if h.value==0 //disable
+        fig=h.userdata
+        execstr("delete(fig)","errcatch")
+        h.userdata=[]
+    else //enable
+        nwin=size(winsid(),"*")
+        if nwin==0 then win=1;else win=max(winsid())+1;end
+        fig=scf(win)
+        a=gca();
+        a.x_label.text=_("frequency (Hz)")
+        a.y_label.text=_("Magnitude (dB)")
+        a.axes_visible="on"
+        a.grid=ones(1,2)*color("gray");
+        xpoly([],[])
+        set(h,"userdata",fig)
+        clearglobal idle
+        wfirUpdate(gcbo)
+    end
+endfunction
+
+function wfirSetSamplingFrequency()
+    Sfreq=gcbo
+    ok=execstr("newfreq="+Sfreq.string,"errcatch")==0
+    ok=ok&newfreq>0
+
+    if ok then
+        Sfreq.ForegroundColor=[0 0 0];
+        curfreq=Sfreq.userdata
+        gui=gcbo.parent
+        H=gui.userdata
+        //H=[Fview,Filtertype,Windowtype,Forderv,Forders,Lcfv,Lcfs,Hcfv,Hcfs,Fpv,Fps,Sfreq,Ftv,Wtv]
+        //update low and high cutoff frequencies
+        r=newfreq/curfreq
+        Lcfv=H(6),Lcfs=H(7),Hcfv=H(8),Hcfs=H(9)
+
+        Lcfs.min=Lcfs.min*r
+        Lcfs.max=Lcfs.max*r
+        Lcfs.value=Lcfs.value*r
+        Lcfv.string=string(evstr(Lcfv.string)*r)
+
+        Hcfs.min=Hcfs.min*r
+        Hcfs.max=Hcfs.max*r
+        Hcfs.value=Hcfs.value*r
+        Hcfv.string=string(evstr(Hcfv.string)*r)
+
+        Sfreq.userdata=newfreq
+        wfirUpdate(gcbo)
+    else
+        Sfreq.ForegroundColor=[1 0 0];
+    end
+
+endfunction
+
+function [ftype,forder,low,high,wtype,fpar,freq_ech]=wfirGetFilterParameters(H)
+    //  low,high,freq_ech are returned in Hz
+    //H=[Fview,Filtertype,Windowtype,Forderv,Forders,Lcfv,Lcfs,Hcfv,Hcfs,Fpv,Fps,Sfreq,Ftv,Wtv]
+    FT=["lp","hp","bp","sb"]
+    WT=["re","tr","hn","hm","kr","ch","ch"]
+    Filtertype=H(2);
+    Windowtype=H(3);
+    Forderv=H(4);
+    Lcfv=H(6);
+    Hcfv=H(8);
+    Fpv=H(10);
+    Sfreq=H(12)
+    ftype=FT(Filtertype.value)
+    wtype=WT(Windowtype.value)
+    forder=evstr(Forderv.string)
+    if or(ftype==["hp" "sb"]) then
+        //force odd ordrer
+        if 2*int(forder/2)==forder then forder=forder+1,end
+    end
+    freq_ech=evstr(Sfreq.string)
+    if ftype=="lp" then
+        low=evstr(Lcfv.string)
+        high=0
+    elseif ftype=="hp" then
+        low=evstr(Hcfv.string)
+        high=0
+    else
+        low=evstr(Lcfv.string)
+        high=evstr(Hcfv.string)
+    end
+    if Windowtype.value==5 //Kaiser
+        fpar=[evstr(Fpv.string) 0]
+    elseif Windowtype.value==6 //Chebychev, main lobe constraint df in [0 0.5]
+        fpar=[-1 evstr(Fpv.string)]
+    elseif Windowtype.value==7  //Chebychev, side lobe constraint dp>0
+        fpar=[evstr(Fpv.string) -1]
+    else
+        fpar=[-1 -1]
+    end
+endfunction
+
+function  wfirUpdate(h)
+    global idle
+    if ~idle then return,end
+    idle=%f
+    gui=h.parent
+    H=gui.userdata
+    //  H=[Fview,Filtertype,Windowtype,Forderv,Forders,Lcfv,Lcfs,Hcfv,Hcfs,Fpv,Fps,Sfreq,Ftv,Wtv]
+    if H(1).value==1 then
+        old=gcf()
+        [ftype,forder,low,high,wtype,fpar,freq_ech]= ...
+        wfirGetFilterParameters(H)
+        cfreq=[low,high]/freq_ech;
+        [filt,wfm,fr]=wfir(ftype,forder,cfreq,wtype,fpar);
+        fig=H(1).userdata
+        if ~is_handle_valid(fig) then
+            //the window has been closed by user
+            fig=scf(max(winsid())+1)
+            a=gca();
+            a.x_label.text=_("frequency (Hz)")
+            a.y_label.text=_("Magnitude (dB)")
+            a.axes_visible="on"
+            a.grid=ones(1,2)*color("gray");
+            xpoly([],[])
+            set(H(1),"userdata",fig)
+        end
+        eps=1d-6;
+        frq=linspace(eps,0.5-eps,10*size(filt,"*"))
+        [frq,repf]=repfreq(syslin("d",poly(filt,"z","c"),1),frq)
+        scf(fig);drawlater()
+        a=gca();p=a.children
+        frq=frq(:)*freq_ech
+        m=20*log10(abs(repf(:)))
+        a.data_bounds=[min(frq) min(m);max(frq) max(m)]
+        p.data=[frq m];
+        drawnow()
+        scf(old)
+        clearglobal idle
+    end
+endfunction
+
+function exprs= wfirGetFilterExprs(H)
+    // H=[Fview,Filtertype,Windowtype,Forderv,Forders,Lcfv,Lcfs,Hcfv,Hcfs,Fpv,Fps,Sfreq,Ftv,Wtv]
+    FT=["lp","hp","bp","sb"]
+    WT=["re","tr","hn","hm","kr","ch","ch"]
+    Filtertype=H(2);
+    Windowtype=H(3);
+    Forderv=H(4);
+    Lcfv=H(6);
+    Hcfv=H(8);
+    Fpv=H(10);
+    Sfreq=H(12)
+    exprs=[Ftv.string;
+    Wtv.string;
+    Forderv.string;
+    Sfreq.string;
+    Lcfv.string
+    Hcfv.string;
+    Fpv.string]
+endfunction
diff --git a/modules/signal_processing/macros/wiener.bin b/modules/signal_processing/macros/wiener.bin
new file mode 100755
index 000000000..a02f26a49
--- /dev/null
+++ b/modules/signal_processing/macros/wiener.bin
diff --git a/modules/signal_processing/macros/wiener.sci b/modules/signal_processing/macros/wiener.sci
new file mode 100755
index 000000000..f6dca1efc
--- /dev/null
+++ b/modules/signal_processing/macros/wiener.sci
@@ -0,0 +1,82 @@
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) INRIA - 1988 - C. Bunks
+//
+// 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.1-en.txt
+
+function [xs,ps,xf,pf]=wiener(y,x0,p0,f,g,h,q,r)
+    //<xs,ps,xf,pf>=wiener(y,x0,p0,f,g,h,q,r)
+    //macro which gives the Wiener estimate using
+    //the forward-backward kalman filter formulation
+    //Input to the macro is:
+    //  f,g,h     :system matrices in the interval [t0,tf]
+    //  q,r       :covariance matrices of dynamics
+    //            :and observation noise
+    //  x0,p0     :initial state estimate and error variance
+    //  y         :observations in the interval [t0,tf]
+    //Output from macro is:
+    //  xs        :Smoothed state estimate
+    //  ps        :Error covariance of smoothed estimate
+    //  xf        :Filtered state estimate
+    //  pf        :Error covariance of filtered estimate
+    //form of y:          [y0,y1,...,yf], and yk is a column m-vector
+    //form matrix inputs: [f0,f1,...,ff], and fk is a nxn matrix
+    //                    [g0,g1,...,gf], and gk is a nxn matrix
+    //                    [h0,h1,...,hf], and hk is a mxn matrix
+    //                    [q0,q1,...,qf], and qk is a nxn matrix
+    //                    [r0,r1,...,rf], and gk is a mxm matrix
+    //form of xs and xf:  [x0,x1,...,xf], and xk is a column n-vector
+    //form of ps and pf:  [p0,p1,...,pf], and pk is a nxn matrix
+    //!
+
+    //obtain the dimensions of xk and yk.
+    //Get the time interval [t0,tf].
+
+    t0=1;
+    [n,x0j]=size(x0);
+    [m,tf]=size(y);
+
+    //run the forward kalman filter taking care to save the
+    //state estimate and forward error covariance matrices
+
+    xf1=x0;
+    pf1=p0;
+    xf=[];
+    pf=[];
+    for k=t0:tf,
+        ind_nk=1+(k-1)*n:k*n;
+        ind_mk=1+(k-1)*m:k*m;
+        yk=y(:,k);
+        fk=f(:,ind_nk);
+        gk=g(:,ind_nk);
+        hk=h(:,ind_nk);
+        qk=q(:,ind_nk);
+        rk=r(:,ind_mk);
+        [x1,p1,x,p]=kalm(yk,x0,p0,fk,gk,hk,qk,rk);
+        xf1=[xf1 x1];
+        pf1=[pf1 p1];
+        xf=[xf x];
+        pf=[pf p];
+        x0=x1;
+        p0=p1;
+    end,
+
+    //run the backward kalman filter to smooth the estimate
+
+    x2=x;
+    p2=p;
+    xs=x2;
+    ps=p2;
+    for k=tf-1:-1:t0,
+        ind_nk=1+(k-1)*n:k*n;
+        ind_nk1=1+k*n:(k+1)*n;
+        ak=pf(:,ind_nk)*f(:,ind_nk)'*pf1(:,ind_nk1)**(-1);
+        x2=xf(:,k)+ak*(x2-xf1(:,k+1));
+        p2=pf(:,ind_nk)+ak*(p2-pf1(:,ind_nk1))*ak';
+        xs=[x2 xs];
+        ps=[p2 ps];
+    end
+endfunction
diff --git a/modules/signal_processing/macros/wigner.bin b/modules/signal_processing/macros/wigner.bin
new file mode 100755
index 000000000..f17135e46
--- /dev/null
+++ b/modules/signal_processing/macros/wigner.bin
diff --git a/modules/signal_processing/macros/wigner.sci b/modules/signal_processing/macros/wigner.sci
new file mode 100755
index 000000000..51cd7e7a8
--- /dev/null
+++ b/modules/signal_processing/macros/wigner.sci
@@ -0,0 +1,44 @@
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) INRIA
+//
+// 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.1-en.txt
+
+function tab=wigner(x,h,deltat,zp)
+    //macro which computes the 'time-frequency' wigner
+    //spectrum of a signal.
+    //
+    //	tab    : wigner spectrum (lines correspond to the time variable)
+    //	x      : analysed signal
+    //	h      : data window
+    //      deltat : analysis time increment (in samples)
+    //      zp     : length of FFT's. %pi/zp gives the frequency increment.
+    //!
+    //
+    //   Initializations
+    //
+    l=prod(size(x));
+    n=prod(size(h));
+    npr=2*n;
+    h=h.*conj(h);tab=[];
+    //
+    //   Analytical signal computation using Hilbert transform
+    //
+    [y,y1]=convol(hilb(127),x);
+    z=x+%i*y;
+    //
+    //   Wigner distribution computation
+    //
+    t=n;
+    while t<=l-n,
+        z1=h.*z(t:t+n-1).*conj(z(t:-1:t-n+1));
+        z1(zp)=0;
+        w=fft(z1,-1);
+        w=2*(2*real(w)-z(t)*z(t)'*ones(w));tab=[tab;w];
+        t=t+deltat;
+    end;
+    tab=tab(:,1:zp);
+endfunction
diff --git a/modules/signal_processing/macros/window.bin b/modules/signal_processing/macros/window.bin
new file mode 100755
index 000000000..6a0687fc9
--- /dev/null
+++ b/modules/signal_processing/macros/window.bin
diff --git a/modules/signal_processing/macros/window.sci b/modules/signal_processing/macros/window.sci
new file mode 100755
index 000000000..a471bc82e
--- /dev/null
+++ b/modules/signal_processing/macros/window.sci
@@ -0,0 +1,148 @@
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) INRIA - 1988 - C. Bunks
+//
+// 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.1-en.txt
+
+function [win_l,cwp]=window(wtype,n,par)
+    //[win_l,cwp]=window(wtype,n,par)
+    //macro which calculates symmetric window
+    // wtype :window type (re,tr,hn,hm,kr,ch)
+    // n     :window length
+    // par   :parameter 2-vector (kaiser window: par(1)=beta>0)
+    //       :                   (chebyshev window:par=[dp,df])
+    //       :                   dp=main lobe width (0<dp<.5)
+    //       :                   df=side lobe height (df>0)
+    // win   :window
+    // cwp   :unspecified Chebyshev window parameter
+    //!
+    WT=["re","tr","hm","hn","kr","ch"]
+
+    if and(wtype<>WT) then
+        error(msprintf(_("%s: Wrong value for input argument #%d: Must be in the set {%s}.\n"),"window",1,strcat(WT,",")))
+    end
+    if type(n)<>1|size(n,"*")<>1|~isreal(n)|size(n,"*")<>1|int(n)<>n|n<1 then
+        error(msprintf(_("%s: Wrong type for input argument #%d: A positive integer expected.\n"),"window",2))
+    end
+
+    if or(wtype==["kr","ch"]) then
+        if type(par)<>1|~isreal(par) then
+            error(msprintf(_("%s: Wrong type for input argument #%d: A %d-by-%d real vector expected.\n"),"window",3,1,2))
+        end
+        if wtype=="kr" then
+            if size(par,"*")==0| size(par,"*")>2 then
+                error(msprintf(_("%s: Wrong size for input argument #%d: A %d-by-%d real vector expected.\n"),"window",3,1,2))
+            end
+            Beta=par(1)
+            if Beta<0 then
+                error(msprintf(_("%s: Input argument #%d must be strictly positive.\n"),"window",3))
+            end
+        else //chebychev
+            if size(par,"*")<>2 then
+                error(msprintf(_("%s: Wrong size for input argument #%d: A %d-by-%d real vector expected.\n"),"window",3,1,2))
+            end
+            dp=par(1);df=par(2)
+            if dp>0 then
+                if df>0 then
+                    error(msprintf(_("%s: Wrong value for input argument #%d: incorrect element #%d\n"),"window",3,2))
+                end
+                if dp>=0.5 then
+                    error(msprintf(_("%s: Wrong value for input argument #%d: incorrect element #%d\n"),"window",3,1))
+                end
+                unknown="df";
+            else
+                if df<=0 then
+                    error(msprintf(_("%s: Wrong value for input argument #%d: incorrect element #%d\n"),"window",3,2))
+                end
+                unknown="dp";
+            end
+        end
+    end
+    cwp=-1;
+
+    //Pre-calculations
+
+    no2=(n-1)/2;
+    xt=(-no2:no2);
+    un=ones(1,n);
+
+    //Select the window type
+
+    select wtype
+    case "re" then           //Rectangular window.
+        win_l=un
+    case "tr" then           //Triangular window.
+        win_l=un-2*abs(xt)/(n+1);
+    case "hm" then           //Hamming window.
+        win_l=.54*un+.46*cos(2*%pi*xt/(n-1));
+    case "hn" then           //Hanning window.
+        win_l=.5*un+.5*cos(2*%pi*xt/(n-1));
+    case "kr" then           //Kaiser window with parameter beta (n,beta)
+        //http://en.wikipedia.org/wiki/Kaiser_window
+        win_l = besseli(0,Beta*sqrt(1-(2*(0:n-1)/(n-1)-1).^2))/besseli(0,Beta);
+    case "ch" then           //Chebyshev window
+        m=(n-1)/2;
+        select unknown
+        case "dp" then,
+            dp=1/cosh((n-1)*acosh(1/cos(%pi*df)));
+            cwp=dp;
+        case "df" then
+            df=real(acos(1/cosh(acosh((1+dp)/dp)/(n-1)))/%pi)
+            cwp=df;
+        end
+
+        //Pre-calculation
+
+        np1=int((n+1)/2);
+        odd=modulo(n,2)==1
+
+        x0=(3-cos(2*%pi*df))/(1+cos(2*%pi*df));
+        alpha=(x0+1)/2;
+        Beta=(x0-1)/2;
+        c2=(n-1)/2;
+
+        //Obtain the frequency values of the Chebyshev window
+        f=(0:n-1)/n;
+        xarg=alpha*cos(2*%pi*f)+Beta;
+
+        //Evaluate Chebyshev polynomials
+        //           / cos(n acos(x),    |x| <= 1
+        //    p(x) = |
+        //           \ cosh(n acosh(x),  |x| > 1
+        pr=zeros(1,n);//real part
+        pi=zeros(1,n); //imaginary part
+        ind = find(xarg<=1); pr(ind)=dp*cos (c2*acos (xarg(ind)));
+        ind = find(xarg>1);  pr(ind)=dp*cosh(c2*acosh(xarg(ind)));
+        pr=real(pr);
+
+        if ~odd then
+            pr=pr.*cos(%pi*f);
+            pi=-pr.*sin(%pi*f);
+            antisym=[1*ones(1:int(n/2)+1),-1*ones(int(n/2)+2:n)];
+            pr=pr.*antisym;
+            pi=pi.*antisym;
+        end
+
+        //Calculate the window coefficients using the inverse DFT
+        twn=2*%pi/n;
+        xj=(0:n-1);
+        for i=1:np1;
+            arg=twn*(i-1)*xj
+            w(i)=sum(pr.*cos(arg)+pi.*sin(arg));
+        end,
+        c1=w(1);
+        w=w/c1;
+        if odd then
+            win_l(np1:n)=w(1:np1);
+            win_l(1:np1-1)=w(np1-1:-1:1);
+        else,
+            win_l(np1+1:n)=w(1:np1);
+            win_l(1:np1)=w(np1:-1:1);
+        end
+        win_l=real(win_l');
+    end
+
+endfunction
diff --git a/modules/signal_processing/macros/xcorr.bin b/modules/signal_processing/macros/xcorr.bin
new file mode 100755
index 000000000..de0e4c996
--- /dev/null
+++ b/modules/signal_processing/macros/xcorr.bin
diff --git a/modules/signal_processing/macros/xcorr.sci b/modules/signal_processing/macros/xcorr.sci
new file mode 100755
index 000000000..a72263b6c
--- /dev/null
+++ b/modules/signal_processing/macros/xcorr.sci
@@ -0,0 +1,124 @@
+// This file is part Scilab
+// Copyright (C) 2012 - INRIA - Serge Steer
+// 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.1-en.txt
+function [c,lagindex]=xcorr(x,varargin)
+    nv=size(varargin)
+    if nv>0&type(varargin(nv))==10 then
+        validemodes=["biased","unbiased","coeff","none"]
+        scalemode=varargin(nv)
+        if and(scalemode<>validemodes) then
+            error(msprintf(_("%s: Wrong value for input argument #%d: Must be in the set {%s}.\n"),...
+            "xcorr",nv+1,strcat(""""+validemodes+"""",",")))
+        end
+        nv=nv-1;
+    else
+        scalemode="none"
+    end
+    //test de validité de x
+    szx=size(x)
+    if type(x)<>1|and(szx>1) then
+        error(msprintf(_("%s: Wrong type for input argument #%d: Real or complex vector expected.\n"),...
+        "xcorr",1))
+    end
+
+    autocorr=%t
+    maxlags=[]
+    if nv==1 then
+        if size(varargin(1),"*")==1 then //xcorr(x,maxlags)
+            autocorr=%t
+            maxlags=int(varargin(1))
+            if type( maxlags)<>1|size(maxlags,"*")>1|~isreal(maxlags)|maxlags<>int(maxlags) then
+                error(msprintf(_("%s: Wrong type for argument #%d: an integer expected.\n"),...
+                "xcorr",2))
+            end
+            if maxlags<1 then
+                error(msprintf(_("%s: Wrong value for argument #%d: the expected value must be greater than %d.\n"),...
+                "xcorr",2,1))
+            end
+        else //xcorr(x,y)
+            autocorr=%f
+            y=varargin(1)
+            if type(y)<>1|and(size(y)>1) then
+                error(msprintf(_("%s: Wrong type for input argument #%d: Real or complex vector expected.\n"),...
+                "xcorr",2))
+            end
+            maxlags=[]
+        end
+    elseif nv==2 then //xcorr(x,y,maxlag)
+        autocorr=%f
+        y=varargin(1)
+        if type(y)<>1|and(size(y)>1) then
+            error(msprintf(_("%s: Wrong type for input argument #%d: Real or complex vector expected.\n"),...
+            "xcorr",2))
+        end
+        maxlags=int(varargin(2))
+        if type( maxlags)<>1|size(maxlags,"*")>1|~isreal(maxlags)|maxlags<>int(maxlags) then
+            error(msprintf(_("%s: Wrong type for argument #%d: an integer expected.\n"),...
+            "xcorr",2))
+        end
+        if maxlags<1 then
+            error(msprintf(_("%s: Wrong value for argument #%d: the expected value must be greater than %d.\n"),...
+            "xcorr",2,1))
+        end
+    end
+
+    if autocorr then //auto correlation
+        x=matrix(x,-1,1);n=size(x,"*")
+        if maxlags==[] then maxlags=n-1,end
+        npad=2^nextpow2(2*n-1)
+        x(npad)=0;
+        t=fft(x);
+        c=ifft(real(t.*conj(t)))
+        if isreal(x) then c=real(c);end
+    else //cross correlation
+        x=matrix(x,-1,1);nx=size(x,1)
+        xx=sum(abs(x).^2)
+        y=matrix(y,-1,1);ny=size(y,1)
+        yy=sum(abs(y).^2)
+        if nx<ny then
+            x(ny)=0;
+        elseif ny<nx then
+            y(nx)=0;
+        end
+        n=max(nx,ny)
+        if maxlags==[] then maxlags=n-1,end
+        npad=2^nextpow2(2*n-1)
+        x(npad)=0;
+        y(npad)=0;
+        c=ifft(fft(x).*conj(fft(y)))
+        if isreal(x)&isreal(y) then c=real(c),end
+    end
+    //extract requested lags
+    padding=zeros(maxlags-n+1,1)
+    if maxlags<n then
+        c=[c($-maxlags+1:$);c(1:maxlags+1);]
+    else
+        padding=zeros(maxlags-n+1,1)
+        c = [padding;
+        c($-n+2:$);
+        c(1:n);
+        padding];
+    end
+    //normalization
+    select scalemode
+    case "biased" then
+        c=c/n
+    case "unbiased" then
+        scale=n-abs(-maxlags:maxlags)
+        scale(scale==0)=1;
+        c=c./scale'
+    case "coeff" then
+        if autocorr then
+            c=c/c(maxlags+1)
+        else
+            c=c/sqrt(xx*yy)
+        end
+    end
+    //give result same orientation as x
+    if szx(1)==1 then c=matrix(c,1,-1),end
+    if argn(1)==2 then lagindex=-maxlags:maxlags,end
+endfunction
diff --git a/modules/signal_processing/macros/xcov.bin b/modules/signal_processing/macros/xcov.bin
new file mode 100755
index 000000000..abc3a1251
--- /dev/null
+++ b/modules/signal_processing/macros/xcov.bin
diff --git a/modules/signal_processing/macros/xcov.sci b/modules/signal_processing/macros/xcov.sci
new file mode 100755
index 000000000..2a744e97d
--- /dev/null
+++ b/modules/signal_processing/macros/xcov.sci
@@ -0,0 +1,76 @@
+// This file is part Scilab
+// Copyright (C) 2012 - INRIA - Serge Steer
+// 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.1-en.txt
+function [c,lagindex]=xcov(x,varargin)
+    nv=size(varargin)
+    if nv>0&type(varargin(nv))==10 then
+        validemodes=["biased","unbiased","coeff","none"]
+        scalemode=varargin(nv)
+        if and(scalemode<>validemodes) then
+            error(msprintf(_("%s: Wrong value for input argument #%d: Must be in the set {%s}.\n"),...
+            "xcov",nv+1,strcat(""""+validemodes+"""",",")))
+        end
+        nv=nv-1;
+    else
+        scalemode="none"
+    end
+    //test de validité de x
+    szx=size(x)
+    if type(x)<>1|and(szx>1) then
+        error(msprintf(_("%s: Wrong type for input argument #%d: Real or complex vector expected.\n"),...
+        "xcov",1))
+    end
+    x=x-mean(x)
+    autocorr=%t
+    maxlags=[]
+    if nv==1 then
+        if size(varargin(1),"*")==1 then //xcov(x,maxlags)
+            autocorr=%t
+            maxlags=int(varargin(1))
+            if type( maxlags)<>1|size(maxlags,"*")>1|~isreal(maxlags)|maxlags<>int(maxlags) then
+                error(msprintf(_("%s: Wrong type for argument #%d: an integer expected.\n"),...
+                "xcov",2))
+            end
+            if maxlags<1 then
+                error(msprintf(_("%s: Wrong value for argument #%d: the expected value must be greater than %d.\n"),...
+                "xcov",2,1))
+            end
+        else //xcov(x,y)
+            autocorr=%f
+            y=varargin(1)
+            if type(y)<>1|and(size(y)>1) then
+                error(msprintf(_("%s: Wrong type for input argument #%d: Real or complex vector expected.\n"),...
+                "xcov",2))
+            end
+            varargin(1)=y-mean(y)
+            maxlags=[]
+        end
+    elseif nv==2 then //xcov(x,y,maxlag)
+        autocorr=%f
+        y=varargin(1)
+        if type(y)<>1|and(size(y)>1) then
+            error(msprintf(_("%s: Wrong type for input argument #%d: Real or complex vector expected.\n"),...
+            "xcov",2))
+        end
+        if type(y)<>1|and(size(y)>1) then
+            error(msprintf(_("%s: Wrong type for input argument #%d: Real or complex vector expected.\n"),...
+            "xcov",2))
+        end
+        varargin(1)=y-mean(y)
+        maxlags=int(varargin(2))
+        if type( maxlags)<>1|size(maxlags,"*")>1|~isreal(maxlags)|maxlags<>int(maxlags) then
+            error(msprintf(_("%s: Wrong type for argument #%d: an integer expected.\n"),...
+            "xcov",2))
+        end
+        if maxlags<1 then
+            error(msprintf(_("%s: Wrong value for argument #%d: the expected value must be greater than %d.\n"),...
+            "xcov",2,1))
+        end
+    end
+    [c,lagindex]=xcorr(x,varargin(:))
+
+endfunction
diff --git a/modules/signal_processing/macros/yulewalk.bin b/modules/signal_processing/macros/yulewalk.bin
new file mode 100755
index 000000000..615ee791f
--- /dev/null
+++ b/modules/signal_processing/macros/yulewalk.bin
diff --git a/modules/signal_processing/macros/yulewalk.sci b/modules/signal_processing/macros/yulewalk.sci
new file mode 100755
index 000000000..fcccba7db
--- /dev/null
+++ b/modules/signal_processing/macros/yulewalk.sci
@@ -0,0 +1,147 @@
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) INRIA
+//
+// 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.1-en.txt
+
+function [Nz,Dz]=yulewalk(Norder, frq, mag)
+    //YULEWALK  filter design using a least-squares method.
+    //	Hz = YULEWALK(N,frq,mag) finds the N-th order iir filter
+    //
+    //		               n-1         n-2
+    //		   B(z)   b(1)z     + b(2)z    + .... + b(n)
+    //	     H(z)= ---- = ---------------------------------
+    //		            n-1       n-2
+    //		   A(z)    z   + a(1)z    + .... + a(n)
+    //
+    //which matches the magnitude frequency response given by vectors F and M.
+    //Vectors frq and mag specify the frequency and magnitude of the desired
+    //frequency response. The frequencies in frq must be between 0.0 and 1.0,
+    //with 1.0 corresponding to half the sample rate. They must be in
+    //increasing order and start with 0.0 and end with 1.0.
+    //
+    // Example: f=[0,0.4,0.4,0.6,0.6,1];H=[0,0,1,1,0,0];Hz=yulewalk(8,f,H);
+    //fs=1000;fhz = f*fs/2;
+    //clf(0);xset('window',0);plot2d(fhz',H');
+    //xtitle('Desired Frequency Response')
+    //[frq,repf]=repfreq(Hz,0:0.001:0.5);
+    //clf(1);xset('window',1);plot2d(fs*frq',abs(repf'));
+    //xtitle('Obtained Frequency Response')
+    //
+    [LHS,RHS]=argn(0);
+    if RHS <>3
+        error(msprintf(gettext("%s: Wrong number of input argument(s): %d expected.\n"),"yulewalk",3));
+    end
+    npt=512;
+    thelap=fix(npt/25);
+
+    [mf,nf]=size(frq);
+    [mm,nn]=size(mag);
+    if mm ~= mf | nn ~= nf
+        error(msprintf(gettext("%s: Incompatible input arguments #%d and #%d: Same sizes expected.\n"),"yulewalk",2,3));
+    end
+    nbrk=max(mf,nf);
+    if mf < nf
+        frq=frq';
+        mag=mag';
+    end
+
+    if abs(frq(1)) > %eps | abs(frq(nbrk) - 1) > %eps
+        error(msprintf(gettext("%s: Wrong values for input argument #%d: The first element must be %s and the last %s.\n"),"yulewalk",2,"0","1"));
+    end
+
+    npt=npt+1;
+    Ht=zeros(1,npt);
+
+    nint=nbrk-1;
+    df=frq(2:nf)-frq(1:nf-1);
+    if (or(df < 0))
+        error(msprintf(gettext("%s: Wrong values for input argument #%d: Elements must be in increasing order.\n"),"yulewalk",2));
+    end
+
+    nb=1;
+    Ht(1)=mag(1);
+    for i=1:nint
+        if df(i) == 0
+            nb=nb - thelap/2;
+            ne=nb + thelap;
+        else
+            ne=int(frq(i+1)*npt);
+        end
+        if (nb < 0 | ne > npt)
+            error(msprintf(gettext("%s: Too abrupt change near end of frequency range.\n"),"yulewalk"));
+        end
+        j=nb:ne;
+        if ne == nb
+            inc=0;
+        else
+            inc=(j-nb)/(ne-nb);
+        end
+        Ht(nb:ne)=inc*mag(i+1) + (1 - inc)*mag(i);
+        nb=ne + 1;
+    end
+    Ht=[Ht Ht(npt-1:-1:2)];
+    n=max(size(Ht));
+    n2=fix((n+1)/2);
+    nb=Norder;
+    nr=4*Norder;
+    nt=0:1:nr-1;
+    R=real(fft(Ht.*Ht,1));
+    R =R(1:nr).*(0.54+0.46*cos(%pi*nt/(nr-1)));
+    Rwindow=[1/2,ones(1,n2-1),zeros(1,n-n2)];
+
+    nr=max(size(R));
+    Rm=toeplitz(R(Norder+1:nr-1),R(Norder+1:-1:2));
+    Rhs=- R(Norder+2:nr);
+    denf=[1 Rhs/Rm'];
+
+    A=polystab(denf);
+    nA=size(A,"*");
+    Dz=poly(A(nA:-1:1),"z","c");
+    Qh=numf([R(1)/2,R(2:nr)],A,Norder);	// compute additive decomposition
+    Qhsci=poly(Qh( size(Qh,"*"):-1:1 ),"z","c")
+    aa=A(:);bb=Qh(:);
+    nna=max(size(aa));nnb=max(size(bb));
+    Ss=2*real((fft([Qh,zeros(1,n-nnb)],-1) ./ fft([A,zeros(1,n-nna)],-1)));
+    hh=fft(exp(fft( Rwindow.*fft(log(Ss),1),-1)),1);
+    impr=filter(1,A,[1 zeros(1,nr-1)]);
+    B=real(hh(1:nr)/toeplitz(impr,[1 zeros(1,nb)])');
+    nB=size(B,"*");
+    Nz=poly(B(nB:-1:1),"z","c");
+    B =real(numf(hh(1:nr),A,nb));
+    if LHS==1 then
+        Nz=syslin("d",Nz/Dz);
+    end
+endfunction
+
+function b=polystab(a);
+    //Utility function for use with yulewalk: polynomial stabilization.
+    //	polystab(A), where A is a vector of polynomial coefficients,
+    //	stabilizes the polynomial with respect to the unit circle;
+    //	roots whose magnitudes are greater than one are reflected
+    //	inside the unit circle.
+    if length(a) == 1, b=a; return, end
+    ll=size(a,"*");
+    ap=poly(a($:-1:1),"s","coeff");
+    v=roots(ap); i=find(v~=0);
+    vs=0.5*(sign(abs(v(i))-1)+1);
+    v(i)=(1-vs).*v(i) + vs./conj(v(i));
+    b=a(1)*coeff(poly(v,"s"));
+    b =b(ll:-1:1);
+    if ~or(imag(a))
+        b=real(b);
+    end
+
+endfunction
+
+function b=numf(h,a,nb)
+    //NUMF	Find numerator B given impulse-response h of B/A and denominator A
+    //NB is the numerator order.  This function is used by yulewalk.
+
+    nh=max(size(h));
+    impr=filter(1,a,[1 zeros(1,nh-1)]);
+    b=h/toeplitz(impr,[1 zeros(1,nb)])';
+endfunction
diff --git a/modules/signal_processing/macros/zpbutt.bin b/modules/signal_processing/macros/zpbutt.bin
new file mode 100755
index 000000000..c87337723
--- /dev/null
+++ b/modules/signal_processing/macros/zpbutt.bin
diff --git a/modules/signal_processing/macros/zpbutt.sci b/modules/signal_processing/macros/zpbutt.sci
new file mode 100755
index 000000000..65c9af9c1
--- /dev/null
+++ b/modules/signal_processing/macros/zpbutt.sci
@@ -0,0 +1,26 @@
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) INRIA - F.D
+// Copyright (C) INRIA - 1996 - C. Bunks
+//
+// 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.1-en.txt
+
+function [pols,gain]=zpbutt(n,omegac)
+    //<pols,gain>=zpbutt(n,omegac)
+    //Computes the poles of a Butterworth analog
+    //filter of order n and cutoff frequency omegac
+    //transfer function H(s) is calculated by
+    //     H(s) = gain/real(poly(pols,'s'))
+    //  n      :Filter order
+    //  omegac :Cut-off frequency in Hertz rd/s
+    //  pols   :Resulting poles of filter
+    //  gain   :Resulting gain of filter
+    //
+    //!
+    angles=ones(1,n)*(%pi/2+%pi/(2*n))+(0:n-1)*%pi/n;
+    pols=omegac*exp(%i*angles);
+    gain=abs((-omegac)^n);
+endfunction
diff --git a/modules/signal_processing/macros/zpch1.bin b/modules/signal_processing/macros/zpch1.bin
new file mode 100755
index 000000000..33bc7d2cf
--- /dev/null
+++ b/modules/signal_processing/macros/zpch1.bin
diff --git a/modules/signal_processing/macros/zpch1.sci b/modules/signal_processing/macros/zpch1.sci
new file mode 100755
index 000000000..d716101e7
--- /dev/null
+++ b/modules/signal_processing/macros/zpch1.sci
@@ -0,0 +1,34 @@
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) INRIA - F.D
+// Copyright (C) INRIA - 1996 - C. Bunks
+//
+// 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.1-en.txt
+
+function [pols,gain]=zpch1(n,epsilon,omegac)
+    //Poles of a Type 1 Chebyshev analog filter
+    //The transfer function is given by :
+    //H(s)=gain/poly(pols,'s')
+    //  n       :Filter order
+    //  epsilon :Ripple in the pass band (0<epsilon<1)
+    //  omegac  :Cut-off frequency in rd/s
+    //  pols    :Resulting filter poles
+    //  gain    :Resulting filter gain
+    //
+    //!
+    Gamma=((1+sqrt(1+epsilon**2))/epsilon)^(1/n);
+    a=omegac*(Gamma-1/Gamma)/2;
+    b=omegac*(Gamma+1/Gamma)/2;
+    v=%pi/(2*n):%pi/n:(2*n-1)/(2*n)*%pi;
+    sigma=-a*sin(v);
+    omega=b*cos(v);
+    pols=sigma+%i*omega;
+    gain=abs(real(prod(pols)));
+    if n==2*int(n/2) then,
+        gain=gain/sqrt(1+epsilon^2);
+    end
+
+endfunction
diff --git a/modules/signal_processing/macros/zpch2.bin b/modules/signal_processing/macros/zpch2.bin
new file mode 100755
index 000000000..8d9596df6
--- /dev/null
+++ b/modules/signal_processing/macros/zpch2.bin
diff --git a/modules/signal_processing/macros/zpch2.sci b/modules/signal_processing/macros/zpch2.sci
new file mode 100755
index 000000000..2e9677ac4
--- /dev/null
+++ b/modules/signal_processing/macros/zpch2.sci
@@ -0,0 +1,42 @@
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) INRIA - F.D
+// Copyright (C) INRIA - 1996 - C. Bunks
+//
+// 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.1-en.txt
+
+function [zers,pols,gain]=zpch2(n,A,omegar)
+    //[zers,pols,gain]=zpch2(n,A,omegar)
+    //Poles and zeros of a type 2 Chebyshev analog filter
+    //gain is the gain of the filter
+    //H(s)=gain*poly(zers,'s')/poly(pols,'s')
+    //  n      :Filter order
+    //  A      :Attenuation in stop band (A>1)
+    //  omegar :Cut-off frequency in rd/s
+    //  zers   :Resulting filter zeros
+    //  pols   :Resulting filter poles
+    //  gain   :Resulting filter gain
+    //
+    //!
+    un=ones(1,n);
+    v=%pi/(2*n)*(1:2:2*n-1);
+    w=exp(%i*v);
+    cosine=real(w);
+    sine=imag(w);
+    n2=int(n/2);
+    if n==2*n2 then,
+        zers=%i*omegar*un./cosine;
+    else,
+        zers=%i*omegar*un(1:n-1)./[cosine(1:n2),cosine(n2+2:n)];
+    end,
+    Gamma=(A+sqrt(A*A-1))**(1/n);
+    alpha=-((Gamma-1/Gamma)/2)*sine;
+    Beta=((Gamma+1/Gamma)/2)*cosine;
+    normal=alpha.*alpha+Beta.*Beta;
+    pols=omegar*(alpha-%i*Beta)./normal;
+    gain=abs(real(prod(pols)/prod(zers)));
+
+endfunction
diff --git a/modules/signal_processing/macros/zpell.bin b/modules/signal_processing/macros/zpell.bin
new file mode 100755
index 000000000..e4b61696c
--- /dev/null
+++ b/modules/signal_processing/macros/zpell.bin
diff --git a/modules/signal_processing/macros/zpell.sci b/modules/signal_processing/macros/zpell.sci
new file mode 100755
index 000000000..ef3b90871
--- /dev/null
+++ b/modules/signal_processing/macros/zpell.sci
@@ -0,0 +1,55 @@
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) INRIA - 1989 - F.Delebecque
+// Copyright (C) INRIA - 1996 - C. Bunks
+//
+// 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.1-en.txt
+
+function [ze,po,gain]=zpell(epsilon,A,omegac,omegar)
+    //[ze,po,gain]=zpell(epsilon,A,omegac,omegar)
+    //Poles and zeros of prototype lowpass elliptic filter
+    //gain is the gain of the filter
+    //  epsilon :Ripple of filter in pass band (0<epsilon<1)
+    //  A       :Attenuation of filter in stop band (A>1)
+    //  omegac  :Pass band cut-off frequency in rd/s
+    //  omegar  :Stop band cut-off frequency in rd/s
+    //  ze      :Resulting zeros of filter
+    //  po      :Resulting poles of filter
+    //  gain    :Resulting gain of filter
+    //
+    //!
+
+    m1=(epsilon*epsilon)/(A*A-1);
+    K1=delip(1,sqrt(m1));
+    K1t=imag(delip(1/sqrt(m1),sqrt(m1)));
+    m=(omegac/omegar)^2;
+    K=delip(1,sqrt(m));
+    Kt=imag(delip(1/sqrt(m),sqrt(m)));
+    n=(K1t*K)/(K1*Kt);
+    order=round(n);
+    u0=-(Kt/K1t)*delip(sqrt(1/(1+epsilon*epsilon)),sqrt(1-m1));
+    even=2*int(order/2);
+    if order<>even then,
+        vmin=2*K/n;
+    else,
+        vmin=K/n;
+    end,
+    v=vmin:(2*K/n):K;
+    un=ones(1:max(size(v)));
+    zlambda=-un*Kt+%i*v;
+    plambda= u0*un+%i*v;
+    ze=%i*imag(%i*omegac*%sn(-%i*zlambda,m));
+    ze=[ze,conj(ze)];
+    po=%i*omegac*%sn(-%i*plambda,m);
+    po=[po,conj(po)];
+    if order<>even then,
+        po=[po,%i*omegac*%sn(-%i*u0,m)];
+    end,
+    gain=abs(real(prod(po))/real(prod(ze)));
+    if order==even then,
+        gain=gain/sqrt(1+epsilon^2);
+    end;
+endfunction