880afac6bac92e9c7a0f32ca747c02d9a72bf088
[scilab.git] / scilab / modules / signal_processing / sci_gateway / c / sci_corr.c
1 /*
2  * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
3  * Copyright (C) 2006 - INRIA - Allan CORNET
4  * Copyright (C) 2009 - Digiteo - Vincent LIARD
5  *
6  * This file must be used under the terms of the CeCILL.
7  * This source file is licensed as described in the file COPYING, which
8  * you should have received as part of this distribution.  The terms
9  * are also available at
10  * http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
11  */
12
13 #include "api_scilab.h"
14 #include "MALLOC.h"
15 #include "gw_signal.h"
16 #include "stack-c.h"
17 #include "Scierror.h"
18 #include "localization.h"
19 #include "../../../statistics/src/c/sci_string_matrix.h"
20
21 extern int C2F(scicorr)(char *id,unsigned long fname_len );
22
23 int corr_default(char *fname, int* _piKey);
24 int corr_updates(char *fname, int* _piKey);
25 int corr_fft(char *fname, int* _piKey);
26
27 /* dispatch to specialized gateways according to mode (default, 'updt', 'fft') */
28 int sci_corr(char *fname, int* _piKey)
29 {
30   char *name;
31   int type;
32   int *arg;
33
34   CheckRhs(2,5);
35   CheckLhs(1,2);
36   getVarAddressFromPosition(_piKey, 1, &arg);
37   getVarType(_piKey, arg, &type);
38         if (type == sci_strings) {
39     name = create_string(_piKey, 1);
40     switch (name[0]) {
41     case 'f': /* fft */
42       CheckRhs(4,5);
43       corr_fft(fname, _piKey);
44       break;
45     case 'u': /* updates */
46       CheckRhs(3,5);
47       corr_updates(fname, _piKey);
48       break;
49     default:
50       Scierror(999, _("%s: Wrong value for input argument #%d: Must be in the set {%s}.\n"), fname, 1, "'fft', 'update'");
51       return 1;
52     }
53     destroy_string(name);
54   }
55   else { /* default */
56     CheckRhs(2,3);
57     corr_default(fname, _piKey);
58   }
59   return 0;
60 }
61
62 int corr_default(char *fname, int* _piKey) {
63
64   extern void C2F(tscccf)(double *x, double *y, int *length,
65                           double *cxy, double *xymean, int *lag, int *error);
66
67   int length, lags_number, error, rows1, cols1, rows2, cols2, rows3, cols3;
68   double *x, *y, *xvar /* crossvariance */, *mean, *data;
69   int *arg;
70   int type;
71
72   getVarAddressFromPosition(_piKey, 1, &arg);
73   getVarDimension(_piKey, arg, &rows1, &cols1);
74   length = rows1 * cols1;
75   getVarType(_piKey, arg, &type);
76   if ((type != sci_matrix) ||
77       (rows1 != 1 && cols1 != 1)) {
78     Scierror(999, _("%s: Wrong size for input argument #%d: A vector expected.\n"), fname, 1);
79     return 1;
80   }
81   getMatrixOfDouble(_piKey, arg, &rows1, &cols1, &x);
82
83   /* interpret 2nd argument according to number of arguments */
84   getVarAddressFromPosition(_piKey, 2, &arg);
85   getMatrixOfDouble(_piKey, arg, &rows2, &cols2, &data);
86   if (Rhs == 2) {
87     if (rows2 != 1 || cols2 != 1) {
88       Scierror(999, _("%s: Wrong size for input argument #%d: A scalar expected.\n"), fname, 2);
89       return 1;
90     }
91     y = x;
92     lags_number = (int)data[0];
93   }
94   else {
95     if (rows1 != rows2 || cols1 != cols2) {
96       Scierror(999, _("%s: Incompatible input arguments #%d and #%d': Same sizes expected.\n"), fname, 1, 2);
97       return 1;
98     }
99     y = data;
100     getVarAddressFromPosition(_piKey, 3, &arg);
101     getMatrixOfDouble(_piKey, arg, &rows3, &cols3, &data);
102     lags_number = (int)data[0];
103   }
104   allocMatrixOfDouble(_piKey, Rhs + 1, 1, lags_number, &xvar);
105   allocMatrixOfDouble(_piKey, Rhs + 2, 1, (x == y) ? 1 : 2, &mean);
106   C2F(tscccf)(x, y, &length, xvar, mean, &lags_number, &error);
107   LhsVar(1) = Rhs + 1;
108   LhsVar(2) = Rhs + 2;
109   PutLhsVar();
110   return error;
111 }
112
113 int corr_fft(char *fname, int* _piKey) {
114
115   extern void C2F(cmpse2)(int *m, int *n, int *mode,
116                           void (*bgetx)(double *, int *, int *),
117                           void (*bgety)(double *, int *, int *),
118                           double *xa, double *xr, double *xi,
119                           double *zr, double *zi, int *error);
120   extern void C2F(bgetx)(double *x, int *increment, int *start_idx);
121   extern void C2F(bgety)(double *x, int *increment, int *start_idx);
122   extern void C2F(setdgetx)(char *name, int *rep);
123   extern void C2F(setdgety)(char *name, int *rep);
124
125   enum {self = 2, cross = 3} mode = self;
126   int error, rows, cols, n, mm, lag;
127   double *r;
128   double *xa, *xr, *xi, *zr, *zi, *data;
129   int return2_length, one = 1;
130   int *arg;
131
132   /* taken from DllmainSignal_processing.c */
133   /* the very f*** global referred to underneath */
134   extern struct {
135     int kgxtop, kgytop, ksec, kisc;
136   } C2F(corradr);
137
138   getVarAddressFromPosition(_piKey, Rhs, &arg);
139   getMatrixOfDouble(_piKey, arg, &rows, &cols, &data);
140   lag = (int)(*data);
141   mm = 2 * lag;
142   getVarAddressFromPosition(_piKey, Rhs - 1, &arg);
143   getMatrixOfDouble(_piKey, arg, &rows, &cols, &data);
144   n = (int)(*data);
145   mode = (Rhs == 4 ? self : cross);
146
147   /** @todo add code to manage calls to compiled functions */
148
149   /* setting the f*** (!@&£*% etc) global parameter for cmpse2 */
150   C2F(corradr).kgxtop = Top - Rhs + 2;
151   C2F(corradr).kgytop = Top - Rhs + 3; /* useless for self correlation but doesn't hurt */
152   C2F(corradr).ksec = Top;
153   C2F(corradr).kisc = Top + 1;
154
155   xa = (double *)MALLOC(mm * sizeof(double));
156   xr = (double *)MALLOC(mm * sizeof(double));
157   xi = (double *)MALLOC(mm * sizeof(double));
158   zr = (double *)MALLOC(mm * sizeof(double));
159   zi = (double *)MALLOC(mm * sizeof(double));
160   {
161     int imode = mode; /* enum* not equiv to int* :'( */
162     C2F(cmpse2)(&mm, &n, &imode,
163                 C2F(bgetx), C2F(bgety),
164                 xa, xr, xi, zr, zi, &error);
165   }
166   allocMatrixOfDouble(_piKey, Rhs + 1, 1, lag, &r);
167   C2F(dcopy)(&lag, xa, &one, r, &one);
168
169   LhsVar(1) = Rhs + 1;
170   return2_length = Rhs - 3;
171   allocMatrixOfDouble(_piKey, Rhs + 2, 1, return2_length, &r);
172   C2F(dcopy)(&return2_length, xr, &one, r, &one);
173   FREE(zi);
174   FREE(zr);
175   FREE(xi);
176   FREE(xr);
177   FREE(xa);
178   LhsVar(2) = Rhs + 2;
179   PutLhsVar();
180   return 0;
181 }
182
183 /* assert: RhsVar must have been checked to hold 4 to 5 elements */
184 int corr_updates(char *fname, int* _piKey) {
185
186   extern void C2F(cmpse3)(int *m, int *n, int *mode, double *x, double *y,
187                           double *xr, double *xi,
188                           double *zr, double *zi,
189                           int *error, int *ichaud, int *nbx);
190
191   enum {self = 0, cross = 1} mode = self;
192   int update = 0;
193   int mfft, wpos;
194   int length, error, rows, cols;
195   double *x, *y,
196     *zr, *zi, *wr, *wi, *x_real,
197     *rr, *ri, y_value;
198   int *arg;
199   int one = 1, nbx = 0;
200   double zero = 0L;
201
202   /* x */
203   getVarAddressFromPosition(_piKey, 2, &arg);
204   getMatrixOfDouble(_piKey, arg, &rows, &cols, &x);
205   length = rows * cols;
206   getVarAddressFromPosition(_piKey, 3, &arg);
207   getVarDimension(_piKey, arg, &rows, &cols);
208   if (length == rows * cols) {
209     mode = cross;
210   }
211   if ((mode == cross && Rhs == 5) ||
212       (mode == self && Rhs == 4)) {
213     update = 1;
214   }
215
216   /* y */
217   if (mode == cross) {
218     CheckRhs(4,5);
219     if (isVarComplex(_piKey, arg)) {
220       Scierror(999, _("%s: Input argument #%d must be real.\n"), fname, 3);
221       return 1;
222     }
223     getMatrixOfDouble(_piKey, arg, &rows, &cols, &y);
224   }
225   else {
226     y_value = 0;
227     y = &y_value;
228   }
229
230   /* w */
231   wpos = (mode == self) ? 3 : 4;
232   /* dimension ought to be a power of 2 (could be checked here) */
233   /* GetVarDimension(wpos, &rows, &cols); */
234   getVarAddressFromPosition(_piKey, wpos, &arg);
235   if (isVarComplex(_piKey, arg)) {
236     getComplexMatrixOfDouble(_piKey, arg, &rows, &cols, &zr, &zi);
237     mfft = rows * cols;
238   }
239   else {
240     getMatrixOfDouble(_piKey, arg, &rows, &cols, &zr);
241     mfft = rows * cols;
242     zi = (double *)MALLOC(mfft * sizeof(double));
243     C2F(dset)(&mfft, &zero, zi, &one);
244   }
245
246   /* xu */
247   wr = (double *)MALLOC(mfft * sizeof(double));
248   wi = (double *)MALLOC(mfft * sizeof(double));
249   C2F(dset)(&mfft, &zero, wr, &one);
250   C2F(dset)(&mfft, &zero, wi, &one);
251   if (update) { /* if xu provided, update it */
252     getVarAddressFromPosition(_piKey, wpos + 1, &arg);
253     getMatrixOfDouble(_piKey, arg, &rows, &cols, &x_real);
254     nbx = rows * cols;
255     if (nbx > mfft) {
256       /* TODO: make error msg more accurate (2nd dim >= dim(arg)) */
257       Scierror(999, _("%s: Wrong size for input argument #%d: A %d-by-%d matrix expected.\n"), fname, 5, 1, 999);
258       return 1;
259     }
260     C2F(dcopy)(&nbx, x_real, &one, wr, &one);
261   }
262   else {
263     nbx = 0;
264   }
265
266   C2F(cmpse3)(&mfft, &length, (int *) &mode, x, y,
267               wr, wi, zr, zi,
268               &error, (int *) &update, &nbx);
269
270   if (error) {
271     Scierror(999, _("%s: Wrong size for input argument #%d: Should be a power of 2"), fname, wpos);
272     return error;
273   }
274   /** @warning x built of data on the stack that's reaffected by matrix allocation
275      -> copied to keep it safe for lhs2 output */
276   C2F(dcopy)(&nbx, &(x[length-nbx]), &one, wr, &one);
277   allocComplexMatrixOfDouble(_piKey, Rhs + 1, 1, mfft, &rr, &ri);
278   C2F(dcopy)(&mfft, zr, &one, rr, &one);
279   C2F(dcopy)(&mfft, zi, &one, ri, &one);
280   LhsVar(1) = Rhs + 1;
281   if (Lhs == 2) {
282     allocMatrixOfDouble(_piKey, Rhs + 2, 1, nbx, &rr);
283     C2F(dcopy)(&nbx, wr, &one, rr, &one);
284     LhsVar(2) = Rhs + 2;
285   }
286   PutLhsVar();
287   FREE(wr);
288   FREE(wi);
289   return error;
290 }