corr("u", ...) 18/15618/12
Anais AUBERT [Tue, 2 Dec 2014 09:16:14 +0000 (10:16 +0100)]
rand('normal');x=rand(1,256);y=-x;
[c,mxy]=corr(x,y,32);
x=x-mxy(1)*ones(x);y=y-mxy(2)*ones(y);  //centring
x1=x(1:128);x2=x(129:256);
y1=y(1:128);y2=y(129:256);
w0=0*ones(1:64);   //32 coeffs
[w1,xu]=corr('u',x1,y1,w0)

Change-Id: Ibb5ba490bf3f2665bf795dc723a07b7c60c22fe5

scilab/modules/signal_processing/sci_gateway/cpp/sci_corr.cpp
scilab/modules/signal_processing/src/cpp/signalprocessingfunctions.hxx

index aca2cdd..fe765bb 100644 (file)
@@ -2,6 +2,7 @@
 * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
 * Copyright (C) 2011 - DIGITEO - Antoine ELIAS
 * Copyright (C) 2014 - Scilab Enterprises - Anais AUBERT
+* Copyright (C) 2014 - Scilab Enterprises - Sylvain GENIN
 *
 * This file must be used under the terms of the CeCILL.
 * This source file is licensed as described in the file COPYING, which
@@ -169,7 +170,219 @@ types::Function::ReturnValue sci_corr(types::typed_list &in, int _iRetCount, typ
         }
         else if (pS->getSize() == 1 && pS->get(0)[0] == L'u')
         {
-            //update
+            types::Double* pDblIn1 = NULL;
+            types::Double* pDblIn2 = NULL;
+            types::Double* pDblIn3 = NULL;
+            types::Double* pDblIn4 = NULL;
+
+            int iErr = 0;
+            int mnx = 0;
+            int mny = 0;
+            int mfft = 0;
+            int nbx = 0;
+            int ichaud = 0;
+            int iMode = 0;
+
+            double* x = NULL;
+            double* xu = NULL;
+            double* xui = NULL;
+            double* w = NULL;
+            double* wi = NULL;
+            double* y = NULL;
+            double* yi = NULL;
+
+            if (in[1]->isDouble() == false)
+            {
+                Scierror(999, _("%s: Wrong type for input argument #%d: Matrix expected.\n"), "corr" , 2);
+                return types::Function::Error;
+            }
+
+            pDblIn1 = in[1]->getAs<types::Double>();
+            if (pDblIn1->isComplex())
+            {
+                Scierror(999, _("%s: Wrong type for input argument #%d: Real matrix expected.\n"), "corr" , 2);
+                return types::Function::Error;
+            }
+
+            mnx = pDblIn1->getRows() * pDblIn1->getCols();
+
+            x = pDblIn1->get();
+
+
+            if (in[2]->isDouble() == false)
+            {
+                Scierror(999, _("%s: Wrong type for input argument #%d: Matrix expected.\n"), "corr" , 3);
+                return types::Function::Error;
+            }
+
+            pDblIn2 = in[2]->getAs<types::Double>();
+
+            mny = pDblIn2->getRows() * pDblIn2->getCols();
+
+            if (mnx == mny)
+            {
+                iMode = 1;
+                if (pDblIn2->isComplex())
+                {
+                    Scierror(999, _("%s: Wrong type for input argument #%d: Real matrix expected.\n"), "corr" , 3);
+                    return types::Function::Error;
+                }
+
+                y = pDblIn2->get();
+            }
+
+            if (iMode == 0)
+            {
+                mfft = mny;
+                if (pDblIn2->isComplex() == false)
+                {
+                    double* wtempo = NULL;
+                    w = new double[pDblIn2->getSize()];
+                    wi = new double[mfft];
+                    memset(wi, 0x00, sizeof(double) * mfft);
+
+                    wtempo = pDblIn2->get();
+                    memcpy(w, wtempo, sizeof(double) * pDblIn2->getSize());
+                }
+                else
+                {
+                    double* wtempo = NULL;
+                    double* witempo = NULL;
+                    w = new double[pDblIn2->getSize()];
+                    wi = new double[pDblIn2->getSize()];
+
+                    wtempo = pDblIn2->getReal();
+                    witempo = pDblIn2->getImg();
+
+                    memcpy(w, wtempo, sizeof(double) * pDblIn2->getSize());
+                    memcpy(wi, witempo, sizeof(double) * pDblIn2->getSize());
+                }
+
+
+                if (in.size() == 4)
+                {
+                    pDblIn3 = in[3]->getAs<types::Double>();
+                    if (pDblIn3->isComplex())
+                    {
+                        Scierror(999, _("%s: Wrong type for input argument #%d: Real matrix expected.\n"), "corr" , 4);
+                        return types::Function::Error;
+                    }
+
+                    xui = new double[mfft * 2];
+                    double* xutempo = NULL;
+                    xutempo = pDblIn3->get();
+                    xu = new double[mfft * 2];
+                    memset(xu, 0x00, sizeof(double) * mfft * 2);
+                    memcpy(xu, xutempo, sizeof(double) * pDblIn3->getSize());
+
+                    nbx =  pDblIn3->getSize();
+                    ichaud = 1;
+                }
+                else
+                {
+                    xu = new double[mfft * 2];
+                    xui = new double[mfft * 2];
+                }
+
+                yi = new double[mny];
+                C2F(cmpse3)(&mfft, &mnx, &iMode, x, yi, xu, xui, w, wi, &iErr, &ichaud, &nbx);
+                if (iErr > 0)
+                {
+                    delete[] xu;
+                    delete[] xui;
+                    delete[] wi;
+                    delete[] w;
+                    Scierror(999, _("fft call : needs power of two!"));
+                    return types::Function::Error;
+                }
+
+            }
+            else
+            {
+                pDblIn3 = in[3]->getAs<types::Double>();
+                mfft  =   pDblIn3->getRows() * pDblIn3->getCols();
+                if (pDblIn3->isComplex() == false)
+                {
+                    wi = new double[mfft];
+                    memset(wi, 0x00, sizeof(double) * mfft);
+
+                    w = new double[pDblIn3->getSize()];
+                    double* wtempo = NULL;
+                    wtempo = pDblIn3->get();
+                    memcpy(w, wtempo, sizeof(double) * pDblIn3->getSize());
+
+                }
+                else
+                {
+                    double* wtempo = NULL;
+                    double* witempo = NULL;
+                    w = new double[pDblIn3->getSize()];
+                    wi = new double[pDblIn3->getSize()];
+
+                    wtempo = pDblIn3->getReal();
+                    witempo = pDblIn3->getImg();
+
+                    memcpy(w, wtempo, sizeof(double) * pDblIn3->getSize());
+                    memcpy(wi, witempo, sizeof(double) * pDblIn3->getSize());
+                }
+                if (in.size() == 5)
+                {
+                    pDblIn4 = in[4]->getAs<types::Double>();
+                    nbx = pDblIn4->getSize();
+                    double* xutempo = NULL;
+                    xutempo = pDblIn4->get();
+                    xu = new double[mfft * 2];
+                    memset(xu, 0x00, sizeof(double) * mfft * 2);
+                    memcpy(xu, xutempo, sizeof(double) * pDblIn4->getSize());
+                    ichaud = 1;
+
+                    xui = new double[mfft * 2];
+                }
+                else
+                {
+                    xu = new double[mfft * 2];
+                    xui = new double[mfft * 2];
+                }
+
+                C2F(cmpse3)(&mfft, &mnx, &iMode, x, y, xu, xui, w, wi, &iErr, &ichaud, &nbx);
+                if (iErr > 0)
+                {
+                    delete[] xu;
+                    delete[] xui;
+                    delete[] wi;
+                    delete[] w;
+                    Scierror(999, _("fft call : needs power of two!"));
+                    return types::Function::Error;
+                }
+
+            }
+
+            types::Double *pDblOut1 = NULL;
+            pDblOut1 = new types::Double(1, mfft, true);
+            pDblOut1->set(w);
+            pDblOut1->setImg(wi);
+            out.push_back(pDblOut1);
+
+            if (_iRetCount == 2)
+            {
+                types::Double *pDblOut2 = NULL;
+                pDblOut2 = new types::Double(1, mfft / 2);
+
+                for (int i = 0; i < mfft / 2; i++)
+                {
+                    xui[i] = x[mnx - mfft / 2 + i];
+                }
+
+                pDblOut2->set(xui);
+                out.push_back(pDblOut2);
+
+            }
+            delete[] w;
+            delete[] wi;
+            delete[] xui;
+            delete[] xu;
+            return types::Function::OK;
+
         }
         else
         {
index 0595c70..a7c09ef 100644 (file)
@@ -29,6 +29,8 @@ extern "C"
 {
     extern void C2F(tscccf)(double *x, double *y, int *length, double *cxy, double *xymean, int *lag, int *error);
     extern void C2F(cmpse2)(int *iSect, int *iTotalSize, int *iMode, void *pXFunction, void *pYFunction, double *xa, double *xr, double *xi, double *zr, double *zi, int *error);
+
+    extern void C2F(cmpse3)(int *mfft, int *mnx, int *iMode, double *x, double *yi, double *xu, double *xui, double * w, double *wi, int *iErr, int *ichaud, int *nbx);
 }
 
 typedef void(*dgetx_f_t)(double*, int*, int*);