signal_processing plugged.
[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 "MALLOC.h"
14 #include "gw_signal.h"
15 #include "stack-c.h"
16 #include "Scierror.h"
17 #include "localization.h"
18
19 #include "api_scilab.h"
20 #include "api_oldstack.h"
21
22 extern void C2F(remez)(int *ngrid, int *nc, int *iext,
23                        double *ad, double *x, double *y,
24                        float *des, float *grid, float *wt,
25                        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,
29                    float *des, float *grid, float *wt, double *output);
30
31 int sci_remez(char *fname, int* _piKey)
32 {
33   /************************************************
34    * Warning : bug 4189                           *
35    * The remez function returns                   *
36    * an array which last element is meaningless.  *
37    * -> sementic of the fortran gw preserved      *
38    * -> watch the curious nc's cooking            *
39    ************************************************/
40
41   int rows, cols, length, ngrid = 0, nc = 0, error = 0;
42   double *output = NULL, *argument = NULL;
43   float *des = NULL, *grid = NULL, *wt = NULL;
44   int *iext;
45   int *p;
46
47   CheckRhs(4,4);
48   CheckLhs(1,1);
49
50   // GetRhsVarMatrixDouble(1, &rows, &cols, &argument);
51   getVarAddressFromPosition(_piKey, 1, &p);
52   getMatrixOfDouble(_piKey, p, &rows, &cols, &argument);
53   iext = (int *)argument;
54   nc = cols * rows;
55   C2F(entier)(&nc, argument, iext);
56
57   // GetRhsVarMatrixDouble(2, &rows, &cols, &argument);
58   getVarAddressFromPosition(_piKey, 2, &p);
59   getMatrixOfDouble(_piKey, p, &rows, &cols, &argument);
60   des = (float *)argument;
61   ngrid = cols * rows;
62   length = rows;
63   C2F(simple)(&ngrid, argument, des);
64
65   // GetRhsVarMatrixDouble(3, &rows, &cols, &argument);
66   getVarAddressFromPosition(_piKey, 3, &p);
67   getMatrixOfDouble(_piKey, p, &rows, &cols, &argument);
68   grid = (float *)argument;
69   C2F(simple)(&ngrid, argument, grid);
70
71   // GetRhsVarMatrixDouble(4, &rows, &cols, &argument);
72   getVarAddressFromPosition(_piKey, 4, &p);
73   getMatrixOfDouble(_piKey, p, &rows, &cols, &argument);
74   wt = (float *)argument;
75   C2F(simple)(&ngrid, argument, wt);
76
77   // iAllocMatrixOfDouble(Rhs + 1, rows, nc - 1, &output);
78   allocMatrixOfDouble(_piKey, Rhs + 1, rows, nc - 1, &output);
79   //createMatrixOfDouble(_piKey, Rhs + 1, rows, nc - 1, output);
80
81   error = remez_buffered(ngrid, nc - 2, iext, des, grid, wt, output);
82   if (error) {
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,
94                    float *des, float *grid, float *wt, double *output) {
95   int one = 1;
96   double *buffer0 = NULL,
97     *buffer1 = NULL, *buffer2 = NULL,
98     *buffer3 = NULL, *buffer4 = NULL,
99     *buffer5 = NULL, *buffer6 = NULL;
100
101   buffer0 = (double *)MALLOC((nc + 2) * sizeof(double));
102   buffer1 = (double *)MALLOC((nc + 2) * sizeof(double));
103   buffer2 = (double *)MALLOC((nc + 2) * sizeof(double));
104   buffer3 = (double *)MALLOC((nc + 2) * sizeof(double));
105   buffer4 = (double *)MALLOC((nc + 2) * sizeof(double));
106   buffer5 = (double *)MALLOC((nc + 2) * sizeof(double));
107   buffer6 = (double *)MALLOC((nc + 2) * sizeof(double));
108   if (buffer0 == NULL ||
109       buffer1 == NULL || buffer2 == NULL ||
110       buffer3 == NULL || buffer4 == NULL ||
111       buffer5 == NULL || buffer6 == NULL) {
112     return 1;
113   }
114
115   C2F(remez)(&ngrid, &nc, iext, buffer1, buffer2, buffer3,
116              des, grid, wt,
117              buffer4, buffer5, buffer6, buffer0);
118   nc++;
119   C2F(dcopy)(&nc, buffer0, &one, output, &one);
120
121   FREE(buffer0);
122   FREE(buffer6);
123   FREE(buffer5);
124   FREE(buffer4);
125   FREE(buffer3);
126   FREE(buffer2);
127   FREE(buffer1);
128
129   return 0;
130 }