GNU Octave  3.8.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
filter.cc
Go to the documentation of this file.
1 /*
2 
3 Copyright (C) 1996-2013 John W. Eaton
4 
5 This file is part of Octave.
6 
7 Octave is free software; you can redistribute it and/or modify it
8 under the terms of the GNU General Public License as published by the
9 Free Software Foundation; either version 3 of the License, or (at your
10 option) any later version.
11 
12 Octave is distributed in the hope that it will be useful, but WITHOUT
13 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
14 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
15 for more details.
16 
17 You should have received a copy of the GNU General Public License
18 along with Octave; see the file COPYING. If not, see
19 <http://www.gnu.org/licenses/>.
20 
21 */
22 
23 // Based on Tony Richardson's filter.m.
24 //
25 // Originally translated to C++ by KH (Kurt.Hornik@wu-wien.ac.at)
26 // with help from Fritz Leisch and Andreas Weingessel on Oct 20, 1994.
27 //
28 // Rewritten to use templates to handle both real and complex cases by
29 // jwe, Wed Nov 1 19:15:29 1995.
30 
31 #ifdef HAVE_CONFIG_H
32 #include <config.h>
33 #endif
34 
35 #include "quit.h"
36 
37 #include "defun.h"
38 #include "error.h"
39 #include "oct-obj.h"
40 
41 #if !defined (CXX_NEW_FRIEND_TEMPLATE_DECL)
42 extern MArray<double>
43 filter (MArray<double>&, MArray<double>&, MArray<double>&, int dim);
44 
45 extern MArray<Complex>
46 filter (MArray<Complex>&, MArray<Complex>&, MArray<Complex>&, int dim);
47 
48 extern MArray<float>
49 filter (MArray<float>&, MArray<float>&, MArray<float>&, int dim);
50 
51 extern MArray<FloatComplex>
52 filter (MArray<FloatComplex>&, MArray<FloatComplex>&, MArray<FloatComplex>&,
53  int dim);
54 #endif
55 
56 template <class T>
57 MArray<T>
58 filter (MArray<T>& b, MArray<T>& a, MArray<T>& x, MArray<T>& si,
59  int dim = 0)
60 {
61  MArray<T> y;
62 
63  octave_idx_type a_len = a.length ();
64  octave_idx_type b_len = b.length ();
65 
66  octave_idx_type ab_len = a_len > b_len ? a_len : b_len;
67 
68  // FIXME: The two lines below should be unecessary because
69  // this template is called with a and b as column vectors
70  // already. However the a.resize line is currently (2011/04/26)
71  // necessary to stop bug #33164.
72  b.resize (dim_vector (ab_len, 1), 0.0);
73  if (a_len > 1)
74  a.resize (dim_vector (ab_len, 1), 0.0);
75 
76  T norm = a (0);
77 
78  if (norm == static_cast<T>(0.0))
79  {
80  error ("filter: the first element of A must be non-zero");
81  return y;
82  }
83 
84  dim_vector x_dims = x.dims ();
85  if (dim < 0 || dim > x_dims.length ())
86  {
87  error ("filter: DIM must be a valid dimension");
88  return y;
89  }
90 
91  octave_idx_type x_len = x_dims(dim);
92 
93  dim_vector si_dims = si.dims ();
94  octave_idx_type si_len = si_dims(0);
95 
96  if (si_len != ab_len - 1)
97  {
98  error ("filter: first dimension of SI must be of length max (length (a), length (b)) - 1");
99  return y;
100  }
101 
102  if (si_dims.length () != x_dims.length ())
103  {
104  error ("filter: dimensionality of SI and X must agree");
105  return y;
106  }
107 
108  for (octave_idx_type i = 1; i < dim; i++)
109  {
110  if (si_dims(i) != x_dims(i-1))
111  {
112  error ("filter: dimensionality of SI and X must agree");
113  return y;
114  }
115  }
116  for (octave_idx_type i = dim+1; i < x_dims.length (); i++)
117  {
118  if (si_dims(i) != x_dims(i))
119  {
120  error ("filter: dimensionality of SI and X must agree");
121  return y;
122  }
123  }
124 
125  if (x_len == 0)
126  return x;
127 
128  if (norm != static_cast<T>(1.0))
129  {
130  a = a / norm;
131  b = b / norm;
132  }
133 
134  if (a_len <= 1 && si_len <= 0)
135  return b(0) * x;
136 
137  y.resize (x_dims, 0.0);
138 
139  int x_stride = 1;
140  for (int i = 0; i < dim; i++)
141  x_stride *= x_dims(i);
142 
143  octave_idx_type x_num = x_dims.numel () / x_len;
144  for (octave_idx_type num = 0; num < x_num; num++)
145  {
146  octave_idx_type x_offset;
147  if (x_stride == 1)
148  x_offset = num * x_len;
149  else
150  {
151  octave_idx_type x_offset2 = 0;
152  x_offset = num;
153  while (x_offset >= x_stride)
154  {
155  x_offset -= x_stride;
156  x_offset2++;
157  }
158  x_offset += x_offset2 * x_stride * x_len;
159  }
160  octave_idx_type si_offset = num * si_len;
161 
162  if (a_len > 1)
163  {
164  T *py = y.fortran_vec ();
165  T *psi = si.fortran_vec ();
166 
167  const T *pa = a.data ();
168  const T *pb = b.data ();
169  const T *px = x.data ();
170 
171  psi += si_offset;
172 
173  for (octave_idx_type i = 0, idx = x_offset;
174  i < x_len;
175  i++, idx += x_stride)
176  {
177  py[idx] = psi[0] + pb[0] * px[idx];
178 
179  if (si_len > 0)
180  {
181  for (octave_idx_type j = 0; j < si_len - 1; j++)
182  {
183  OCTAVE_QUIT;
184 
185  psi[j] = psi[j+1] - pa[j+1] * py[idx] + pb[j+1] * px[idx];
186  }
187 
188  psi[si_len-1] = pb[si_len] * px[idx] - pa[si_len] * py[idx];
189  }
190  else
191  {
192  OCTAVE_QUIT;
193 
194  psi[0] = pb[si_len] * px[idx] - pa[si_len] * py[idx];
195  }
196  }
197  }
198  else if (si_len > 0)
199  {
200  T *py = y.fortran_vec ();
201  T *psi = si.fortran_vec ();
202 
203  const T *pb = b.data ();
204  const T *px = x.data ();
205 
206  psi += si_offset;
207 
208  for (octave_idx_type i = 0, idx = x_offset;
209  i < x_len;
210  i++, idx += x_stride)
211  {
212  py[idx] = psi[0] + pb[0] * px[idx];
213 
214  if (si_len > 1)
215  {
216  for (octave_idx_type j = 0; j < si_len - 1; j++)
217  {
218  OCTAVE_QUIT;
219 
220  psi[j] = psi[j+1] + pb[j+1] * px[idx];
221  }
222 
223  psi[si_len-1] = pb[si_len] * px[idx];
224  }
225  else
226  {
227  OCTAVE_QUIT;
228 
229  psi[0] = pb[1] * px[idx];
230  }
231  }
232  }
233  }
234 
235  return y;
236 }
237 
238 #if !defined (CXX_NEW_FRIEND_TEMPLATE_DECL)
239 extern MArray<double>
240 filter (MArray<double>&, MArray<double>&, MArray<double>&,
241  MArray<double>&, int dim);
242 
243 extern MArray<Complex>
244 filter (MArray<Complex>&, MArray<Complex>&, MArray<Complex>&,
245  MArray<Complex>&, int dim);
246 
247 extern MArray<float>
248 filter (MArray<float>&, MArray<float>&, MArray<float>&,
249  MArray<float>&, int dim);
250 
251 extern MArray<FloatComplex>
252 filter (MArray<FloatComplex>&, MArray<FloatComplex>&, MArray<FloatComplex>&,
253  MArray<FloatComplex>&, int dim);
254 #endif
255 
256 template <class T>
257 MArray<T>
258 filter (MArray<T>& b, MArray<T>& a, MArray<T>& x, int dim = -1)
259 {
260  dim_vector x_dims = x.dims ();
261 
262  if (dim < 0)
263  {
264  // Find first non-singleton dimension
265  while (dim < x_dims.length () && x_dims(dim) <= 1)
266  dim++;
267 
268  // All dimensions singleton, pick first dimension
269  if (dim == x_dims.length ())
270  dim = 0;
271  }
272  else if (dim < 0 || dim > x_dims.length ())
273  {
274  error ("filter: DIM must be a valid dimension");
275  return MArray<T> ();
276  }
277 
278  octave_idx_type a_len = a.length ();
279  octave_idx_type b_len = b.length ();
280 
281  octave_idx_type si_len = (a_len > b_len ? a_len : b_len) - 1;
282  dim_vector si_dims = x.dims ();
283  for (int i = dim; i > 0; i--)
284  si_dims(i) = si_dims(i-1);
285  si_dims(0) = si_len;
286 
287  MArray<T> si (si_dims, T (0.0));
288 
289  return filter (b, a, x, si, dim);
290 }
291 
292 DEFUN (filter, args, nargout,
293  "-*- texinfo -*-\n\
294 @deftypefn {Built-in Function} {y =} filter (@var{b}, @var{a}, @var{x})\n\
295 @deftypefnx {Built-in Function} {[@var{y}, @var{sf}] =} filter (@var{b}, @var{a}, @var{x}, @var{si})\n\
296 @deftypefnx {Built-in Function} {[@var{y}, @var{sf}] =} filter (@var{b}, @var{a}, @var{x}, [], @var{dim})\n\
297 @deftypefnx {Built-in Function} {[@var{y}, @var{sf}] =} filter (@var{b}, @var{a}, @var{x}, @var{si}, @var{dim})\n\
298 Return the solution to the following linear, time-invariant difference\n\
299 equation:\n\
300 @tex\n\
301 $$\n\
302 \\sum_{k=0}^N a_{k+1} y_{n-k} = \\sum_{k=0}^M b_{k+1} x_{n-k}, \\qquad\n\
303  1 \\le n \\le P\n\
304 $$\n\
305 @end tex\n\
306 @ifnottex\n\
307 @c Set example in small font to prevent overfull line\n\
308 \n\
309 @smallexample\n\
310 @group\n\
311  N M\n\
312 SUM a(k+1) y(n-k) = SUM b(k+1) x(n-k) for 1<=n<=length(x)\n\
313 k=0 k=0\n\
314 @end group\n\
315 @end smallexample\n\
316 \n\
317 @end ifnottex\n\
318 \n\
319 @noindent\n\
320 where\n\
321 @ifnottex\n\
322 N=length(a)-1 and M=length(b)-1.\n\
323 @end ifnottex\n\
324 @tex\n\
325 $a \\in \\Re^{N-1}$, $b \\in \\Re^{M-1}$, and $x \\in \\Re^P$.\n\
326 @end tex\n\
327 The result is calculated over the first non-singleton dimension of @var{x}\n\
328 or over @var{dim} if supplied.\n\
329 \n\
330 An equivalent form of the equation is:\n\
331 @tex\n\
332 $$\n\
333 y_n = -\\sum_{k=1}^N c_{k+1} y_{n-k} + \\sum_{k=0}^M d_{k+1} x_{n-k}, \\qquad\n\
334  1 \\le n \\le P\n\
335 $$\n\
336 @end tex\n\
337 @ifnottex\n\
338 @c Set example in small font to prevent overfull line\n\
339 \n\
340 @smallexample\n\
341 @group\n\
342  N M\n\
343 y(n) = - SUM c(k+1) y(n-k) + SUM d(k+1) x(n-k) for 1<=n<=length(x)\n\
344  k=1 k=0\n\
345 @end group\n\
346 @end smallexample\n\
347 \n\
348 @end ifnottex\n\
349 \n\
350 @noindent\n\
351 where\n\
352 @ifnottex\n\
353  c = a/a(1) and d = b/a(1).\n\
354 @end ifnottex\n\
355 @tex\n\
356 $c = a/a_1$ and $d = b/a_1$.\n\
357 @end tex\n\
358 \n\
359 If the fourth argument @var{si} is provided, it is taken as the\n\
360 initial state of the system and the final state is returned as\n\
361 @var{sf}. The state vector is a column vector whose length is\n\
362 equal to the length of the longest coefficient vector minus one.\n\
363 If @var{si} is not supplied, the initial state vector is set to all\n\
364 zeros.\n\
365 \n\
366 In terms of the Z Transform, y is the result of passing the discrete-\n\
367 time signal x through a system characterized by the following rational\n\
368 system function:\n\
369 @tex\n\
370 $$\n\
371 H(z) = {\\displaystyle\\sum_{k=0}^M d_{k+1} z^{-k}\n\
372  \\over 1 + \\displaystyle\\sum_{k+1}^N c_{k+1} z^{-k}}\n\
373 $$\n\
374 @end tex\n\
375 @ifnottex\n\
376 \n\
377 @example\n\
378 @group\n\
379  M\n\
380  SUM d(k+1) z^(-k)\n\
381  k=0\n\
382 H(z) = ---------------------\n\
383  N\n\
384  1 + SUM c(k+1) z^(-k)\n\
385  k=1\n\
386 @end group\n\
387 @end example\n\
388 \n\
389 @end ifnottex\n\
390 @seealso{filter2, fftfilt, freqz}\n\
391 @end deftypefn")
392 {
393  octave_value_list retval;
394 
395  int nargin = args.length ();
396 
397  if (nargin < 3 || nargin > 5)
398  {
399  print_usage ();
400  return retval;
401  }
402 
403  const char *errmsg = "filter: arguments a and b must be vectors";
404 
405  int dim;
406  dim_vector x_dims = args(2).dims ();
407 
408  if (nargin == 5)
409  {
410  dim = args(4).nint_value () - 1;
411  if (dim < 0 || dim >= x_dims.length ())
412  {
413  error ("filter: DIM must be a valid dimension");
414  return retval;
415  }
416  }
417  else
418  {
419  // Find first non-singleton dimension
420  dim = 0;
421  while (dim < x_dims.length () && x_dims(dim) <= 1)
422  dim++;
423 
424  // All dimensions singleton, pick first dimension
425  if (dim == x_dims.length ())
426  dim = 0;
427  }
428 
429  bool isfloat = (args(0).is_single_type ()
430  || args(1).is_single_type ()
431  || args(2).is_single_type ()
432  || (nargin >= 4 && args(3).is_single_type ()));
433 
434  if (args(0).is_complex_type ()
435  || args(1).is_complex_type ()
436  || args(2).is_complex_type ()
437  || (nargin >= 4 && args(3).is_complex_type ()))
438  {
439  if (isfloat)
440  {
441  FloatComplexColumnVector b (args(0).float_complex_vector_value ());
442  FloatComplexColumnVector a (args(1).float_complex_vector_value ());
443 
444  FloatComplexNDArray x (args(2).float_complex_array_value ());
445 
446  if (! error_state)
447  {
448  FloatComplexNDArray si;
449 
450  if (nargin == 3 || args(3).is_empty ())
451  {
452  octave_idx_type a_len = a.length ();
453  octave_idx_type b_len = b.length ();
454 
455  octave_idx_type si_len = (a_len > b_len ? a_len : b_len) - 1;
456 
457  dim_vector si_dims = x.dims ();
458  for (int i = dim; i > 0; i--)
459  si_dims(i) = si_dims(i-1);
460  si_dims(0) = si_len;
461 
462  si.resize (si_dims, 0.0);
463  }
464  else
465  {
466  si = args(3).float_complex_array_value ();
467 
468  if (si.is_vector () && x.is_vector ())
469  si = si.reshape (dim_vector (si.numel (), 1));
470  }
471 
472  if (! error_state)
473  {
474  FloatComplexNDArray y (filter (b, a, x, si, dim));
475 
476  if (nargout == 2)
477  retval(1) = si;
478 
479  retval(0) = y;
480  }
481  else
482  error (errmsg);
483  }
484  else
485  error (errmsg);
486  }
487  else
488  {
489  ComplexColumnVector b (args(0).complex_vector_value ());
490  ComplexColumnVector a (args(1).complex_vector_value ());
491 
492  ComplexNDArray x (args(2).complex_array_value ());
493 
494  if (! error_state)
495  {
496  ComplexNDArray si;
497 
498  if (nargin == 3 || args(3).is_empty ())
499  {
500  octave_idx_type a_len = a.length ();
501  octave_idx_type b_len = b.length ();
502 
503  octave_idx_type si_len = (a_len > b_len ? a_len : b_len) - 1;
504 
505  dim_vector si_dims = x.dims ();
506  for (int i = dim; i > 0; i--)
507  si_dims(i) = si_dims(i-1);
508  si_dims(0) = si_len;
509 
510  si.resize (si_dims, 0.0);
511  }
512  else
513  {
514  si = args(3).complex_array_value ();
515 
516  if (si.is_vector () && x.is_vector ())
517  si = si.reshape (dim_vector (si.numel (), 1));
518  }
519 
520  if (! error_state)
521  {
522  ComplexNDArray y (filter (b, a, x, si, dim));
523 
524  if (nargout == 2)
525  retval(1) = si;
526 
527  retval(0) = y;
528  }
529  else
530  error (errmsg);
531  }
532  else
533  error (errmsg);
534  }
535  }
536  else
537  {
538  if (isfloat)
539  {
540  FloatColumnVector b (args(0).float_vector_value ());
541  FloatColumnVector a (args(1).float_vector_value ());
542 
543  FloatNDArray x (args(2).float_array_value ());
544 
545  if (! error_state)
546  {
547  FloatNDArray si;
548 
549  if (nargin == 3 || args(3).is_empty ())
550  {
551  octave_idx_type a_len = a.length ();
552  octave_idx_type b_len = b.length ();
553 
554  octave_idx_type si_len = (a_len > b_len ? a_len : b_len) - 1;
555 
556  dim_vector si_dims = x.dims ();
557  for (int i = dim; i > 0; i--)
558  si_dims(i) = si_dims(i-1);
559  si_dims(0) = si_len;
560 
561  si.resize (si_dims, 0.0);
562  }
563  else
564  {
565  si = args(3).float_array_value ();
566 
567  if (si.is_vector () && x.is_vector ())
568  si = si.reshape (dim_vector (si.numel (), 1));
569  }
570 
571  if (! error_state)
572  {
573  FloatNDArray y (filter (b, a, x, si, dim));
574 
575  if (nargout == 2)
576  retval(1) = si;
577 
578  retval(0) = y;
579  }
580  else
581  error (errmsg);
582  }
583  else
584  error (errmsg);
585  }
586  else
587  {
588  ColumnVector b (args(0).vector_value ());
589  ColumnVector a (args(1).vector_value ());
590 
591  NDArray x (args(2).array_value ());
592 
593  if (! error_state)
594  {
595  NDArray si;
596 
597  if (nargin == 3 || args(3).is_empty ())
598  {
599  octave_idx_type a_len = a.length ();
600  octave_idx_type b_len = b.length ();
601 
602  octave_idx_type si_len = (a_len > b_len ? a_len : b_len) - 1;
603 
604  dim_vector si_dims = x.dims ();
605  for (int i = dim; i > 0; i--)
606  si_dims(i) = si_dims(i-1);
607  si_dims(0) = si_len;
608 
609  si.resize (si_dims, 0.0);
610  }
611  else
612  {
613  si = args(3).array_value ();
614 
615  if (si.is_vector () && x.is_vector ())
616  si = si.reshape (dim_vector (si.numel (), 1));
617  }
618 
619  if (! error_state)
620  {
621  NDArray y (filter (b, a, x, si, dim));
622 
623  if (nargout == 2)
624  retval(1) = si;
625 
626  retval(0) = y;
627  }
628  else
629  error (errmsg);
630  }
631  else
632  error (errmsg);
633  }
634  }
635 
636  return retval;
637 }
638 
639 template MArray<double>
640 filter (MArray<double>&, MArray<double>&, MArray<double>&,
641  MArray<double>&, int dim);
642 
643 template MArray<double>
644 filter (MArray<double>&, MArray<double>&, MArray<double>&, int dim);
645 
646 template MArray<Complex>
647 filter (MArray<Complex>&, MArray<Complex>&, MArray<Complex>&,
648  MArray<Complex>&, int dim);
649 
650 template MArray<Complex>
651 filter (MArray<Complex>&, MArray<Complex>&, MArray<Complex>&, int dim);
652 
653 template MArray<float>
654 filter (MArray<float>&, MArray<float>&, MArray<float>&,
655  MArray<float>&, int dim);
656 
657 template MArray<float>
658 filter (MArray<float>&, MArray<float>&, MArray<float>&, int dim);
659 
660 template MArray<FloatComplex>
661 filter (MArray<FloatComplex>&, MArray<FloatComplex>&, MArray<FloatComplex>&,
662  MArray<FloatComplex>&, int dim);
663 
664 template MArray<FloatComplex>
665 filter (MArray<FloatComplex>&, MArray<FloatComplex>&, MArray<FloatComplex>&,
666  int dim);
667 
668 /*
669 %!shared a, b, x, r
670 %!test
671 %! a = [1 1];
672 %! b = [1 1];
673 %! x = zeros (1,10); x(1) = 1;
674 %! assert (filter (b, [1], x ), [1 1 0 0 0 0 0 0 0 0]);
675 %! assert (filter (b, [1], x.'), [1 1 0 0 0 0 0 0 0 0].');
676 %! assert (filter (b.', [1], x ), [1 1 0 0 0 0 0 0 0 0] );
677 %! assert (filter (b.', [1], x.'), [1 1 0 0 0 0 0 0 0 0].');
678 %! assert (filter ([1], a, x ), [+1 -1 +1 -1 +1 -1 +1 -1 +1 -1] );
679 %! assert (filter ([1], a, x.'), [+1 -1 +1 -1 +1 -1 +1 -1 +1 -1].');
680 %! assert (filter ([1], a.', x ), [+1 -1 +1 -1 +1 -1 +1 -1 +1 -1] );
681 %! assert (filter ([1], a.', x.'), [+1 -1 +1 -1 +1 -1 +1 -1 +1 -1].');
682 %! assert (filter (b, a, x ), [1 0 0 0 0 0 0 0 0 0] );
683 %! assert (filter (b.', a, x ), [1 0 0 0 0 0 0 0 0 0] );
684 %! assert (filter (b, a.', x ), [1 0 0 0 0 0 0 0 0 0] );
685 %! assert (filter (b.', a, x ), [1 0 0 0 0 0 0 0 0 0] );
686 %! assert (filter (b, a, x.'), [1 0 0 0 0 0 0 0 0 0].');
687 %! assert (filter (b.', a, x.'), [1 0 0 0 0 0 0 0 0 0].');
688 %! assert (filter (b, a.', x.'), [1 0 0 0 0 0 0 0 0 0].');
689 %! assert (filter (b.', a, x.'), [1 0 0 0 0 0 0 0 0 0].');
690 
691 %!test
692 %! r = sqrt (1/2) * (1+i);
693 %! a = a*r;
694 %! b = b*r;
695 %! assert (filter (b, [1], x ), r*[1 1 0 0 0 0 0 0 0 0] );
696 %! assert (filter (b, [1], r*x ), r*r*[1 1 0 0 0 0 0 0 0 0] );
697 %! assert (filter (b, [1], x.' ), r*[1 1 0 0 0 0 0 0 0 0].' );
698 %! assert (filter (b, a, x ), [1 0 0 0 0 0 0 0 0 0] );
699 %! assert (filter (b, a, r*x ), r*[1 0 0 0 0 0 0 0 0 0] );
700 
701 %!shared a, b, x, y, so
702 %!test
703 %! a = [1,1];
704 %! b = [1,1];
705 %! x = zeros (1,10); x(1) = 1;
706 %! [y, so] = filter (b, [1], x, [-1]);
707 %! assert (y, [0 1 0 0 0 0 0 0 0 0]);
708 %! assert (so, 0);
709 
710 %!test
711 %! x = zeros (10,3); x(1,1) = -1; x(1,2) = 1;
712 %! y0 = zeros (10,3); y0(1:2,1) = -1; y0(1:2,2) = 1;
713 %! y = filter (b, [1], x);
714 %! assert (y, y0);
715 
716 %!test
717 %! a = [1,1];
718 %! b=[1,1];
719 %! x = zeros (4,4,2); x(1,1:4,1) = +1; x(1,1:4,2) = -1;
720 %! y0 = zeros (4,4,2); y0(1:2,1:4,1) = +1; y0(1:2,1:4,2) = -1;
721 %! y = filter (b, [1], x);
722 %! assert (y, y0);
723 
724 %!assert (filter (1, ones (10,1) / 10, []), [])
725 %!assert (filter (1, ones (10,1) / 10, zeros (0,10)), zeros (0,10))
726 %!assert (filter (1, ones (10,1) / 10, single (1:5)), repmat (single (10), 1, 5))
727 
728 %% Test using initial conditions
729 %!assert (filter ([1, 1, 1], [1, 1], [1 2], [1, 1]), [2 2])
730 %!assert (filter ([1, 1, 1], [1, 1], [1 2], [1, 1]'), [2 2])
731 %!assert (filter ([1, 3], [1], [1 2; 3 4; 5 6], [4, 5]), [5 7; 6 10; 14 18])
732 %!error (filter ([1, 3], [1], [1 2; 3 4; 5 6], [4, 5]'))
733 %!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])
734 
735 ## Test of DIM parameter
736 %!test
737 %! x = ones (2, 1, 3, 4);
738 %! x(1,1,:,:) = [1 2 3 4; 5 6 7 8; 9 10 11 12];
739 %! 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];
740 %! y0 = reshape (y0, size (x));
741 %! y = filter ([1 1 1], 1, x, [], 3);
742 %! assert (y, y0);
743 */
Generated on Mon Dec 30 2013 03:04:23 for GNU Octave by   doxygen 1.8.1.2