fe765bb5d2608d89d3e347e04ed30f5e221ce507
[scilab.git] / 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