GNU Octave 7.1.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
filter.cc
Go to the documentation of this file.
1////////////////////////////////////////////////////////////////////////
2//
3// Copyright (C) 1996-2022 The Octave Project Developers
4//
5// See the file COPYRIGHT.md in the top-level directory of this
6// distribution or <https://octave.org/copyright/>.
7//
8// This file is part of Octave.
9//
10// Octave is free software: you can redistribute it and/or modify it
11// under the terms of the GNU General Public License as published by
12// the Free Software Foundation, either version 3 of the License, or
13// (at your option) any later version.
14//
15// Octave is distributed in the hope that it will be useful, but
16// WITHOUT ANY WARRANTY; without even the implied warranty of
17// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18// GNU General Public License for more details.
19//
20// You should have received a copy of the GNU General Public License
21// along with Octave; see the file COPYING. If not, see
22// <https://www.gnu.org/licenses/>.
23//
24////////////////////////////////////////////////////////////////////////
25
26// Based on Tony Richardson's filter.m.
27//
28// Originally translated to C++ by KH (Kurt.Hornik@wu-wien.ac.at)
29// with help from Fritz Leisch and Andreas Weingessel on Oct 20, 1994.
30//
31// Rewritten to use templates to handle both real and complex cases by
32// jwe, Wed Nov 1 19:15:29 1995.
33
34#if defined (HAVE_CONFIG_H)
35# include "config.h"
36#endif
37
38#include "quit.h"
39
40#include "defun.h"
41#include "error.h"
42#include "ovl.h"
43
44OCTAVE_NAMESPACE_BEGIN
45
46template <typename T>
47MArray<T>
48filter (MArray<T>& b, MArray<T>& a, MArray<T>& x, MArray<T>& si,
49 int dim = 0)
50{
51 MArray<T> y;
52
53 octave_idx_type a_len = a.numel ();
54 octave_idx_type b_len = b.numel ();
55
56 octave_idx_type ab_len = (a_len > b_len ? a_len : b_len);
57
58 // FIXME: The two lines below should be unnecessary because
59 // this template is called with a and b as column vectors
60 // already. However the a.resize line is currently (2011/04/26)
61 // necessary to stop bug #33164.
62 b.resize (dim_vector (ab_len, 1), 0.0);
63 if (a_len > 1)
64 a.resize (dim_vector (ab_len, 1), 0.0);
65
66 T norm = a (0);
67
68 if (norm == static_cast<T> (0.0))
69 error ("filter: the first element of A must be nonzero");
70
71 dim_vector x_dims = x.dims ();
72 if (dim < 0 || dim > x_dims.ndims ())
73 error ("filter: DIM must be a valid dimension");
74
75 octave_idx_type x_len = x_dims(dim);
76
77 dim_vector si_dims = si.dims ();
78 octave_idx_type si_len = si_dims(0);
79
80 if (si_len != ab_len - 1)
81 error ("filter: first dimension of SI must be of length max (length (a), length (b)) - 1");
82
83 if (si_dims.ndims () != x_dims.ndims ())
84 error ("filter: dimensionality of SI and X must agree");
85
86 for (octave_idx_type i = 1; i < dim; i++)
87 {
88 if (si_dims(i) != x_dims(i-1))
89 error ("filter: dimensionality of SI and X must agree");
90 }
91 for (octave_idx_type i = dim+1; i < x_dims.ndims (); i++)
92 {
93 if (si_dims(i) != x_dims(i))
94 error ("filter: dimensionality of SI and X must agree");
95 }
96
97 if (x_len == 0)
98 return x;
99
100 if (norm != static_cast<T> (1.0))
101 {
102 a /= norm;
103 b /= norm;
104 }
105
106 if (a_len <= 1 && si_len <= 0)
107 return b(0) * x;
108
109 y.resize (x_dims, 0.0);
110
111 octave_idx_type x_stride = 1;
112 for (int i = 0; i < dim; i++)
113 x_stride *= x_dims(i);
114
115 octave_idx_type x_num = x_dims.numel () / x_len;
116 for (octave_idx_type num = 0; num < x_num; num++)
117 {
118 octave_idx_type x_offset;
119 if (x_stride == 1)
120 x_offset = num * x_len;
121 else
122 {
123 x_offset = num;
124 octave_idx_type n_strides = num / x_stride;
125 x_offset += n_strides * x_stride * (x_len - 1);
126 }
127 octave_idx_type si_offset = num * si_len;
128
129 if (a_len > 1)
130 {
131 T *py = y.fortran_vec ();
132 T *psi = si.fortran_vec ();
133
134 const T *pa = a.data ();
135 const T *pb = b.data ();
136 const T *px = x.data ();
137
138 psi += si_offset;
139
140 for (octave_idx_type i = 0, idx = x_offset;
141 i < x_len;
142 i++, idx += x_stride)
143 {
144 py[idx] = psi[0] + pb[0] * px[idx];
145
146 if (si_len > 0)
147 {
148 for (octave_idx_type j = 0; j < si_len - 1; j++)
149 {
150 octave_quit ();
151
152 psi[j] = psi[j+1] - pa[j+1] * py[idx] + pb[j+1] * px[idx];
153 }
154
155 psi[si_len-1] = pb[si_len] * px[idx] - pa[si_len] * py[idx];
156 }
157 else
158 {
159 octave_quit ();
160
161 psi[0] = pb[si_len] * px[idx] - pa[si_len] * py[idx];
162 }
163 }
164 }
165 else if (si_len > 0)
166 {
167 T *py = y.fortran_vec ();
168 T *psi = si.fortran_vec ();
169
170 const T *pb = b.data ();
171 const T *px = x.data ();
172
173 psi += si_offset;
174
175 for (octave_idx_type i = 0, idx = x_offset;
176 i < x_len;
177 i++, idx += x_stride)
178 {
179 py[idx] = psi[0] + pb[0] * px[idx];
180
181 if (si_len > 1)
182 {
183 for (octave_idx_type j = 0; j < si_len - 1; j++)
184 {
185 octave_quit ();
186
187 psi[j] = psi[j+1] + pb[j+1] * px[idx];
188 }
189
190 psi[si_len-1] = pb[si_len] * px[idx];
191 }
192 else
193 {
194 octave_quit ();
195
196 psi[0] = pb[1] * px[idx];
197 }
198 }
199 }
200 }
201
202 return y;
203}
204
205template <typename T>
206MArray<T>
207filter (MArray<T>& b, MArray<T>& a, MArray<T>& x, int dim = -1)
208{
209 dim_vector x_dims = x.dims ();
210
211 if (dim < 0)
212 dim = x_dims.first_non_singleton ();
213 else if (dim > x_dims.ndims ())
214 error ("filter: DIM must be a valid dimension");
215
216 octave_idx_type a_len = a.numel ();
217 octave_idx_type b_len = b.numel ();
218
219 octave_idx_type si_len = (a_len > b_len ? a_len : b_len) - 1;
220 dim_vector si_dims = x.dims ();
221 for (int i = dim; i > 0; i--)
222 si_dims(i) = si_dims(i-1);
223 si_dims(0) = si_len;
224
225 MArray<T> si (si_dims, T (0.0));
226
227 return filter (b, a, x, si, dim);
228}
229
230DEFUN (filter, args, ,
231 doc: /* -*- texinfo -*-
232@deftypefn {} {@var{y} =} filter (@var{b}, @var{a}, @var{x})
233@deftypefnx {} {[@var{y}, @var{sf}] =} filter (@var{b}, @var{a}, @var{x}, @var{si})
234@deftypefnx {} {[@var{y}, @var{sf}] =} filter (@var{b}, @var{a}, @var{x}, [], @var{dim})
235@deftypefnx {} {[@var{y}, @var{sf}] =} filter (@var{b}, @var{a}, @var{x}, @var{si}, @var{dim})
236Apply a 1-D digital filter to the data @var{x}.
237
238@code{filter} returns the solution to the following linear, time-invariant
239difference equation:
240@tex
241$$
242\sum_{k=0}^N a_{k+1} y_{n-k} = \sum_{k=0}^M b_{k+1} x_{n-k}, \qquad
243 1 \le n \le P
244$$
245@end tex
246@ifnottex
247@c Set example in small font to prevent overfull line
248
249@smallexample
250@group
251 N M
252SUM a(k+1) y(n-k) = SUM b(k+1) x(n-k) for 1<=n<=length(x)
253k=0 k=0
254@end group
255@end smallexample
256
257@end ifnottex
258
259@noindent
260where
261@ifnottex
262N=length(a)-1 and M=length(b)-1.
263@end ifnottex
264@tex
265$a \in \Re^{N-1}$, $b \in \Re^{M-1}$, and $x \in \Re^P$.
266@end tex
267The result is calculated over the first non-singleton dimension of @var{x}
268or over @var{dim} if supplied.
269
270An equivalent form of the equation is:
271@tex
272$$
273y_n = -\sum_{k=1}^N c_{k+1} y_{n-k} + \sum_{k=0}^M d_{k+1} x_{n-k}, \qquad
274 1 \le n \le P
275$$
276@end tex
277@ifnottex
278@c Set example in small font to prevent overfull line
279
280@smallexample
281@group
282 N M
283y(n) = - SUM c(k+1) y(n-k) + SUM d(k+1) x(n-k) for 1<=n<=length(x)
284 k=1 k=0
285@end group
286@end smallexample
287
288@end ifnottex
289
290@noindent
291where
292@ifnottex
293 c = a/a(1) and d = b/a(1).
294@end ifnottex
295@tex
296$c = a/a_1$ and $d = b/a_1$.
297@end tex
298
299If the fourth argument @var{si} is provided, it is taken as the
300initial state of the system and the final state is returned as
301@var{sf}. The state vector is a column vector whose length is
302equal to the length of the longest coefficient vector minus one.
303If @var{si} is not supplied, the initial state vector is set to all
304zeros.
305
306In terms of the Z Transform, @var{y} is the result of passing the
307discrete-time signal @var{x} through a system characterized by the following
308rational system function:
309@tex
310$$
311H(z) = {\displaystyle\sum_{k=0}^M d_{k+1} z^{-k}
312 \over 1 + \displaystyle\sum_{k+1}^N c_{k+1} z^{-k}}
313$$
314@end tex
315@ifnottex
316
317@example
318@group
319 M
320 SUM d(k+1) z^(-k)
321 k=0
322H(z) = ---------------------
323 N
324 1 + SUM c(k+1) z^(-k)
325 k=1
326@end group
327@end example
328
329@end ifnottex
330@seealso{filter2, fftfilt, freqz}
331@end deftypefn */)
332{
333 int nargin = args.length ();
334
335 if (nargin < 3 || nargin > 5)
336 print_usage ();
337
338 int dim;
339 dim_vector x_dims = args(2).dims ();
340
341 if (nargin == 5)
342 {
343 dim = args(4).nint_value () - 1;
344 if (dim < 0 || dim >= x_dims.ndims ())
345 error ("filter: DIM must be a valid dimension");
346 }
347 else
348 dim = x_dims.first_non_singleton ();
349
350 octave_value_list retval;
351
352 const char *a_b_errmsg = "filter: A and B must be vectors";
353 const char *x_si_errmsg = "filter: X and SI must be arrays";
354
355 bool isfloat = (args(0).is_single_type ()
356 || args(1).is_single_type ()
357 || args(2).is_single_type ()
358 || (nargin >= 4 && args(3).is_single_type ()));
359
360 if (args(0).iscomplex ()
361 || args(1).iscomplex ()
362 || args(2).iscomplex ()
363 || (nargin >= 4 && args(3).iscomplex ()))
364 {
365 if (isfloat)
366 {
367 FloatComplexColumnVector b = args(0).xfloat_complex_vector_value (a_b_errmsg);
368 FloatComplexColumnVector a = args(1).xfloat_complex_vector_value (a_b_errmsg);
369 FloatComplexNDArray x = args(2).xfloat_complex_array_value (x_si_errmsg);
370
371 FloatComplexNDArray si;
372
373 if (nargin == 3 || args(3).isempty ())
374 {
375 octave_idx_type a_len = a.numel ();
376 octave_idx_type b_len = b.numel ();
377
378 octave_idx_type si_len = (a_len > b_len ? a_len : b_len) - 1;
379
380 dim_vector si_dims = x.dims ();
381 for (int i = dim; i > 0; i--)
382 si_dims(i) = si_dims(i-1);
383 si_dims(0) = si_len;
384
385 si.resize (si_dims, 0.0);
386 }
387 else
388 {
389 si = args(3).xfloat_complex_array_value (x_si_errmsg);
390
391 if (si.isvector () && x.isvector ())
392 si = si.reshape (dim_vector (si.numel (), 1));
393 }
394
395 FloatComplexNDArray y (filter (b, a, x, si, dim));
396
397 retval = ovl (y, si);
398 }
399 else
400 {
401 ComplexColumnVector b = args(0).xcomplex_vector_value (a_b_errmsg);
402 ComplexColumnVector a = args(1).xcomplex_vector_value (a_b_errmsg);
403
404 ComplexNDArray x = args(2).xcomplex_array_value (x_si_errmsg);
405
406 ComplexNDArray si;
407
408 if (nargin == 3 || args(3).isempty ())
409 {
410 octave_idx_type a_len = a.numel ();
411 octave_idx_type b_len = b.numel ();
412
413 octave_idx_type si_len = (a_len > b_len ? a_len : b_len) - 1;
414
415 dim_vector si_dims = x.dims ();
416 for (int i = dim; i > 0; i--)
417 si_dims(i) = si_dims(i-1);
418 si_dims(0) = si_len;
419
420 si.resize (si_dims, 0.0);
421 }
422 else
423 {
424 si = args(3).xcomplex_array_value (x_si_errmsg);
425
426 if (si.isvector () && x.isvector ())
427 si = si.reshape (dim_vector (si.numel (), 1));
428 }
429
430 ComplexNDArray y (filter (b, a, x, si, dim));
431
432 retval = ovl (y, si);
433 }
434 }
435 else
436 {
437 if (isfloat)
438 {
439 FloatColumnVector b = args(0).xfloat_vector_value (a_b_errmsg);
440 FloatColumnVector a = args(1).xfloat_vector_value (a_b_errmsg);
441
442 FloatNDArray x = args(2).xfloat_array_value (x_si_errmsg);
443
444 FloatNDArray si;
445
446 if (nargin == 3 || args(3).isempty ())
447 {
448 octave_idx_type a_len = a.numel ();
449 octave_idx_type b_len = b.numel ();
450
451 octave_idx_type si_len = (a_len > b_len ? a_len : b_len) - 1;
452
453 dim_vector si_dims = x.dims ();
454 for (int i = dim; i > 0; i--)
455 si_dims(i) = si_dims(i-1);
456 si_dims(0) = si_len;
457
458 si.resize (si_dims, 0.0);
459 }
460 else
461 {
462 si = args(3).xfloat_array_value (x_si_errmsg);
463
464 if (si.isvector () && x.isvector ())
465 si = si.reshape (dim_vector (si.numel (), 1));
466 }
467
468 FloatNDArray y (filter (b, a, x, si, dim));
469
470 retval = ovl (y, si);
471 }
472 else
473 {
474 ColumnVector b = args(0).xvector_value (a_b_errmsg);
475 ColumnVector a = args(1).xvector_value (a_b_errmsg);
476
477 NDArray x = args(2).xarray_value (x_si_errmsg);
478
479 NDArray si;
480
481 if (nargin == 3 || args(3).isempty ())
482 {
483 octave_idx_type a_len = a.numel ();
484 octave_idx_type b_len = b.numel ();
485
486 octave_idx_type si_len = (a_len > b_len ? a_len : b_len) - 1;
487
488 dim_vector si_dims = x.dims ();
489 for (int i = dim; i > 0; i--)
490 si_dims(i) = si_dims(i-1);
491 si_dims(0) = si_len;
492
493 si.resize (si_dims, 0.0);
494 }
495 else
496 {
497 si = args(3).xarray_value (x_si_errmsg);
498
499 if (si.isvector () && x.isvector ())
500 si = si.reshape (dim_vector (si.numel (), 1));
501 }
502
503 NDArray y (filter (b, a, x, si, dim));
504
505 retval = ovl (y, si);
506 }
507 }
508
509 return retval;
510}
511
512template MArray<double>
513filter (MArray<double>&, MArray<double>&, MArray<double>&,
514 MArray<double>&, int dim);
515
516template MArray<double>
517filter (MArray<double>&, MArray<double>&, MArray<double>&, int dim);
518
519template MArray<Complex>
520filter (MArray<Complex>&, MArray<Complex>&, MArray<Complex>&,
521 MArray<Complex>&, int dim);
522
523template MArray<Complex>
524filter (MArray<Complex>&, MArray<Complex>&, MArray<Complex>&, int dim);
525
526template MArray<float>
527filter (MArray<float>&, MArray<float>&, MArray<float>&,
528 MArray<float>&, int dim);
529
530template MArray<float>
531filter (MArray<float>&, MArray<float>&, MArray<float>&, int dim);
532
533template MArray<FloatComplex>
534filter (MArray<FloatComplex>&, MArray<FloatComplex>&, MArray<FloatComplex>&,
535 MArray<FloatComplex>&, int dim);
536
537template MArray<FloatComplex>
538filter (MArray<FloatComplex>&, MArray<FloatComplex>&, MArray<FloatComplex>&,
539 int dim);
540
541/*
542%!shared a, b, x, r
543%!test
544%! a = [1 1];
545%! b = [1 1];
546%! x = zeros (1,10); x(1) = 1;
547%! assert (filter (b, [1], x ), [1 1 0 0 0 0 0 0 0 0]);
548%! assert (filter (b, [1], x.'), [1 1 0 0 0 0 0 0 0 0].');
549%! assert (filter (b.', [1], x ), [1 1 0 0 0 0 0 0 0 0] );
550%! assert (filter (b.', [1], x.'), [1 1 0 0 0 0 0 0 0 0].');
551%! assert (filter ([1], a, x ), [+1 -1 +1 -1 +1 -1 +1 -1 +1 -1] );
552%! assert (filter ([1], a, x.'), [+1 -1 +1 -1 +1 -1 +1 -1 +1 -1].');
553%! assert (filter ([1], a.', x ), [+1 -1 +1 -1 +1 -1 +1 -1 +1 -1] );
554%! assert (filter ([1], a.', x.'), [+1 -1 +1 -1 +1 -1 +1 -1 +1 -1].');
555%! assert (filter (b, a, x ), [1 0 0 0 0 0 0 0 0 0] );
556%! assert (filter (b.', a, x ), [1 0 0 0 0 0 0 0 0 0] );
557%! assert (filter (b, a.', x ), [1 0 0 0 0 0 0 0 0 0] );
558%! assert (filter (b.', a, x ), [1 0 0 0 0 0 0 0 0 0] );
559%! assert (filter (b, a, x.'), [1 0 0 0 0 0 0 0 0 0].');
560%! assert (filter (b.', a, x.'), [1 0 0 0 0 0 0 0 0 0].');
561%! assert (filter (b, a.', x.'), [1 0 0 0 0 0 0 0 0 0].');
562%! assert (filter (b.', a, x.'), [1 0 0 0 0 0 0 0 0 0].');
563
564%!test
565%! r = sqrt (1/2) * (1+i);
566%! a = a*r;
567%! b = b*r;
568%! assert (filter (b, [1], x ), r*[1 1 0 0 0 0 0 0 0 0] );
569%! assert (filter (b, [1], r*x ), r*r*[1 1 0 0 0 0 0 0 0 0] );
570%! assert (filter (b, [1], x.' ), r*[1 1 0 0 0 0 0 0 0 0].' );
571%! assert (filter (b, a, x ), [1 0 0 0 0 0 0 0 0 0] );
572%! assert (filter (b, a, r*x ), r*[1 0 0 0 0 0 0 0 0 0] );
573
574%!shared a, b, x, y, so
575%!test
576%! a = [1,1];
577%! b = [1,1];
578%! x = zeros (1,10); x(1) = 1;
579%! [y, so] = filter (b, [1], x, [-1]);
580%! assert (y, [0 1 0 0 0 0 0 0 0 0]);
581%! assert (so, 0);
582
583%!test
584%! x = zeros (10,3); x(1,1) = -1; x(1,2) = 1;
585%! y0 = zeros (10,3); y0(1:2,1) = -1; y0(1:2,2) = 1;
586%! y = filter (b, [1], x);
587%! assert (y, y0);
588
589%!test
590%! a = [1,1];
591%! b=[1,1];
592%! x = zeros (4,4,2); x(1,1:4,1) = +1; x(1,1:4,2) = -1;
593%! y0 = zeros (4,4,2); y0(1:2,1:4,1) = +1; y0(1:2,1:4,2) = -1;
594%! y = filter (b, [1], x);
595%! assert (y, y0);
596
597%!assert (filter (1, ones (10,1) / 10, []), [])
598%!assert (filter (1, ones (10,1) / 10, zeros (0,10)), zeros (0,10))
599%!assert (filter (1, ones (10,1) / 10, single (1:5)), repmat (single (10), 1, 5))
600
601## Test using initial conditions
602%!assert (filter ([1, 1, 1], [1, 1], [1 2], [1, 1]), [2 2])
603%!assert (filter ([1, 1, 1], [1, 1], [1 2], [1, 1]'), [2 2])
604%!assert (filter ([1, 3], [1], [1 2; 3 4; 5 6], [4, 5]), [5 7; 6 10; 14 18])
605%!error filter ([1, 3], [1], [1 2; 3 4; 5 6], [4, 5]')
606%!assert (filter ([1, 3, 2], [1], [1 2; 3 4; 5 6], [1 0 0; 1 0 0], 2), [2 6; 3 13; 5 21])
607
608## Test of DIM parameter
609%!test
610%! x = ones (2, 1, 3, 4);
611%! x(1,1,:,:) = [1 2 3 4; 5 6 7 8; 9 10 11 12];
612%! y0 = [1 1 6 2 15 3 2 1 8 2 18 3 3 1 10 2 21 3 4 1 12 2 24 3];
613%! y0 = reshape (y0, size (x));
614%! y = filter ([1 1 1], 1, x, [], 3);
615%! assert (y, y0);
616*/
617
618OCTAVE_NAMESPACE_END
Array::numel
octave_idx_type numel(void) const
Number of elements in the array.
Definition: Array.h:411
Array::dims
const dim_vector & dims(void) const
Return a const-reference so that dims ()(i) works efficiently.
Definition: Array.h:487
Array::isvector
bool isvector(void) const
Size of the specified dimension.
Definition: Array.h:609
Array::resize
OCTARRAY_API void resize(const dim_vector &dv, const T &rfv)
Size of the specified dimension.
Definition: Array.cc:1010
Array::data
const T * data(void) const
Size of the specified dimension.
Definition: Array.h:616
Array::fortran_vec
OCTARRAY_API T * fortran_vec(void)
Size of the specified dimension.
Definition: Array.cc:1744
MArray
Template for N-dimensional array classes with like-type math operators.
Definition: MArray.h:63
MArray::reshape
MArray< T > reshape(const dim_vector &new_dims) const
Definition: MArray.h:87
dim_vector
Vector representing the dimensions (size) of an Array.
Definition: dim-vector.h:94
dim_vector::numel
octave_idx_type numel(int n=0) const
Number of elements that a matrix with this dimensions would have.
Definition: dim-vector.h:335
dim_vector::first_non_singleton
int first_non_singleton(int def=0) const
Definition: dim-vector.h:444
dim_vector::ndims
octave_idx_type ndims(void) const
Number of dimensions.
Definition: dim-vector.h:257
octave_idx_type
octave_value_list
Definition: ovl.h:44
print_usage
OCTINTERP_API void print_usage(void)
Definition: defun-int.h:72
DEFUN
#define DEFUN(name, args_name, nargout_name, doc)
Macro to define a builtin function.
Definition: defun.h:56
error
void error(const char *fmt,...)
Definition: error.cc:980
filter
OCTAVE_NAMESPACE_BEGIN MArray< T > filter(MArray< T > &b, MArray< T > &a, MArray< T > &x, MArray< T > &si, int dim=0)
Definition: filter.cc:48
norm
double norm(const ColumnVector &v)
Definition: graphics.cc:5940
x
F77_RET_T const F77_DBLE * x
Definition: lo-slatec-proto.h:38
octave::math::psi
double psi(double z)
Definition: lo-specfun.cc:2099
ovl
octave_value_list ovl(const OV_Args &... args)
Construct an octave_value_list with less typing.
Definition: ovl.h:211