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                         return types::Function::Error;
 266                     }
 267
 268                     xui = new double[mfft * 2];
 269                     double* xutempo = NULL;
 270                     xutempo = pDblIn3->get();
 271                     xu = new double[mfft * 2];
 272                     memset(xu, 0x00, sizeof(double) * mfft * 2);
 273                     memcpy(xu, xutempo, sizeof(double) * pDblIn3->getSize());
 274
 275                     nbx =  pDblIn3->getSize();
 276                     ichaud = 1;
 277                 }
 278                 else
 279                 {
 280                     xu = new double[mfft * 2];
 281                     xui = new double[mfft * 2];
 282                 }
 283
 284                 yi = new double[mny];
 285                 C2F(cmpse3)(&mfft, &mnx, &iMode, x, yi, xu, xui, w, wi, &iErr, &ichaud, &nbx);
 286                 if (iErr > 0)
 287                 {
 288                     delete[] xu;
 289                     delete[] xui;
 290                     delete[] wi;
 291                     delete[] w;
 292                     Scierror(999, _("fft call : needs power of two!"));
 293                     return types::Function::Error;
 294                 }
 295
 296             }
 297             else
 298             {
 299                 pDblIn3 = in[3]->getAs<types::Double>();
 300                 mfft  =   pDblIn3->getRows() * pDblIn3->getCols();
 301                 if (pDblIn3->isComplex() == false)
 302                 {
 303                     wi = new double[mfft];
 304                     memset(wi, 0x00, sizeof(double) * mfft);
 305
 306                     w = new double[pDblIn3->getSize()];
 307                     double* wtempo = NULL;
 308                     wtempo = pDblIn3->get();
 309                     memcpy(w, wtempo, sizeof(double) * pDblIn3->getSize());
 310
 311                 }
 312                 else
 313                 {
 314                     double* wtempo = NULL;
 315                     double* witempo = NULL;
 316                     w = new double[pDblIn3->getSize()];
 317                     wi = new double[pDblIn3->getSize()];
 318
 319                     wtempo = pDblIn3->getReal();
 320                     witempo = pDblIn3->getImg();
 321
 322                     memcpy(w, wtempo, sizeof(double) * pDblIn3->getSize());
 323                     memcpy(wi, witempo, sizeof(double) * pDblIn3->getSize());
 324                 }
 325                 if (in.size() == 5)
 326                 {
 327                     pDblIn4 = in[4]->getAs<types::Double>();
 328                     nbx = pDblIn4->getSize();
 329                     double* xutempo = NULL;
 330                     xutempo = pDblIn4->get();
 331                     xu = new double[mfft * 2];
 332                     memset(xu, 0x00, sizeof(double) * mfft * 2);
 333                     memcpy(xu, xutempo, sizeof(double) * pDblIn4->getSize());
 334                     ichaud = 1;
 335
 336                     xui = new double[mfft * 2];
 337                 }
 338                 else
 339                 {
 340                     xu = new double[mfft * 2];
 341                     xui = new double[mfft * 2];
 342                 }
 343
 344                 C2F(cmpse3)(&mfft, &mnx, &iMode, x, y, xu, xui, w, wi, &iErr, &ichaud, &nbx);
 345                 if (iErr > 0)
 346                 {
 347                     delete[] xu;
 348                     delete[] xui;
 349                     delete[] wi;
 350                     delete[] w;
 351                     Scierror(999, _("fft call : needs power of two!"));
 352                     return types::Function::Error;
 353                 }
 354
 355             }
 356
 357             types::Double *pDblOut1 = NULL;
 358             pDblOut1 = new types::Double(1, mfft, true);
 359             pDblOut1->set(w);
 360             pDblOut1->setImg(wi);
 361             out.push_back(pDblOut1);
 362
 363             if (_iRetCount == 2)
 364             {
 365                 types::Double *pDblOut2 = NULL;
 366                 pDblOut2 = new types::Double(1, mfft / 2);
 367
 368                 for (int i = 0; i < mfft / 2; i++)
 369                 {
 370                     xui[i] = x[mnx - mfft / 2 + i];
 371                 }
 372
 373                 pDblOut2->set(xui);
 374                 out.push_back(pDblOut2);
 375
 376             }
 377             delete[] w;
 378             delete[] wi;
 379             delete[] xui;
 380             delete[] xu;
 381             return types::Function::OK;
 382
 383         }
 384         else
 385         {
 386             //error
 387             Scierror(999, _("%s: Wrong value for input argument #%d: Must be in the set {%s}.\n"), "corr", 1, "'fft', 'update'");
 388             return types::Function::Error;
 389         }
 390     }
 391     else
 392     {
 393         //usual case [cov,mean]=corr(x,[y],nlags)
 394         int iErr                        = 0;
 395         int iCorrelation                = 0;
 396         types::Double* pDblX            = NULL;
 397         types::Double* pDblY            = NULL;
 398         types::Double* pDblCorrelation  = NULL;
 399         types::Double* pDblMean         = NULL;
 400         int iSize                       = 0;
 401         double pdblMean[2];
 402
 403         //check input parameters
 404         if (in.size() < 2 || in.size() > 3)
 405         {
 406             Scierror(77, _("%s: Wrong number of input argument(s): %d to %d expected.\n"), "corr", 2, 3);
 407             return types::Function::Error;
 408         }
 409
 410         //get last parameter nlags
 411         int iPos = (int)(in.size() - 1);
 412         if (in[iPos]->isDouble() == false || in[iPos]->getAs<types::Double>()->isScalar() == false)
 413         {
 414             Scierror(999, _("%s: Wrong type for input argument #%d: A scalar expected.\n"), "corr", iPos + 1);
 415             return types::Function::Error;
 416         }
 417
 418         iCorrelation = (int)in[iPos]->getAs<types::Double>()->get(0);
 419
 420         if (in.size() == 3)
 421         {
 422             if (in[1]->isDouble() == false)
 423             {
 424                 Scierror(999, _("%s: Wrong type for input argument #%d: Matrix expected.\n"), "corr" , 2);
 425                 return types::Function::Error;
 426             }
 427
 428             pDblY = in[1]->getAs<types::Double>();
 429
 430             if (in[0]->isDouble() == false)
 431             {
 432                 Scierror(999, _("%s: Wrong type for input argument #%d: Matrix expected.\n"), "corr" , 1);
 433                 return types::Function::Error;
 434             }
 435
 436             pDblX = in[0]->getAs<types::Double>();
 437
 438             if (pDblX->getSize() != pDblY->getSize())
 439             {
 440                 Scierror(60, _("%s: Wrong size for argument: Incompatible dimensions.\n"), "corr");
 441                 return types::Function::Error;
 442             }
 443         }
 444         else
 445         {
 446             if (in[0]->isDouble() == false)
 447             {
 448                 Scierror(999, _("%s: Wrong type for input argument #%d: Matrix expected.\n"), "corr" , 1);
 449                 return types::Function::Error;
 450             }
 451
 452             pDblX = in[0]->getAs<types::Double>();
 453             pDblY = pDblX;
 454         }
 455
 456         iSize = pDblX->getSize();
 457         pDblCorrelation = new types::Double(1, iCorrelation);
 458         C2F(tscccf)(pDblX->get(), pDblY->get(), &iSize, pDblCorrelation->get(), pdblMean, &iCorrelation, &iErr);
 459         if (iErr == -1)
 460         {
 461             delete pDblCorrelation;
 462             Scierror(999, _("%s: Too many coefficients are required.\n"), "corr");
 463             return types::Function::Error;
 464         }
 465
 466         out.push_back(pDblCorrelation);
 467
 468         if (_iRetCount == 2)
 469         {
 470             if (in.size() == 3)
 471             {
 472                 pDblMean = new types::Double(1, 2);
 473             }
 474             else
 475             {
 476                 pDblMean = new types::Double(1, 1);
 477             }
 478
 479             pDblMean->set(pdblMean);
 480             out.push_back(pDblMean);
 481         }
 482     }
 483
 484     return types::Function::OK;
 485 }
 486