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