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