944c6545359c948e3804e0134733ebbf6503707b
[scilab.git] / scilab / modules / signal_processing / sci_gateway / c / sci_remez.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
20 extern void C2F(remez)(int *ngrid, int *nc, int *iext,
21                        double *ad, double *x, double *y,
22                        float *des, float *grid, float *wt,
23                        double *a, double *p, double *q, double *alpha);
24
25 /* allocates required buffers and calls remez */
26 int remez_buffered(int ngrid, int nc, int *iext,
27                    float *des, float *grid, float *wt, double *output);
28
29 int sci_remez(char *fname, int* _piKey)
30 {
31   /************************************************
32    * Warning : bug 4189                           *
33    * The remez function returns                   *
34    * an array which last element is meaningless.  *
35    * -> sementic of the fortran gw preserved      *
36    * -> watch the curious nc's cooking            *
37    ************************************************/
38
39   int rows, cols, length, ngrid = 0, nc = 0, error = 0;
40   double *output = NULL, *argument = NULL;
41   float *des = NULL, *grid = NULL, *wt = NULL;
42   int *iext;
43   int *p;
44
45   CheckRhs(4,4);
46   CheckLhs(1,1);
47
48   // GetRhsVarMatrixDouble(1, &rows, &cols, &argument);
49   getVarAddressFromPosition(_piKey, 1, &p);
50   getMatrixOfDouble(_piKey, p, &rows, &cols, &argument);
51   iext = (int *)argument;
52   nc = cols * rows;
53   C2F(entier)(&nc, argument, iext);
54
55   // GetRhsVarMatrixDouble(2, &rows, &cols, &argument);
56   getVarAddressFromPosition(_piKey, 2, &p);
57   getMatrixOfDouble(_piKey, p, &rows, &cols, &argument);
58   des = (float *)argument;
59   ngrid = cols * rows;
60   length = rows;
61   C2F(simple)(&ngrid, argument, des);
62
63   // GetRhsVarMatrixDouble(3, &rows, &cols, &argument);
64   getVarAddressFromPosition(_piKey, 3, &p);
65   getMatrixOfDouble(_piKey, p, &rows, &cols, &argument);
66   grid = (float *)argument;
67   C2F(simple)(&ngrid, argument, grid);
68
69   // GetRhsVarMatrixDouble(4, &rows, &cols, &argument);
70   getVarAddressFromPosition(_piKey, 4, &p);
71   getMatrixOfDouble(_piKey, p, &rows, &cols, &argument);
72   wt = (float *)argument;
73   C2F(simple)(&ngrid, argument, wt);
74
75   // iAllocMatrixOfDouble(Rhs + 1, rows, nc - 1, &output);
76   allocMatrixOfDouble(_piKey, Rhs + 1, rows, nc - 1, &output);
77   //createMatrixOfDouble(_piKey, Rhs + 1, rows, nc - 1, output);
78
79   error = remez_buffered(ngrid, nc - 2, iext, des, grid, wt, output);
80   if (error) {
81     Scierror(999, _("%s : Memory allocation error.\n"), fname);
82     return 1;
83   }
84
85   LhsVar(1) = Rhs + 1;
86   PutLhsVar();
87
88   return 0;
89 }
90
91 int remez_buffered(int ngrid, int nc, int *iext,
92                    float *des, float *grid, float *wt, double *output) {
93   int one = 1;
94   double *buffer0 = NULL,
95     *buffer1 = NULL, *buffer2 = NULL,
96     *buffer3 = NULL, *buffer4 = NULL,
97     *buffer5 = NULL, *buffer6 = NULL;
98
99   buffer0 = (double *)MALLOC((nc + 2) * sizeof(double));
100   buffer1 = (double *)MALLOC((nc + 2) * sizeof(double));
101   buffer2 = (double *)MALLOC((nc + 2) * sizeof(double));
102   buffer3 = (double *)MALLOC((nc + 2) * sizeof(double));
103   buffer4 = (double *)MALLOC((nc + 2) * sizeof(double));
104   buffer5 = (double *)MALLOC((nc + 2) * sizeof(double));
105   buffer6 = (double *)MALLOC((nc + 2) * sizeof(double));
106   if (buffer0 == NULL ||
107       buffer1 == NULL || buffer2 == NULL ||
108       buffer3 == NULL || buffer4 == NULL ||
109       buffer5 == NULL || buffer6 == NULL) {
110     return 1;
111   }
112
113   C2F(remez)(&ngrid, &nc, iext, buffer1, buffer2, buffer3,
114              des, grid, wt,
115              buffer4, buffer5, buffer6, buffer0);
116   nc++;
117   C2F(dcopy)(&nc, buffer0, &one, output, &one);
118
119   FREE(buffer0);
120   FREE(buffer6);
121   FREE(buffer5);
122   FREE(buffer4);
123   FREE(buffer3);
124   FREE(buffer2);
125   FREE(buffer1);
126
127   return 0;
128 }