* Bug 15907 fixed: filter was corrupting its input state array
[scilab.git] / scilab / modules / signal_processing / sci_gateway / cpp / sci_filter.cpp
1 /*
2  * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
3  * Copyright (C) 2017 - ESI - Antoine ELIAS
4  *
5  * Copyright (C) 2012 - 2016 - Scilab Enterprises
6  *
7  * This file is hereby licensed under the terms of the GNU GPL v2.0,
8  * pursuant to article 5.3.4 of the CeCILL v.2.1.
9  * This file was originally licensed under the terms of the CeCILL v2.1,
10  * and continues to be available under such terms.
11  * For more information, see the COPYING file which you should have received
12  * along with this program.
13  *
14  */
15
16 /*
17
18 Copyright (C) 1996-2017 John W. Eaton
19
20 This file is part of Octave.
21
22 Octave is free software; you can redistribute it and/or modify it
23 under the terms of the GNU General Public License as published by
24 the Free Software Foundation; either version 3 of the License, or
25 (at your option) any later version.
26
27 Octave is distributed in the hope that it will be useful, but
28 WITHOUT ANY WARRANTY; without even the implied warranty of
29 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
30 GNU General Public License for more details.
31
32 You should have received a copy of the GNU General Public License
33 along with Octave; see the file COPYING.  If not, see
34 <http://www.gnu.org/licenses/>.
35
36 */
37
38 // Based on Tony Richardson's filter.m.
39 //
40 // Originally translated to C++ by KH (Kurt.Hornik@wu-wien.ac.at)
41 // with help from Fritz Leisch and Andreas Weingessel on Oct 20, 1994.
42 //
43 // Rewritten to use templates to handle both real and complex cases by
44 // jwe, Wed Nov  1 19:15:29 1995.
45
46 #include "signal_gw.hxx"
47 #include "double.hxx"
48 #include "polynom.hxx"
49 #include "function.hxx"
50 #include "overload.hxx"
51
52 extern "C"
53 {
54 #include "sci_malloc.h"
55 #include "localization.h"
56 #include "Scierror.h"
57 }
58
59 const char fname[] = "filter";
60
61 static void clean_filter(types::Double* a, bool alloc_a, types::Double* b, bool alloc_b)
62 {
63     if (alloc_a)
64     {
65         a->killMe();
66     }
67
68     if (alloc_b)
69     {
70         b->killMe();
71     }
72 }
73
74 static types::Double* filter(types::Double* b, types::Double* a, types::Double* x, types::Double* si)
75 {
76     bool alloc_a = false;
77     bool alloc_b = false;
78
79     int size_a = a->getSize();
80     int size_b = b->getSize();
81
82     int len = std::max(size_a, size_b);
83
84     int dims[2] = {1, len};
85
86     //uniformize
87     types::Double* new_a = a->resize(dims, 2)->getAs<types::Double>();
88     if (new_a != a)
89     {
90         alloc_a = true;
91     }
92
93     types::Double* new_b = b->resize(dims, 2)->getAs<types::Double>();
94     if (new_b != b)
95     {
96         alloc_b = true;
97     }
98
99     double* pa = new_a->get();
100     double* pb = new_b->get();
101     double* px = x->get();
102
103     double norm = pa[0];
104
105     if (norm == 0)
106     {
107         clean_filter(new_a, alloc_a, new_b, alloc_b);
108         Scierror(999, _("%s: Wrong value for input argument #%d: First element must not be %s.\n"), fname, 2, "0");
109         return nullptr;
110     }
111
112     int size_x = x->getSize();
113
114     int rows_si = si->getRows();
115
116     if (rows_si != len - 1)
117     {
118         clean_filter(new_a, alloc_a, new_b, alloc_b);
119         Scierror(999, "%s: first dimension of SI must be of length max (length (a), length (b)) - 1 (%d) (%d)", fname, rows_si, size_x);
120         return nullptr;
121     }
122
123     if (size_x == 0)
124     {
125         clean_filter(new_a, alloc_a, new_b, alloc_b);
126         return x;
127     }
128
129     if (norm != 1)
130     {
131         if (alloc_a == false)
132         {
133             new_a = new_a->clone();
134         }
135
136         if (alloc_b == false)
137         {
138             new_b = new_b->clone();
139         }
140
141         pa = new_a->get();
142         pb = new_b->get();
143         alloc_a = true;
144         alloc_b = true;
145
146         for (int i = 0; i < len; ++i)
147         {
148             pa[i] /= norm;
149             pb[i] /= norm;
150         }
151     }
152
153     if (size_a <= 1 && rows_si <= 0)
154     {
155         types::Double* ret = x->clone();
156         double* pret = ret->get();
157         int s = ret->getSize();
158         for (int i = 0; i < s; ++i)
159         {
160             pret[i] *= pb[0];
161         }
162
163         clean_filter(new_a, alloc_a, new_b, alloc_b);
164         return ret;
165     }
166
167     types::Double* y = new types::Double(x->getRows(), x->getCols());
168
169     double* py = y->get();
170     double* psi = si->get();
171     if (size_a > 1)
172     {
173         if (rows_si > 0)
174         {
175             for (int i = 0; i < size_x; ++i)
176             {
177                 py[i] = psi[0] + pb[0] * px[i];
178                 for (int j = 0; j < rows_si - 1; ++j)
179                 {
180                     psi[j] = psi[j + 1] - pa[j + 1] * py[i] + pb[j + 1] * px[i];
181                 }
182
183                 psi[rows_si - 1] = pb[rows_si] * px[i] - pa[rows_si] * py[i];
184             }
185         }
186         else
187         {
188             for (int i = 0; i < size_x; ++i)
189             {
190                 py[i] = psi[0] + pb[0] * px[i];
191                 psi[0] = pb[rows_si] * px[i] - pa[rows_si] * py[i];
192             }
193
194         }
195     }
196     else
197     {
198         for (int i = 0; i < size_x; ++i)
199         {
200             py[i] = psi[0] + pb[0] * px[i];
201
202             if (rows_si > 1)
203             {
204                 for (int j = 0; j < rows_si - 1; j++)
205                 {
206                     psi[j] = psi[j + 1] + pb[j + 1] * px[i];
207                 }
208
209                 psi[rows_si - 1] = pb[rows_si] * px[i];
210             }
211             else
212             {
213                 psi[0] = pb[1] * px[i];
214             }
215         }
216     }
217
218
219     clean_filter(new_a, alloc_a, new_b, alloc_b);
220     return y;
221 }
222 /*--------------------------------------------------------------------------*/
223 types::Function::ReturnValue sci_filter(types::typed_list &in, int _iRetCount, types::typed_list &out)
224 {
225     int iRhs = (int)in.size();
226
227     types::Double* b = nullptr;
228     types::Double* a = nullptr;
229     types::Double* x = nullptr;
230     types::Double* si = nullptr;
231
232     if (iRhs < 3 || iRhs > 4)
233     {
234         Scierror(999, _("%s: Wrong number of input arguments: %d to %d expected.\n"), fname, 3, 4);
235         return types::Function::Error;
236     }
237
238     types::InternalType::ScilabId type_b = in[0]->getId();
239     types::InternalType::ScilabId type_a = in[1]->getId();
240     types::InternalType::ScilabId type_x = in[2]->getId();
241     types::InternalType::ScilabId type_si = types::InternalType::IdDouble;
242
243     if (iRhs > 3)
244     {
245         type_si = in[3]->getId();
246     }
247
248     if (type_b  != types::InternalType::IdDouble && type_b  != types::InternalType::IdScalarDouble ||
249         type_a  != types::InternalType::IdDouble && type_a  != types::InternalType::IdScalarDouble ||
250         type_x  != types::InternalType::IdDouble && type_x  != types::InternalType::IdScalarDouble ||
251         type_si != types::InternalType::IdDouble && type_si != types::InternalType::IdScalarDouble)
252     {
253         return Overload::call(L"%_filter", in, _iRetCount, out);
254     }
255
256     b = in[0]->getAs<types::Double>();
257     a = in[1]->getAs<types::Double>();
258     x = in[2]->getAs<types::Double>();
259
260     if (iRhs > 3)
261     {
262         si = in[3]->getAs<types::Double>()->clone();
263     }
264     else
265     {
266         int size_a = a->getSize();
267         int size_b = b->getSize();
268         int len = std::max(size_a, size_b) - 1;
269         int leadDim = x->getRows() == 1 ? 1 : 0;
270
271         int si_dims[2] = {x->getRows(), x->getCols()};
272         for (int i = leadDim; i > 0; i--)
273         {
274             si_dims[i] = si_dims[i - 1];
275         }
276
277         si_dims[0] = len;
278
279         si = new types::Double(2, si_dims);
280         si->setZeros();
281     }
282
283     if (b->isVector() == false)
284     {
285         Scierror(999, _("%s: Wrong size for input argument #%d: Vector expected.\n"), fname, 1);
286         return types::Function::Error;
287     }
288
289     if (a->isVector() == false)
290     {
291         Scierror(999, _("%s: Wrong size for input argument #%d: Vector expected.\n"), fname, 2);
292         return types::Function::Error;
293     }
294
295     if (x->isVector() == false)
296     {
297         Scierror(999, _("%s: Wrong size for input argument #%d: Vector expected.\n"), fname, 3);
298         return types::Function::Error;
299     }
300
301     if (iRhs > 3 && si->isVector() == false)
302     {
303         Scierror(999, _("%s: Wrong size for input argument #%d: Vector expected.\n"), fname, 4);
304         return types::Function::Error;
305     }
306
307     types::Double* ret = filter(b, a, x, si);
308     if (_iRetCount != 2)
309     {
310         si->killMe();
311     }
312
313     if (ret == nullptr)
314     {
315         return types::Function::Error;
316     }
317
318     out.push_back(ret);
319     if (_iRetCount == 2)
320     {
321         out.push_back(si);
322     }
323     return types::Function::OK;
324 }