lib file was not embedded in the binary version
[scilab.git] / scilab / modules / signal_processing / sci_gateway / c / sci_syredi.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
14 #include "api_scilab.h"
15 #include "gw_signal.h"
16 #include "MALLOC.h"
17 #include "stack-c.h"
18 #include "Scierror.h"
19 #include "localization.h"
20 #include <math.h>
21
22 extern void C2F(syredi)(int *maxdeg, int *ityp, int *iapro,
23                         double *om, double *adelp, double *adels,
24                         /* outputs */
25                         int *deg_count, int *zeros_count,
26                         double *fact,
27                         double *b2, double *b1, double *b0,
28                         double *c1, double *c0,
29                         double *zzr, double *zzi,
30                         double *zpr, double *zpi,
31                         int *ierr,
32                         /* working buffers */
33                         double *spr, double *spi,
34                         double *pren, double *pimn,
35                         double *zm, double *sm, double *rom,
36                         /* v-- doutful types but whatever... */
37                         double *nzero, double *nze);
38
39 enum filter_type {
40   low_pass = 1,
41   high_pass = 2,
42   band_pass = 3,
43   stop_band = 4};
44
45 enum design_type {
46   butterworth = 1,
47   elliptic = 2,
48   chebytchev1 = 3,
49   chebytchev2 = 4};
50
51 int syredi_buffered(enum filter_type filter, enum design_type design,
52                     double cutoff_frequencies[4],
53                     double ripple_passband, double ripple_stopband,
54                     /* outputs */
55                     int *deg_count, int *zeros_count, double *fact,
56                     double b2[], double b1[], double b0[],
57                     double c1[], double c0[],
58                     double zzr[], double zzi[],
59                     double zpr[], double zpi[]);
60
61 void questionable_copy(int count,
62                        double input_real[], double input_imaginary[],
63                        double output_real[], double output_imaginary[]);
64
65 int coherent_cutoffs(double cutoffs[], enum filter_type kind);
66
67 int is_sorted_ascending(double values[], int length);
68
69 double maximum(double values[], int count);
70 double minimum(double values[], int count);
71
72 int sci_syredi(char *fname, int* _piKey)
73 {
74   enum filter_type filter;
75   enum design_type design;
76   int rows, cols, one = 1;
77   double *argument = NULL, *cutoff_frequencies = NULL;
78   double ripple_passband, ripple_stopband;
79
80   /* output variables */
81   int deg_count, zeros_count, error;
82   double *fact, *b0, *b1, *b2, *c0, *c1;
83   double
84     *zeros_real = NULL, *zeros_imaginary = NULL,
85     *poles_real = NULL, *poles_imaginary = NULL;
86
87   /* output buffers used for computation, ultimately copied to output variables */
88   #define OUTPUT_BUFFERS_COUNT 9
89   int output_buffers_count = OUTPUT_BUFFERS_COUNT;
90   int output_buffers_sizes[OUTPUT_BUFFERS_COUNT] =
91     {32, 32, 32, 32, 32, 64, 64, 64, 64};
92   double *output_buffers[OUTPUT_BUFFERS_COUNT];
93   int output_cursor;
94   #undef OUTPUT_BUFFERS_COUNT
95
96   int *p;
97
98   CheckRhs(5,5);
99   CheckLhs(8,8);
100
101   /* arg1: filter kind */
102   // GetRhsVarMatrixDouble(1, &rows, &cols, &argument);
103   getVarAddressFromPosition(_piKey, 1, &p);
104   getMatrixOfDouble(_piKey, p, &rows, &cols, &argument);
105   filter = (int)argument[0];
106
107   /* arg2: approximation kind */
108   // GetRhsVarMatrixDouble(2, &rows, &cols, &argument);
109   getVarAddressFromPosition(_piKey, 2, &p);
110   getMatrixOfDouble(_piKey, p, &rows, &cols, &argument);
111   design = (int)argument[0];
112
113   /* arg3: cutoff frequencies */
114   // GetRhsVarMatrixDouble(3, &rows, &cols, &cutoff_frequencies);
115   getVarAddressFromPosition(_piKey, 3, &p);
116   getMatrixOfDouble(_piKey, p, &rows, &cols, &cutoff_frequencies);
117   if (rows != 1 || cols != 4) {
118     Scierror(999, _("%s: Wrong size for input argument #%d: A %d-by-%d array expected.\n"), fname, 3, 1, 4);
119     return 1;
120   }
121   if (minimum(cutoff_frequencies, 4) < 0 || maximum(cutoff_frequencies, 4) > M_PI) {
122     Scierror(999, _("%s: Wrong value for input argument #%d: Must be in the interval [%s, %s].\n"), fname, 4, "0", "%pi");
123     return 1;
124   }
125   if (((filter == low_pass || filter == high_pass)
126        && !is_sorted_ascending(cutoff_frequencies, 2)) ||
127       ((filter == band_pass || filter == stop_band)
128        && !is_sorted_ascending(cutoff_frequencies, 4))) {
129     Scierror(999, _("%s: Wrong values for input argument #%d: Elements must be in increasing order.\n"), fname, 3);
130     return 1;
131   }
132
133   /* arg4: passband's ripple */
134   // GetRhsVarMatrixDouble(4, &rows, &cols, &argument);
135   getVarAddressFromPosition(_piKey, 4, &p);
136   getMatrixOfDouble(_piKey, p, &rows, &cols, &argument);
137   ripple_passband = argument[0];
138   if (ripple_passband <= 0 || ripple_passband >= 1) {
139     Scierror(999, _("%s: Wrong value for input argument #%d: Must be in the interval [%s, %s].\n"), fname, 4, "0", "1");
140     return 1;
141   }
142
143   /* arg5: stopband's ripple */
144   // GetRhsVarMatrixDouble(5, &rows, &cols, &argument);
145   getVarAddressFromPosition(_piKey, 5, &p);
146   getMatrixOfDouble(_piKey, p, &rows, &cols, &argument);
147   ripple_stopband = argument[0];
148   if (ripple_stopband <= 0 || ripple_stopband >= 1) {
149     Scierror(999, _("%s: Wrong value for input argument #%d: Must be in the interval [%s, %s].\n"), fname, 5, "0", "1");
150     return 1;
151   }
152
153   for (output_cursor = 0 ; output_cursor < output_buffers_count ; ++output_cursor) {
154     output_buffers[output_cursor] = (double *)MALLOC(output_buffers_sizes[output_cursor] * sizeof(double));
155     if (output_buffers[output_cursor] == NULL) {
156       Scierror(999, _("%s : Memory allocation error.\n"), fname);
157       // TODO: free buffers (elegantly, possibly use an abstraction for this)
158       return 1;
159     }
160   }
161
162   // iAllocMatrixOfDouble(Rhs + 1, 1, 1, &fact);
163   allocMatrixOfDouble(_piKey, Rhs + 1, 1, 1, &fact);
164   createMatrixOfDouble(_piKey, Rhs + 1, 1, 1, fact);
165
166   error = syredi_buffered(/* inputs */
167                           filter, design,
168                           cutoff_frequencies, ripple_passband, ripple_stopband,
169                           /* outputs */
170                           &zeros_count, &deg_count, fact,
171                           output_buffers[0], output_buffers[1],
172                           output_buffers[2], output_buffers[3],
173                           output_buffers[4], output_buffers[5],
174                           output_buffers[6], output_buffers[7],
175                           output_buffers[8]);
176   if (error) {
177     Scierror(999, _("%s : Memory allocation error.\n"), fname);
178     return 1;
179   }
180
181   /* ret1: fact */
182   LhsVar(1) = Rhs + 1;
183
184   /* ret2: b2 */
185   // iAllocMatrixOfDouble(Rhs + 2, 1, deg_count, &b2);
186   allocMatrixOfDouble(_piKey, Rhs + 2, 1, deg_count, &b2);
187   //createMatrixOfDouble(Rhs + 2, 1, deg_count, b2);
188   C2F(dcopy)(&deg_count, output_buffers[0], &one, b2, &one);
189   LhsVar(2) = Rhs + 2;
190
191   /* ret3: b1 */
192   // iAllocMatrixOfDouble(Rhs + 3, 1, deg_count, &b1);
193   allocMatrixOfDouble(_piKey, Rhs + 3, 1, deg_count, &b1);
194   //createMatrixOfDouble(Rhs + 3, 1, deg_count, b1);
195   C2F(dcopy)(&deg_count, output_buffers[1], &one, b1, &one);
196   LhsVar(3) = Rhs + 3;
197
198   /* ret4: b0 */
199   // iAllocMatrixOfDouble(Rhs + 4, 1, deg_count, &b0);
200   allocMatrixOfDouble(_piKey, Rhs + 4, 1, deg_count, &b0);
201   //createMatrixOfDouble(Rhs + 4, 1, deg_count, b0);
202   C2F(dcopy)(&deg_count, output_buffers[2], &one, b0, &one);
203   LhsVar(4) = Rhs + 4;
204
205   /* ret5: c1 */
206   // iAllocMatrixOfDouble(Rhs + 5, 1, deg_count, &c1);
207   allocMatrixOfDouble(_piKey, Rhs + 5, 1, deg_count, &c1);
208   //createMatrixOfDouble(Rhs + 5, 1, deg_count, c1);
209   C2F(dcopy)(&deg_count, output_buffers[3], &one, c1, &one);
210   LhsVar(5) = Rhs + 5;
211
212   /* ret6: c0 */
213   // iAllocMatrixOfDouble(Rhs + 6, 1, deg_count, &c0);
214   allocMatrixOfDouble(_piKey, Rhs + 6, 1, deg_count, &c0);
215   //createMatrixOfDouble(Rhs + 6, 1, deg_count, c0);
216   C2F(dcopy)(&deg_count, output_buffers[4], &one, c0, &one);
217   LhsVar(6) = Rhs + 6;
218
219   /* ret7: zeros */
220   // iAllocComplexMatrixOfDouble(Rhs + 7, 1, zeros_count, &zeros_real, &zeros_imaginary);
221   allocComplexMatrixOfDouble(_piKey, Rhs + 7, 1, zeros_count, &zeros_real, &zeros_imaginary);
222   //createComplexMatrixOfDouble(Rhs + 7, 1, zeros_count, zeros_real, zeros_imaginary);
223   questionable_copy(zeros_count, output_buffers[5], output_buffers[6],
224                     zeros_real, zeros_imaginary);
225   LhsVar(7) = Rhs + 7;
226
227   /* ret8: poles */
228   // iAllocComplexMatrixOfDouble(Rhs + 8, 1, zeros_count, &poles_real, &poles_imaginary);
229   allocComplexMatrixOfDouble(_piKey, Rhs + 8, 1, zeros_count, &poles_real, &poles_imaginary);
230   //createComplexMatrixOfDouble(Rhs + 8, 1, zeros_count, poles_real, poles_imaginary);
231   questionable_copy(zeros_count, output_buffers[7], output_buffers[8], poles_real, poles_imaginary);
232   LhsVar(8) = Rhs + 8;
233
234   for (output_cursor-- ; output_cursor >= 0 ; output_cursor--) {
235     FREE(output_buffers[output_cursor]);
236     output_buffers[output_cursor] = NULL;
237   }
238
239   PutLhsVar();
240   return error;
241 }
242
243 int syredi_buffered(/* inputs */
244                     enum filter_type filter, enum design_type design,
245                     double cutoff_frequencies[4], double ripple_passband, double ripple_stopband,
246                     /* outputs */
247                     int *zeros_count, int *deg_count, double *fact,
248                     double b2[], double b1[], double b0[],
249                     double c1[], double c0[],
250                     double zzr[], double zzi[],
251                     double zpr[], double zpi[])
252 {
253   #define BUFFERS_COUNT 9
254   const int buffers_count = BUFFERS_COUNT;
255   int buffers_sizes[BUFFERS_COUNT] = {64, 64, 64, 64, 256, 256, 16, 64, 64};
256   double *buffers[BUFFERS_COUNT];
257   #undef BUFFERS_COUNT
258   int maxdeg = 64, error = 1, i;
259
260   for (i = 0 ; i < buffers_count ; ++i) {
261     buffers[i] = (double *)MALLOC(buffers_sizes[i] * sizeof(double));
262     if (buffers[i] == NULL) {
263       break;
264     }
265   }
266   if (i == buffers_count) {
267     C2F(syredi)(/* inputs */
268                 &maxdeg, (int *) &filter, (int *) &design,
269                 cutoff_frequencies, &ripple_passband, &ripple_stopband,
270                 /* outputs */
271                 zeros_count, deg_count, fact,
272                 b2, b1, b0, c1, c0,
273                 zzr, zzi, zpr, zpi,
274                 &error,
275                 /* buffers */
276                 buffers[0], buffers[1], buffers[2], buffers[3],
277                 buffers[4], buffers[5],
278                 buffers[6],
279                 buffers[7], buffers[8]);
280   }
281   for (i-- ; i >= 0 ; i--) {
282     FREE(buffers[i]);
283     buffers[i] = NULL;
284   }
285   return error;
286 }
287
288 /* copies input to output adding conjugates */
289 /* questionable until the reason of input data loss is clarified */
290 /* translated from sci_syredi.f */
291 void questionable_copy(int length,
292                        double input_real[], double input_imaginary[],
293                        double output_real[], double output_imaginary[])
294 {
295   int i, j;
296   for (i = 0, j = 0 ; j < length ; ++i, ++j) {
297     output_real[j] = input_real[i];
298     output_imaginary[j] = input_imaginary[i];
299     /* doubles compared with == as in fortran's source */
300     if (input_imaginary[i] != 0L) {
301       ++j;
302       output_real[j] = input_real[i];
303       output_imaginary[j] = - input_imaginary[i];
304     }
305   }
306 }
307
308 /* checks whether the values array (with length elements) is ascendingly sorted */
309 int is_sorted_ascending(double values[], int length) {
310   int i;
311   for (i = 0 ; i < length - 1 ; ++i) {
312     if (values[i] > values[i + 1]) {
313       return 0;
314     }
315   }
316   return 1;
317 }
318
319 /* returns the maximum of the values array (with length elements) */
320 double maximum(double values[], int length) {
321   int i;
322   double max;
323   for (i = 0, max = values[0] ; i < length ; ++i) {
324     max = Max(values[i], max);
325   }
326   return max;
327 }
328
329 /* returns the minimum of the values array (with length elements) */
330 double minimum(double values[], int length) {
331   int i;
332   double min;
333   for (i = 0, min = values[0] ; i < length ; ++i) {
334     min = Min(values[i], min);
335   }
336   return min;
337 }