scilab/modules/signal_processing/sci_gateway/cpp/sci_corr.cpp

   1 /*
   2 * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
   3 * Copyright (C) 2011 - DIGITEO - Antoine ELIAS
   4 * Copyright (C) 2014 - Scilab Enterprises - Anais AUBERT
   5 * Copyright (C) 2014 - Scilab Enterprises - Sylvain GENIN
   6 *
   7  * Copyright (C) 2012 - 2016 - Scilab Enterprises
   8  *
   9  * This file is hereby licensed under the terms of the GNU GPL v2.0,
  10  * pursuant to article 5.3.4 of the CeCILL v.2.1.
  11  * This file was originally licensed under the terms of the CeCILL v2.1,
  12  * and continues to be available under such terms.
  13  * For more information, see the COPYING file which you should have received
  14  * along with this program.
  15 *
  16 */
  17 /*--------------------------------------------------------------------------*/
  18
  19 #include "signal_gw.hxx"
  20 #include "double.hxx"
  21 #include "string.hxx"
  22 #include "internal.hxx"
  23 #include "function.hxx"
  24 #include "signalprocessingfunctions.hxx"
  25
  26 extern "C"
  27 {
  28 #include "localization.h"
  29 #include "Scierror.h"
  30 #include "sciprint.h"
  31
  32 }
  33
  34 /*--------------------------------------------------------------------------*/
  35 types::Function::ReturnValue sci_corr(types::typed_list &in, int _iRetCount, types::typed_list &out)
  36 {
  37     //check input parameters
  38     if (in.size() < 2 || in.size() > 5)
  39     {
  40         Scierror(77, _("%s: Wrong number of input argument(s): %d to %d expected.\n"), "corr", 2, 5);
  41         return types::Function::Error;
  42     }
  43
  44     //call format
  45     if (in[0]->isString())
  46     {
  47         types::String* pS = in[0]->getAs<types::String>();
  48         if (pS->getSize() == 1 && pS->get(0)[0] == L'f')
  49         {
  50             //[cov,mean]=corr('fft',xmacro,[ymacro],n,sect)
  51
  52             int iErr        = 0;
  53             int iSect       = 0;
  54             int iOutSize    = 0;
  55             int iTotalSize  = 0;
  56             int iSize       = 0;
  57             int iMode       = 0;
  58
  59             double* xa = NULL;
  60             double* xi = NULL;
  61             double* xr = NULL;
  62             double* zr = NULL;
  63             double* zi = NULL;
  64
  65             char *dx = NULL;
  66             char *dy = NULL;
  67             bool bOK = false;
  68
  69             //check input parameters
  70             if (in.size() < 4 || in.size() > 5)
  71             {
  72                 Scierror(77, _("%s: Wrong number of input argument(s): %d to %d expected.\n"), "corr", 4, 5);
  73                 return types::Function::Error;
  74             }
  75
  76             //get parameter sect
  77             int iPos = (int)(in.size() - 1);
  78             if (in[iPos]->isDouble() == false || in[iPos]->getAs<types::Double>()->isScalar() == false)
  79             {
  80                 Scierror(999, _("%s: Wrong type for input argument #%d: A scalar expected.\n"), "corr", iPos + 1);
  81                 return types::Function::Error;
  82             }
  83
  84             iOutSize = (int)in[iPos]->getAs<types::Double>()->get(0);
  85             iSect = iOutSize * 2;
  86
  87             //get parameter n
  88             iPos--;
  89             if (in[iPos]->isDouble() == false || in[iPos]->getAs<types::Double>()->isScalar() == false)
  90             {
  91                 Scierror(999, _("%s: Wrong type for input argument #%d: A scalar expected.\n"), "corr", iPos + 1);
  92                 return types::Function::Error;
  93             }
  94
  95             iTotalSize = (int)in[iPos]->getAs<types::Double>()->get(0);
  96
  97             Signalprocessingfunctions* spFunctionsManager = new Signalprocessingfunctions(L"corr");
  98             Signalprocessing::addSignalprocessingfunctions(spFunctionsManager);
  99
 100             //get xmacro
 101             if (in[1]->isCallable())
 102             {
 103                 spFunctionsManager->setDgetx(in[1]->getAs<types::Callable>());
 104             }
 105             else if (in[1]->isString())
 106             {
 107                 spFunctionsManager->setDgetx(in[1]->getAs<types::String>());
 108             }
 109             else
 110             {
 111                 Scierror(999, _("%s: Wrong type for input argument #%d: A scalar expected.\n"), "corr", iPos + 1);
 112                 return types::Function::Error;
 113             }
 114
 115             iMode = 2;
 116
 117             if (in.size() == 5)
 118             {
 119                 //get ymacro
 120                 if (in[2]->isCallable())
 121                 {
 122                     spFunctionsManager->setDgety(in[2]->getAs<types::Callable>());
 123                 }
 124                 else if (in[2]->isString())
 125                 {
 126                     spFunctionsManager->setDgety(in[2]->getAs<types::String>());
 127                 }
 128                 else
 129                 {
 130                     Scierror(999, _("%s: Wrong type for input argument #%d: A scalar expected.\n"), "corr", iPos + 2);
 131                     return types::Function::Error;
 132                 }
 133
 134                 iMode = 3;
 135             }
 136
 137             xa = new double[iSect];
 138             xr = new double[iSect];
 139             xi = new double[iSect];
 140             zr = new double[iSect / 2 + 1];
 141             zi = new double[iSect / 2 + 1];
 142             C2F(cmpse2)(&iSect, &iTotalSize, &iMode, (void*) dgetx_f, (void*) dgety_f, xa, xr, xi, zr, zi, &iErr);
 143
 144             delete[] xi;
 145             delete[] zr;
 146             delete[] zi;
 147
 148             if (iErr > 0)
 149             {
 150                 delete[] xa;
 151                 delete[] xr;
 152                 Scierror(999, _("fft call : needs power of two!"));
 153                 return types::Function::Error;
 154             }
 155
 156             types::Double *pDblOut1 = new types::Double(1, iOutSize);
 157             pDblOut1->set(xa);
 158             delete[] xa;
 159             out.push_back(pDblOut1);
 160
 161             types::Double *pDblOut2 = new types::Double(1, iMode - 1);
 162             pDblOut2->set(xr);
 163             delete[] xr;
 164             out.push_back(pDblOut2);
 165
 166             return types::Function::OK;
 167         }
 168         else if (pS->getSize() == 1 && pS->get(0)[0] == L'u')
 169         {
 170             types::Double* pDblIn1 = NULL;
 171             types::Double* pDblIn2 = NULL;
 172             types::Double* pDblIn3 = NULL;
 173             types::Double* pDblIn4 = NULL;
 174
 175             int iErr = 0;
 176             int mnx = 0;
 177             int mny = 0;
 178             int mfft = 0;
 179             int nbx = 0;
 180             int ichaud = 0;
 181             int iMode = 0;
 182
 183             double* x = NULL;
 184             double* xu = NULL;
 185             double* xui = NULL;
 186             double* w = NULL;
 187             double* wi = NULL;
 188             double* y = NULL;
 189             double* yi = NULL;
 190
 191             if (in[1]->isDouble() == false)
 192             {
 193                 Scierror(999, _("%s: Wrong type for input argument #%d: Matrix expected.\n"), "corr" , 2);
 194                 return types::Function::Error;
 195             }
 196
 197             pDblIn1 = in[1]->getAs<types::Double>();
 198             if (pDblIn1->isComplex())
 199             {
 200                 Scierror(999, _("%s: Wrong type for input argument #%d: Real matrix expected.\n"), "corr" , 2);
 201                 return types::Function::Error;
 202             }
 203
 204             mnx = pDblIn1->getRows() * pDblIn1->getCols();
 205
 206             x = pDblIn1->get();
 207
 208
 209             if (in[2]->isDouble() == false)
 210             {
 211                 Scierror(999, _("%s: Wrong type for input argument #%d: Matrix expected.\n"), "corr" , 3);
 212                 return types::Function::Error;
 213             }
 214
 215             pDblIn2 = in[2]->getAs<types::Double>();
 216
 217             mny = pDblIn2->getRows() * pDblIn2->getCols();
 218
 219             if (mnx == mny)
 220             {
 221                 iMode = 1;
 222                 if (pDblIn2->isComplex())
 223                 {
 224                     Scierror(999, _("%s: Wrong type for input argument #%d: Real matrix expected.\n"), "corr" , 3);
 225                     return types::Function::Error;
 226                 }
 227
 228                 y = pDblIn2->get();
 229             }
 230
 231             if (iMode == 0)
 232             {
 233                 mfft = mny;
 234                 if (pDblIn2->isComplex() == false)
 235                 {
 236                     double* wtempo = NULL;
 237                     w = new double[pDblIn2->getSize()];
 238                     wi = new double[mfft];
 239                     memset(wi, 0x00, sizeof(double) * mfft);
 240
 241                     wtempo = pDblIn2->get();
 242                     memcpy(w, wtempo, sizeof(double) * pDblIn2->getSize());
 243                 }
 244                 else
 245                 {
 246                     double* wtempo = NULL;
 247                     double* witempo = NULL;
 248                     w = new double[pDblIn2->getSize()];
 249                     wi = new double[pDblIn2->getSize()];
 250
 251                     wtempo = pDblIn2->getReal();
 252                     witempo = pDblIn2->getImg();
 253
 254                     memcpy(w, wtempo, sizeof(double) * pDblIn2->getSize());
 255                     memcpy(wi, witempo, sizeof(double) * pDblIn2->getSize());
 256                 }
 257
 258
 259                 if (in.size() == 4)
 260                 {
 261                     pDblIn3 = in[3]->getAs<types::Double>();
 262                     if (pDblIn3->isComplex())
 263                     {
 264                         Scierror(999, _("%s: Wrong type for input argument #%d: Real matrix expected.\n"), "corr" , 4);
 265                         delete[] wi;
 266                         delete[] w;
 267                         return types::Function::Error;
 268                     }
 269
 270                     xui = new double[mfft * 2];
 271                     double* xutempo = NULL;
 272                     xutempo = pDblIn3->get();
 273                     xu = new double[mfft * 2];
 274                     memset(xu, 0x00, sizeof(double) * mfft * 2);
 275                     memcpy(xu, xutempo, sizeof(double) * pDblIn3->getSize());
 276
 277                     nbx =  pDblIn3->getSize();
 278                     ichaud = 1;
 279                 }
 280                 else
 281                 {
 282                     xu = new double[mfft * 2];
 283                     xui = new double[mfft * 2];
 284                 }
 285
 286                 yi = new double[mny];
 287                 C2F(cmpse3)(&mfft, &mnx, &iMode, x, yi, xu, xui, w, wi, &iErr, &ichaud, &nbx);
 288                 if (iErr > 0)
 289                 {
 290                     delete[] xu;
 291                     delete[] xui;
 292                     delete[] wi;
 293                     delete[] w;
 294                     Scierror(999, _("fft call : needs power of two!"));
 295                     return types::Function::Error;
 296                 }
 297
 298             }
 299             else
 300             {
 301                 pDblIn3 = in[3]->getAs<types::Double>();
 302                 mfft  =   pDblIn3->getRows() * pDblIn3->getCols();
 303                 if (pDblIn3->isComplex() == false)
 304                 {
 305                     wi = new double[mfft];
 306                     memset(wi, 0x00, sizeof(double) * mfft);
 307
 308                     w = new double[pDblIn3->getSize()];
 309                     double* wtempo = NULL;
 310                     wtempo = pDblIn3->get();
 311                     memcpy(w, wtempo, sizeof(double) * pDblIn3->getSize());
 312
 313                 }
 314                 else
 315                 {
 316                     double* wtempo = NULL;
 317                     double* witempo = NULL;
 318                     w = new double[pDblIn3->getSize()];
 319                     wi = new double[pDblIn3->getSize()];
 320
 321                     wtempo = pDblIn3->getReal();
 322                     witempo = pDblIn3->getImg();
 323
 324                     memcpy(w, wtempo, sizeof(double) * pDblIn3->getSize());
 325                     memcpy(wi, witempo, sizeof(double) * pDblIn3->getSize());
 326                 }
 327                 if (in.size() == 5)
 328                 {
 329                     pDblIn4 = in[4]->getAs<types::Double>();
 330                     nbx = pDblIn4->getSize();
 331                     double* xutempo = NULL;
 332                     xutempo = pDblIn4->get();
 333                     xu = new double[mfft * 2];
 334                     memset(xu, 0x00, sizeof(double) * mfft * 2);
 335                     memcpy(xu, xutempo, sizeof(double) * pDblIn4->getSize());
 336                     ichaud = 1;
 337
 338                     xui = new double[mfft * 2];
 339                 }
 340                 else
 341                 {
 342                     xu = new double[mfft * 2];
 343                     xui = new double[mfft * 2];
 344                 }
 345
 346                 C2F(cmpse3)(&mfft, &mnx, &iMode, x, y, xu, xui, w, wi, &iErr, &ichaud, &nbx);
 347                 if (iErr > 0)
 348                 {
 349                     delete[] xu;
 350                     delete[] xui;
 351                     delete[] wi;
 352                     delete[] w;
 353                     Scierror(999, _("fft call : needs power of two!"));
 354                     return types::Function::Error;
 355                 }
 356
 357             }
 358
 359             types::Double *pDblOut1 = NULL;
 360             pDblOut1 = new types::Double(1, mfft, true);
 361             pDblOut1->set(w);
 362             pDblOut1->setImg(wi);
 363             out.push_back(pDblOut1);
 364
 365             if (_iRetCount == 2)
 366             {
 367                 types::Double *pDblOut2 = NULL;
 368                 pDblOut2 = new types::Double(1, mfft / 2);
 369
 370                 for (int i = 0; i < mfft / 2; i++)
 371                 {
 372                     xui[i] = x[mnx - mfft / 2 + i];
 373                 }
 374
 375                 pDblOut2->set(xui);
 376                 out.push_back(pDblOut2);
 377
 378             }
 379             delete[] w;
 380             delete[] wi;
 381             delete[] xui;
 382             delete[] xu;
 383             return types::Function::OK;
 384
 385         }
 386         else
 387         {
 388             //error
 389             Scierror(999, _("%s: Wrong value for input argument #%d: Must be in the set {%s}.\n"), "corr", 1, "'fft', 'update'");
 390             return types::Function::Error;
 391         }
 392     }
 393     else
 394     {
 395         //usual case [cov,mean]=corr(x,[y],nlags)
 396         int iErr                        = 0;
 397         int iCorrelation                = 0;
 398         types::Double* pDblX            = NULL;
 399         types::Double* pDblY            = NULL;
 400         types::Double* pDblCorrelation  = NULL;
 401         types::Double* pDblMean         = NULL;
 402         int iSize                       = 0;
 403         double pdblMean[2];
 404
 405         //check input parameters
 406         if (in.size() < 2 || in.size() > 3)
 407         {
 408             Scierror(77, _("%s: Wrong number of input argument(s): %d to %d expected.\n"), "corr", 2, 3);
 409             return types::Function::Error;
 410         }
 411
 412         //get last parameter nlags
 413         int iPos = (int)(in.size() - 1);
 414         if (in[iPos]->isDouble() == false || in[iPos]->getAs<types::Double>()->isScalar() == false)
 415         {
 416             Scierror(999, _("%s: Wrong type for input argument #%d: A scalar expected.\n"), "corr", iPos + 1);
 417             return types::Function::Error;
 418         }
 419
 420         iCorrelation = (int)in[iPos]->getAs<types::Double>()->get(0);
 421
 422         if (in.size() == 3)
 423         {
 424             if (in[1]->isDouble() == false)
 425             {
 426                 Scierror(999, _("%s: Wrong type for input argument #%d: Matrix expected.\n"), "corr" , 2);
 427                 return types::Function::Error;
 428             }
 429
 430             pDblY = in[1]->getAs<types::Double>();
 431
 432             if (in[0]->isDouble() == false)
 433             {
 434                 Scierror(999, _("%s: Wrong type for input argument #%d: Matrix expected.\n"), "corr" , 1);
 435                 return types::Function::Error;
 436             }
 437
 438             pDblX = in[0]->getAs<types::Double>();
 439
 440             if (pDblX->getSize() != pDblY->getSize())
 441             {
 442                 Scierror(60, _("%s: Wrong size for argument: Incompatible dimensions.\n"), "corr");
 443                 return types::Function::Error;
 444             }
 445         }
 446         else
 447         {
 448             if (in[0]->isDouble() == false)
 449             {
 450                 Scierror(999, _("%s: Wrong type for input argument #%d: Matrix expected.\n"), "corr" , 1);
 451                 return types::Function::Error;
 452             }
 453
 454             pDblX = in[0]->getAs<types::Double>();
 455             pDblY = pDblX;
 456         }
 457
 458         iSize = pDblX->getSize();
 459         pDblCorrelation = new types::Double(1, iCorrelation);
 460         C2F(tscccf)(pDblX->get(), pDblY->get(), &iSize, pDblCorrelation->get(), pdblMean, &iCorrelation, &iErr);
 461         if (iErr == -1)
 462         {
 463             delete pDblCorrelation;
 464             Scierror(999, _("%s: Too many coefficients are required.\n"), "corr");
 465             return types::Function::Error;
 466         }
 467
 468         out.push_back(pDblCorrelation);
 469
 470         if (_iRetCount == 2)
 471         {
 472             if (in.size() == 3)
 473             {
 474                 pDblMean = new types::Double(1, 2);
 475             }
 476             else
 477             {
 478                 pDblMean = new types::Double(1, 1);
 479             }
 480
 481             pDblMean->set(pdblMean);
 482             out.push_back(pDblMean);
 483         }
 484     }
 485
 486     return types::Function::OK;
 487 }
 488