66148a70cbb04cab3125d37522837565bfdebf3a
[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  * Copyright (C) 2012 - 2016 - Scilab Enterprises
7  *
8  * This file is hereby licensed under the terms of the GNU GPL v2.0,
9  * pursuant to article 5.3.4 of the CeCILL v.2.1.
10  * This file was originally licensed under the terms of the CeCILL v2.1,
11  * and continues to be available under such terms.
12  * For more information, see the COPYING file which you should have received
13  * along with this program.
14  *
15  */
16
17 #include "gw_signal.h"
18 #include "api_scilab.h"
19 #include "sci_malloc.h"
20 #include "Scierror.h"
21 #include "localization.h"
22 #include "elem_common.h"
23
24 extern void C2F(remez) (int *ngrid, int *nc, int *iext,
25                         double *ad, double *x, double *y, float *des, float *grid, float *wt, double *a, double *p, double *q, double *alpha);
26
27 /* allocates required buffers and calls remez */
28 int remez_buffered(int ngrid, int nc, int *iext, float *des, float *grid, float *wt, double *output);
29
30 int sci_remez(char *fname, void *pvApiCtx)
31 {
32     /************************************************
33      * Warning : bug 4189                           *
34      * The remez function returns                   *
35      * an array which last element is meaningless.  *
36      * -> sementic of the fortran gw preserved      *
37      * -> watch the curious nc's cooking            *
38      ************************************************/
39
40     int rows, cols, length, ngrid = 0, nc = 0, error = 0;
41     double *output = NULL, *argument = NULL;
42     float *des = NULL, *grid = NULL, *wt = NULL;
43     int *iext;
44     int *p;
45
46     CheckRhs(4, 4);
47     CheckLhs(1, 1);
48
49     // GetRhsVarMatrixDouble(1, &rows, &cols, &argument);
50     getVarAddressFromPosition(pvApiCtx, 1, &p);
51     getMatrixOfDouble(pvApiCtx, p, &rows, &cols, &argument);
52     iext = (int *)argument;
53     nc = cols * rows;
54     C2F(entier) (&nc, argument, iext);
55
56     // GetRhsVarMatrixDouble(2, &rows, &cols, &argument);
57     getVarAddressFromPosition(pvApiCtx, 2, &p);
58     getMatrixOfDouble(pvApiCtx, p, &rows, &cols, &argument);
59     des = (float *)argument;
60     ngrid = cols * rows;
61     length = rows;
62     C2F(simple) (&ngrid, argument, des);
63
64     // GetRhsVarMatrixDouble(3, &rows, &cols, &argument);
65     getVarAddressFromPosition(pvApiCtx, 3, &p);
66     getMatrixOfDouble(pvApiCtx, p, &rows, &cols, &argument);
67     grid = (float *)argument;
68     C2F(simple) (&ngrid, argument, grid);
69
70     // GetRhsVarMatrixDouble(4, &rows, &cols, &argument);
71     getVarAddressFromPosition(pvApiCtx, 4, &p);
72     getMatrixOfDouble(pvApiCtx, p, &rows, &cols, &argument);
73     wt = (float *)argument;
74     C2F(simple) (&ngrid, argument, wt);
75
76     // iAllocMatrixOfDouble(Rhs + 1, rows, nc - 1, &output);
77     allocMatrixOfDouble(pvApiCtx, Rhs + 1, rows, nc - 1, &output);
78     //createMatrixOfDouble(pvApiCtx, Rhs + 1, rows, nc - 1, output);
79
80     error = remez_buffered(ngrid, nc - 2, iext, des, grid, wt, output);
81     if (error)
82     {
83         Scierror(999, _("%s : Memory allocation error.\n"), fname);
84         return 1;
85     }
86
87     LhsVar(1) = Rhs + 1;
88     PutLhsVar();
89
90     return 0;
91 }
92
93 int remez_buffered(int ngrid, int nc, int *iext, float *des, float *grid, float *wt, double *output)
94 {
95     int one = 1;
96     double *buffer0 = NULL, *buffer1 = NULL, *buffer2 = NULL, *buffer3 = NULL, *buffer4 = NULL, *buffer5 = NULL, *buffer6 = NULL;
97
98     buffer0 = (double *)MALLOC((nc + 2) * sizeof(double));
99     buffer1 = (double *)MALLOC((nc + 2) * sizeof(double));
100     buffer2 = (double *)MALLOC((nc + 2) * sizeof(double));
101     buffer3 = (double *)MALLOC((nc + 2) * sizeof(double));
102     buffer4 = (double *)MALLOC((nc + 2) * sizeof(double));
103     buffer5 = (double *)MALLOC((nc + 2) * sizeof(double));
104     buffer6 = (double *)MALLOC((nc + 2) * sizeof(double));
105     if (buffer0 == NULL || buffer1 == NULL || buffer2 == NULL || buffer3 == NULL || buffer4 == NULL || buffer5 == NULL || buffer6 == NULL)
106     {
107         return 1;
108     }
109
110     C2F(remez) (&ngrid, &nc, iext, buffer1, buffer2, buffer3, des, grid, wt, buffer4, buffer5, buffer6, buffer0);
111     nc++;
112     C2F(dcopy) (&nc, buffer0, &one, output, &one);
113
114     FREE(buffer0);
115     FREE(buffer6);
116     FREE(buffer5);
117     FREE(buffer4);
118     FREE(buffer3);
119     FREE(buffer2);
120     FREE(buffer1);
121
122     return 0;
123 }