72a0f9e5cabbdedbb2049ec36f7b5a6f943423ff
[scilab.git] / scilab / modules / signal_processing / sci_gateway / cpp / sci_syredi.cpp
1 /*
2 * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
3 * Copyright (C) 2011 - DIGITEO - Antoine ELIAS
4 *
5 * This file must be used under the terms of the CeCILL.
6 * This source file is licensed as described in the file COPYING, which
7 * you should have received as part of this distribution.  The terms
8 * are also available at
9 * http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
10 *
11 */
12 /*--------------------------------------------------------------------------*/
13
14 #include "signal_gw.hxx"
15 #include "double.hxx"
16 #include "function.hxx"
17
18 extern "C"
19 {
20 #include <math.h>
21 #include "sci_malloc.h"
22 #include "localization.h"
23 #include "Scierror.h"
24 #include "core_math.h"
25
26     extern void C2F(syredi)(int *maxdeg, int *ityp, int *iapro,
27                             double *om, double *adelp, double *adels,
28                             /* outputs */
29                             int *deg_count, int *zeros_count,
30                             double *fact,
31                             double *b2, double *b1, double *b0,
32                             double *c1, double *c0,
33                             double *zzr, double *zzi,
34                             double *zpr, double *zpi,
35                             int *ierr,
36                             /* working buffers */
37                             double *spr, double *spi,
38                             double *pren, double *pimn,
39                             double *zm, double *sm, double *rom,
40                             /* v-- doutful types but whatever... */
41                             double *nzero, double *nze);
42 }
43
44 enum filter_type
45 {
46     low_pass = 1,
47     high_pass = 2,
48     band_pass = 3,
49     stop_band = 4
50 };
51
52 enum design_type
53 {
54     butterworth = 1,
55     elliptic = 2,
56     chebytchev1 = 3,
57     chebytchev2 = 4
58 };
59
60 #define OUT_COUNT 18
61
62 //local functions
63 double maximum(double* _pDblVal, int _iSize);
64 double minimum(double* _pDblVal, int _iSize);
65 bool isSortedAscending(double* _pdblVal, int _iSize);
66 void reshapeFilters(types::Double* _pDblInR, types::Double* _pDblInI, types::Double* _pDblOut);
67
68 /*--------------------------------------------------------------------------*/
69 types::Function::ReturnValue sci_syredi(types::typed_list &in, int _iRetCount, types::typed_list &out)
70 {
71     int iMaxDeg                         = 64;
72     types::Double *pDblType             = NULL;
73     filter_type iType;
74
75     types::Double *pDblAppro            = NULL;
76     design_type iAppro;
77
78     types::Double *pDblCutOff           = NULL;
79
80     types::Double *pDblDeltaP           = NULL;
81     double dblDeltaP                    = 0;
82
83     types::Double *pDblDeltaS           = NULL;
84     double dblDeltaS                    = 0;
85
86     types::Double* pDblOut[OUT_COUNT];
87     int piOutSize[OUT_COUNT]            = {32, 32, 32, 32, 32, 64, 64, 64, 64, 64, 64, 64, 64, 256, 256, 16, 64, 64};
88
89     int iErr                            = 0;
90     int iZeroCount                      = 0;
91     int iDegCount                       = 0;
92
93     types::Double* pDblFact             = NULL;
94     double dblFact                      = 0;
95
96     types::Double* pOut                 = NULL;
97
98     //check input parameters
99     if (in.size() != 5)
100     {
101         Scierror(77, _("%s: Wrong number of input argument(s): %d expected.\n"), "syredi", 5);
102         return types::Function::Error;
103     }
104
105     //check 1st input parameter : filter type ( 1 int )
106     pDblType = in[0]->getAs<types::Double>();
107     if (in[0]->isDouble() == false || pDblType->isScalar() == false || pDblType->isComplex() == true)
108     {
109         Scierror(999, _("%s: Wrong type for argument %d: Real scalar expected.\n"), "syredi", 1);
110         return types::Function::Error;
111     }
112
113     iType = (filter_type)(int)in[0]->getAs<types::Double>()->get(0);
114
115     //check 2nd input parameter : approximation type ( 1 int )
116     pDblAppro = in[1]->getAs<types::Double>();
117     if (in[1]->isDouble() == false || pDblAppro->isScalar() == false || pDblAppro->isComplex() == true)
118     {
119         Scierror(999, _("%s: Wrong type for input argument #%d: Real scalar expected.\n"), "syredi", 2);
120         return types::Function::Error;
121     }
122
123     iAppro = (design_type)(int)in[1]->getAs<types::Double>()->get(0);
124
125     //check 3rd input parameter : cuttof frequencies ( 4-row vector )
126     pDblCutOff = in[2]->getAs<types::Double>();
127     if (in[2]->isDouble() == false || pDblCutOff->getSize() != 4 || pDblCutOff->getCols() != 4)
128     {
129         Scierror(999, _("%s: Wrong size for input argument #%d: A %d-by-%d array expected.\n"), "syredi", 3, 1, 4);
130         return types::Function::Error;
131     }
132
133     if (minimum(pDblCutOff->get(), pDblCutOff->getSize()) < 0 || minimum(pDblCutOff->get(), pDblCutOff->getSize()) > M_PI)
134     {
135         Scierror(999, _("%s: Wrong value for input argument #%d: Must be in the interval [%s, %s].\n"), "syredi", 3, "0", "%pi");
136         return types::Function::Error;
137     }
138
139     if ((iType == low_pass || iType == high_pass) && isSortedAscending(pDblCutOff->get(), 2) == false)
140     {
141         Scierror(999, _("%s: Wrong values for input argument #%d: Elements must be in increasing order.\n"), "syredi", 3);
142         return types::Function::Error;
143     }
144
145     if ((iType == band_pass || iType == stop_band) && isSortedAscending(pDblCutOff->get(), 4) == false)
146     {
147         Scierror(999, _("%s: Wrong values for input argument #%d: Elements must be in increasing order.\n"), "syredi", 3);
148         return types::Function::Error;
149     }
150
151     //check 4th input parameter : ripple in passband ( 0 < deltap < 1 )
152     pDblDeltaP = in[3]->getAs<types::Double>();
153     if (in[3]->isDouble() == false || pDblAppro->isScalar() == false || pDblAppro->isComplex() == true)
154     {
155         Scierror(999, _("%s: Wrong type for input argument #%d: Real scalar expected.\n"), "syredi", 4);
156         return types::Function::Error;
157     }
158
159     dblDeltaP = pDblDeltaP->get(0);
160
161     //check 5th input parameter : ripple in stopband ( 0 < deltas < 1 )
162     pDblDeltaS = in[4]->getAs<types::Double>();
163     if (in[4]->isDouble() == false || pDblDeltaS->isScalar() == false || pDblDeltaS->isComplex() == true)
164     {
165         Scierror(999, _("%s: Wrong type for input argument #%d: Real scalar expected.\n"), "syredi", 5);
166         return types::Function::Error;
167     }
168
169     dblDeltaS = pDblDeltaS->get(0);
170
171
172     //alloc temporary variables
173     for (int i = 0 ; i < OUT_COUNT ; i++)
174     {
175         pDblOut[i] = new types::Double(1, piOutSize[i]);
176     }
177
178     //call math function
179     C2F(syredi)(&iMaxDeg, (int*)&iType, (int*)&iAppro, pDblCutOff->get(), &dblDeltaP,
180                 &dblDeltaS, &iZeroCount, &iDegCount, &dblFact,
181                 pDblOut[0]->get(), pDblOut[1]->get(), pDblOut[2]->get(), pDblOut[3]->get(),
182                 pDblOut[4]->get(), pDblOut[5]->get(), pDblOut[6]->get(), pDblOut[7]->get(),
183                 pDblOut[8]->get(), &iErr, pDblOut[9]->get(), pDblOut[10]->get(), pDblOut[11]->get(),
184                 pDblOut[12]->get(), pDblOut[13]->get(), pDblOut[14]->get(), pDblOut[15]->get(),
185                 pDblOut[16]->get(), pDblOut[17]->get());
186
187     if (iErr)
188     {
189         if (iErr == -7)
190         {
191             Scierror(999, _("%s: specs => invalid order filter.\n"), "syredi");
192             return types::Function::Error;
193         }
194         else if (iErr == -9)
195         {
196             Scierror(999, _("%s: specs => too high order filter.\n"), "syredi");
197             return types::Function::Error;
198         }
199         else
200         {
201             Scierror(999, _("%s: error in function syredi.\n"), "syredi");
202             return types::Function::Error;
203         }
204     }
205
206     //1st output : fact
207     pDblFact = new types::Double(dblFact);
208     out.push_back(pDblFact);
209
210     //2nd output : b2
211     pOut = new types::Double(1, iDegCount);
212     pOut->set(pDblOut[0]->get());
213     out.push_back(pOut);
214
215     //3rd output : b1
216     pOut = new types::Double(1, iDegCount);
217     pOut->set(pDblOut[1]->get());
218     out.push_back(pOut);
219
220     //4th output : b0
221     pOut = new types::Double(1, iDegCount);
222     pOut->set(pDblOut[2]->get());
223     out.push_back(pOut);
224
225     //5th output : c1
226     pOut = new types::Double(1, iDegCount);
227     pOut->set(pDblOut[3]->get());
228     out.push_back(pOut);
229
230     //6th output : c0
231     pOut = new types::Double(1, iDegCount);
232     pOut->set(pDblOut[4]->get());
233     out.push_back(pOut);
234
235     //7th output : zeros
236     pOut = new types::Double(1, iZeroCount, true);
237     reshapeFilters(pDblOut[5], pDblOut[6], pOut);
238     out.push_back(pOut);
239
240     //8th output : poles
241     pOut = new types::Double(1, iZeroCount, true);
242     reshapeFilters(pDblOut[7], pDblOut[8], pOut);
243     out.push_back(pOut);
244
245     //clear temporary variables
246     for (int i = 0 ; i < OUT_COUNT ; i++)
247     {
248         delete pDblOut[i];
249     }
250
251     return types::Function::OK;
252 }
253
254 /* returns the maximum of the values array (with length elements) */
255 double maximum(double* _pDblVal, int _iSize)
256 {
257     double dblMax = 0;
258     if (_iSize < 1)
259     {
260         return dblMax;
261     }
262
263     dblMax = _pDblVal[0];
264     for (int i = 1 ; i < _iSize ; i++)
265     {
266         dblMax = Max(_pDblVal[i], dblMax);
267     }
268     return dblMax;
269 }
270
271 /* returns the minimum of the values array (with length elements) */
272 double minimum(double* _pDblVal, int _iSize)
273 {
274     double dblMin = 0;
275     if (_iSize < 1)
276     {
277         return dblMin;
278     }
279
280     dblMin = _pDblVal[0];
281     for (int i = 1 ; i < _iSize ; i++)
282     {
283         dblMin = Min(_pDblVal[i], dblMin);
284     }
285     return dblMin;
286 }
287
288 bool isSortedAscending(double* _pdblVal, int _iSize)
289 {
290     for (int i = 1 ; i < _iSize ; i++)
291     {
292         if (_pdblVal[i - 1] > _pdblVal[i])
293         {
294             return false;
295         }
296     }
297     return true;
298 }
299
300 int syredi_buffered(/* inputs */
301     enum filter_type filter, enum design_type design,
302     double cutoff_frequencies[4], double ripple_passband, double ripple_stopband,
303     /* outputs */
304     int *zeros_count, int *deg_count, double *fact,
305     double b2[], double b1[], double b0[],
306     double c1[], double c0[],
307     double zzr[], double zzi[],
308     double zpr[], double zpi[])
309 {
310 #define BUFFERS_COUNT 9
311     const int buffers_count = BUFFERS_COUNT;
312     int buffers_sizes[BUFFERS_COUNT] = {64, 64, 64, 64, 256, 256, 16, 64, 64};
313     double *buffers[BUFFERS_COUNT];
314 #undef BUFFERS_COUNT
315     int maxdeg = 64, error = 1, i;
316
317     for (i = 0 ; i < buffers_count ; ++i)
318     {
319         buffers[i] = (double *)MALLOC(buffers_sizes[i] * sizeof(double));
320         if (buffers[i] == NULL)
321         {
322             break;
323         }
324     }
325     if (i == buffers_count)
326     {
327         C2F(syredi)(/* inputs */
328             &maxdeg, (int *) &filter, (int *) &design,
329             cutoff_frequencies, &ripple_passband, &ripple_stopband,
330             /* outputs */
331             zeros_count, deg_count, fact,
332             b2, b1, b0, c1, c0,
333             zzr, zzi, zpr, zpi,
334             &error,
335             /* buffers */
336             buffers[0], buffers[1], buffers[2], buffers[3],
337             buffers[4], buffers[5],
338             buffers[6],
339             buffers[7], buffers[8]);
340     }
341     for (i-- ; i >= 0 ; i--)
342     {
343         FREE(buffers[i]);
344         buffers[i] = NULL;
345     }
346     return error;
347 }
348
349 void reshapeFilters(types::Double* _pDblInR, types::Double* _pDblInI, types::Double* _pDblOut)
350 {
351     int iSize           = _pDblOut->getSize();
352     double* pdblInR     = _pDblInR->getReal();
353     double* pdblInI     = _pDblInI->getReal();
354     double* pdblOutR    = _pDblOut->getReal();
355     double* pdblOutI    = _pDblOut->getImg();
356
357     for (int i = 0, j = 0 ; j < iSize ; i++, j++)
358     {
359         if (pdblInI[i] == 0)
360         {
361             pdblOutR[j] = pdblInR[i];
362             pdblOutI[j] = 0;
363         }
364         else
365         {
366             pdblOutR[j] = pdblInR[i];
367             pdblOutI[j] = pdblInI[i];
368             j++;
369             pdblOutR[j] = pdblInR[i];
370             pdblOutI[j] = - pdblInI[i];
371         }
372     }
373 }