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     bool alloc_si = false;
 228     types::Double* b = nullptr;
 229     types::Double* a = nullptr;
 230     types::Double* x = nullptr;
 231     types::Double* si = nullptr;
 232
 233     if (iRhs < 3 || iRhs > 4)
 234     {
 235         Scierror(999, _("%s: Wrong number of input arguments: %d to %d expected.\n"), fname, 3, 4);
 236         return types::Function::Error;
 237     }
 238
 239     types::InternalType::ScilabId type_b = in[0]->getId();
 240     types::InternalType::ScilabId type_a = in[1]->getId();
 241     types::InternalType::ScilabId type_x = in[2]->getId();
 242     types::InternalType::ScilabId type_si = types::InternalType::IdDouble;
 243
 244     if (iRhs > 3)
 245     {
 246         type_si = in[3]->getId();
 247     }
 248
 249     if (type_b  != types::InternalType::IdDouble && type_b  != types::InternalType::IdScalarDouble ||
 250         type_a  != types::InternalType::IdDouble && type_a  != types::InternalType::IdScalarDouble ||
 251         type_x  != types::InternalType::IdDouble && type_x  != types::InternalType::IdScalarDouble ||
 252         type_si != types::InternalType::IdDouble && type_si != types::InternalType::IdScalarDouble)
 253     {
 254         return Overload::call(L"%_filter", in, _iRetCount, out);
 255     }
 256
 257     b = in[0]->getAs<types::Double>();
 258     a = in[1]->getAs<types::Double>();
 259     x = in[2]->getAs<types::Double>();
 260
 261     if (iRhs > 3)
 262     {
 263         si = in[3]->getAs<types::Double>();
 264     }
 265     else
 266     {
 267         alloc_si = true;
 268         int size_a = a->getSize();
 269         int size_b = b->getSize();
 270         int len = std::max(size_a, size_b) - 1;
 271         int leadDim = x->getRows() == 1 ? 1 : 0;
 272
 273         int si_dims[2] = {x->getRows(), x->getCols()};
 274         for (int i = leadDim; i > 0; i--)
 275         {
 276             si_dims[i] = si_dims[i - 1];
 277         }
 278
 279         si_dims[0] = len;
 280
 281         si = new types::Double(2, si_dims);
 282         si->setZeros();
 283     }
 284
 285     if (b->isVector() == false)
 286     {
 287         Scierror(999, _("%s: Wrong size for input argument #%d: Vector expected.\n"), fname, 1);
 288         return types::Function::Error;
 289     }
 290
 291     if (a->isVector() == false)
 292     {
 293         Scierror(999, _("%s: Wrong size for input argument #%d: Vector expected.\n"), fname, 2);
 294         return types::Function::Error;
 295     }
 296
 297     if (x->isVector() == false)
 298     {
 299         Scierror(999, _("%s: Wrong size for input argument #%d: Vector expected.\n"), fname, 3);
 300         return types::Function::Error;
 301     }
 302
 303     if (iRhs > 3 && si->isVector() == false)
 304     {
 305         Scierror(999, _("%s: Wrong size for input argument #%d: Vector expected.\n"), fname, 4);
 306         return types::Function::Error;
 307     }
 308
 309     types::Double* ret = filter(b, a, x, si);
 310     if (_iRetCount != 2 && alloc_si)
 311     {
 312         si->killMe();
 313     }
 314
 315     if (ret == nullptr)
 316     {
 317         return types::Function::Error;
 318     }
 319
 320     out.push_back(ret);
 321     if (_iRetCount == 2)
 322     {
 323         out.push_back(si);
 324     }
 325     return types::Function::OK;
 326 }