annotate src/DLD-FUNCTIONS/filter.cc @ 5823:080c08b192d8

[project @ 2006-05-19 05:32:17 by jwe]
author jwe
date Fri, 19 May 2006 05:32:19 +0000
parents 2a16423e4aa0
children 2533c8c41aad
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
2928
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
1 /*
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
2
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
3 Copyright (C) 1996, 1997 John W. Eaton
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
4
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
5 This file is part of Octave.
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
6
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
7 Octave is free software; you can redistribute it and/or modify it
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
8 under the terms of the GNU General Public License as published by the
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
9 Free Software Foundation; either version 2, or (at your option) any
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
10 later version.
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
11
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
12 Octave is distributed in the hope that it will be useful, but WITHOUT
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
13 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
14 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
15 for more details.
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
16
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
17 You should have received a copy of the GNU General Public License
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
18 along with Octave; see the file COPYING. If not, write to the Free
5307
4c8a2e4e0717 [project @ 2005-04-26 19:24:27 by jwe]
jwe
parents: 5275
diff changeset
19 Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
4c8a2e4e0717 [project @ 2005-04-26 19:24:27 by jwe]
jwe
parents: 5275
diff changeset
20 02110-1301, USA.
2928
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
21
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
22 */
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
23
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
24 // Based on Tony Richardson's filter.m.
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
25 //
5428
2a16423e4aa0 [project @ 2005-08-23 18:38:27 by jwe]
jwe
parents: 5307
diff changeset
26 // Originally translated to C++ by KH (Kurt.Hornik@wu-wien.ac.at)
2928
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
27 // with help from Fritz Leisch and Andreas Weingessel on Oct 20, 1994.
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
28 //
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
29 // Rewritten to use templates to handle both real and complex cases by
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
30 // jwe, Wed Nov 1 19:15:29 1995.
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
31
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
32 #ifdef HAVE_CONFIG_H
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
33 #include <config.h>
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
34 #endif
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
35
4153
6b96ce9f5743 [project @ 2002-11-06 20:38:49 by jwe]
jwe
parents: 3419
diff changeset
36 #include "quit.h"
6b96ce9f5743 [project @ 2002-11-06 20:38:49 by jwe]
jwe
parents: 3419
diff changeset
37
2928
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
38 #include "defun-dld.h"
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
39 #include "error.h"
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
40 #include "oct-obj.h"
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
41
3109
fcb2f5063e7f [project @ 1997-11-20 08:52:42 by jwe]
jwe
parents: 3107
diff changeset
42 #if !defined (CXX_NEW_FRIEND_TEMPLATE_DECL)
4844
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
43 extern MArrayN<double>
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
44 filter (MArray<double>&, MArray<double>&, MArrayN<double>&, int dim);
2928
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
45
4844
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
46 extern MArrayN<Complex>
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
47 filter (MArray<Complex>&, MArray<Complex>&, MArrayN<Complex>&, int dim);
3107
a8dcfbf87ea3 [project @ 1997-11-19 22:34:59 by jwe]
jwe
parents: 3097
diff changeset
48 #endif
2928
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
49
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
50 template <class T>
4844
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
51 MArrayN<T>
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
52 filter (MArray<T>& b, MArray<T>& a, MArrayN<T>& x, MArrayN<T>& si,
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
53 int dim = 0)
2928
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
54 {
4844
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
55 MArrayN<T> y;
2928
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
56
5275
23b37da9fd5b [project @ 2005-04-08 16:07:35 by jwe]
jwe
parents: 5114
diff changeset
57 octave_idx_type a_len = a.length ();
23b37da9fd5b [project @ 2005-04-08 16:07:35 by jwe]
jwe
parents: 5114
diff changeset
58 octave_idx_type b_len = b.length ();
2928
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
59
5275
23b37da9fd5b [project @ 2005-04-08 16:07:35 by jwe]
jwe
parents: 5114
diff changeset
60 octave_idx_type ab_len = a_len > b_len ? a_len : b_len;
2928
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
61
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
62 b.resize (ab_len, 0.0);
4844
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
63 if (a_len > 1)
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
64 a.resize (ab_len, 0.0);
2928
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
65
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
66 T norm = a (0);
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
67
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
68 if (norm == 0.0)
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
69 {
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
70 error ("filter: the first element of a must be non-zero");
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
71 return y;
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
72 }
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
73
4844
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
74 dim_vector x_dims = x.dims ();
5056
9379d8716407 [project @ 2004-10-22 14:58:05 by jwe]
jwe
parents: 5055
diff changeset
75 if (dim < 0 || dim > x_dims.length ())
4844
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
76 {
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
77 error ("filter: filtering over invalid dimension");
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
78 return y;
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
79 }
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
80
5275
23b37da9fd5b [project @ 2005-04-08 16:07:35 by jwe]
jwe
parents: 5114
diff changeset
81 octave_idx_type x_len = x_dims(dim);
4844
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
82
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
83 dim_vector si_dims = si.dims ();
5275
23b37da9fd5b [project @ 2005-04-08 16:07:35 by jwe]
jwe
parents: 5114
diff changeset
84 octave_idx_type si_len = si_dims(0);
4844
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
85
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
86 if (si_len != ab_len - 1)
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
87 {
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
88 error ("filter: first dimension of si must be of length max (length (a), length (b)) - 1");
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
89 return y;
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
90 }
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
91
5059
38dd53c8935e [project @ 2004-11-02 02:42:25 by jwe]
jwe
parents: 5056
diff changeset
92 if (si_dims.length () != x_dims.length ())
4844
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
93 {
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
94 error ("filter: dimensionality of si and x must agree");
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
95 return y;
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
96 }
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
97
5275
23b37da9fd5b [project @ 2005-04-08 16:07:35 by jwe]
jwe
parents: 5114
diff changeset
98 octave_idx_type si_dim = 0;
23b37da9fd5b [project @ 2005-04-08 16:07:35 by jwe]
jwe
parents: 5114
diff changeset
99 for (octave_idx_type i = 0; i < x_dims.length (); i++)
4844
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
100 {
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
101 if (i == dim)
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
102 continue;
5055
51a4406317e9 [project @ 2004-10-22 14:46:52 by jwe]
jwe
parents: 4901
diff changeset
103
51a4406317e9 [project @ 2004-10-22 14:46:52 by jwe]
jwe
parents: 4901
diff changeset
104 if (x_dims(i) == 1)
51a4406317e9 [project @ 2004-10-22 14:46:52 by jwe]
jwe
parents: 4901
diff changeset
105 continue;
51a4406317e9 [project @ 2004-10-22 14:46:52 by jwe]
jwe
parents: 4901
diff changeset
106
5059
38dd53c8935e [project @ 2004-11-02 02:42:25 by jwe]
jwe
parents: 5056
diff changeset
107 if (si_dims(++si_dim) != x_dims(i))
4844
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
108 {
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
109 error ("filter: dimensionality of si and x must agree");
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
110 return y;
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
111 }
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
112 }
2928
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
113
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
114 if (norm != 1.0)
4844
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
115 {
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
116 a = a / norm;
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
117 b = b / norm;
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
118 }
2928
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
119
5056
9379d8716407 [project @ 2004-10-22 14:58:05 by jwe]
jwe
parents: 5055
diff changeset
120 if (a_len <= 1 && si_len <= 0)
4844
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
121 return b(0) * x;
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
122
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
123 y.resize (x_dims, 0.0);
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
124
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
125 int x_stride = 1;
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
126 for (int i = 0; i < dim; i++)
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
127 x_stride *= x_dims(i);
2928
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
128
5275
23b37da9fd5b [project @ 2005-04-08 16:07:35 by jwe]
jwe
parents: 5114
diff changeset
129 octave_idx_type x_num = x_dims.numel () / x_len;
23b37da9fd5b [project @ 2005-04-08 16:07:35 by jwe]
jwe
parents: 5114
diff changeset
130 for (octave_idx_type num = 0; num < x_num; num++)
4844
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
131 {
5275
23b37da9fd5b [project @ 2005-04-08 16:07:35 by jwe]
jwe
parents: 5114
diff changeset
132 octave_idx_type x_offset;
4844
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
133 if (x_stride == 1)
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
134 x_offset = num * x_len;
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
135 else
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
136 {
5275
23b37da9fd5b [project @ 2005-04-08 16:07:35 by jwe]
jwe
parents: 5114
diff changeset
137 octave_idx_type x_offset2 = 0;
4844
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
138 x_offset = num;
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
139 while (x_offset >= x_stride)
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
140 {
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
141 x_offset -= x_stride;
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
142 x_offset2++;
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
143 }
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
144 x_offset += x_offset2 * x_stride * x_len;
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
145 }
5275
23b37da9fd5b [project @ 2005-04-08 16:07:35 by jwe]
jwe
parents: 5114
diff changeset
146 octave_idx_type si_offset = num * si_len;
2928
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
147
4844
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
148 if (a_len > 1)
2928
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
149 {
5114
cda8c0a823c5 [project @ 2005-01-21 17:02:21 by jwe]
jwe
parents: 5059
diff changeset
150 T *py = y.fortran_vec ();
cda8c0a823c5 [project @ 2005-01-21 17:02:21 by jwe]
jwe
parents: 5059
diff changeset
151 T *psi = si.fortran_vec ();
cda8c0a823c5 [project @ 2005-01-21 17:02:21 by jwe]
jwe
parents: 5059
diff changeset
152
cda8c0a823c5 [project @ 2005-01-21 17:02:21 by jwe]
jwe
parents: 5059
diff changeset
153 const T *pa = a.data ();
cda8c0a823c5 [project @ 2005-01-21 17:02:21 by jwe]
jwe
parents: 5059
diff changeset
154 const T *pb = b.data ();
cda8c0a823c5 [project @ 2005-01-21 17:02:21 by jwe]
jwe
parents: 5059
diff changeset
155 const T *px = x.data ();
cda8c0a823c5 [project @ 2005-01-21 17:02:21 by jwe]
jwe
parents: 5059
diff changeset
156
cda8c0a823c5 [project @ 2005-01-21 17:02:21 by jwe]
jwe
parents: 5059
diff changeset
157 psi += si_offset;
cda8c0a823c5 [project @ 2005-01-21 17:02:21 by jwe]
jwe
parents: 5059
diff changeset
158
5275
23b37da9fd5b [project @ 2005-04-08 16:07:35 by jwe]
jwe
parents: 5114
diff changeset
159 for (octave_idx_type i = 0, idx = x_offset; i < x_len; i++, idx += x_stride)
2928
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
160 {
5114
cda8c0a823c5 [project @ 2005-01-21 17:02:21 by jwe]
jwe
parents: 5059
diff changeset
161 py[idx] = psi[0] + pb[0] * px[idx];
cda8c0a823c5 [project @ 2005-01-21 17:02:21 by jwe]
jwe
parents: 5059
diff changeset
162
cda8c0a823c5 [project @ 2005-01-21 17:02:21 by jwe]
jwe
parents: 5059
diff changeset
163 if (si_len > 0)
cda8c0a823c5 [project @ 2005-01-21 17:02:21 by jwe]
jwe
parents: 5059
diff changeset
164 {
5275
23b37da9fd5b [project @ 2005-04-08 16:07:35 by jwe]
jwe
parents: 5114
diff changeset
165 for (octave_idx_type j = 0; j < si_len - 1; j++)
5114
cda8c0a823c5 [project @ 2005-01-21 17:02:21 by jwe]
jwe
parents: 5059
diff changeset
166 {
cda8c0a823c5 [project @ 2005-01-21 17:02:21 by jwe]
jwe
parents: 5059
diff changeset
167 OCTAVE_QUIT;
cda8c0a823c5 [project @ 2005-01-21 17:02:21 by jwe]
jwe
parents: 5059
diff changeset
168
cda8c0a823c5 [project @ 2005-01-21 17:02:21 by jwe]
jwe
parents: 5059
diff changeset
169 psi[j] = psi[j+1] - pa[j+1] * py[idx] + pb[j+1] * px[idx];
cda8c0a823c5 [project @ 2005-01-21 17:02:21 by jwe]
jwe
parents: 5059
diff changeset
170 }
cda8c0a823c5 [project @ 2005-01-21 17:02:21 by jwe]
jwe
parents: 5059
diff changeset
171
cda8c0a823c5 [project @ 2005-01-21 17:02:21 by jwe]
jwe
parents: 5059
diff changeset
172 psi[si_len-1] = pb[si_len] * px[idx] - pa[si_len] * py[idx];
cda8c0a823c5 [project @ 2005-01-21 17:02:21 by jwe]
jwe
parents: 5059
diff changeset
173 }
cda8c0a823c5 [project @ 2005-01-21 17:02:21 by jwe]
jwe
parents: 5059
diff changeset
174 else
cda8c0a823c5 [project @ 2005-01-21 17:02:21 by jwe]
jwe
parents: 5059
diff changeset
175 {
cda8c0a823c5 [project @ 2005-01-21 17:02:21 by jwe]
jwe
parents: 5059
diff changeset
176 OCTAVE_QUIT;
cda8c0a823c5 [project @ 2005-01-21 17:02:21 by jwe]
jwe
parents: 5059
diff changeset
177
cda8c0a823c5 [project @ 2005-01-21 17:02:21 by jwe]
jwe
parents: 5059
diff changeset
178 psi[0] = pb[si_len] * px[idx] - pa[si_len] * py[idx];
cda8c0a823c5 [project @ 2005-01-21 17:02:21 by jwe]
jwe
parents: 5059
diff changeset
179 }
cda8c0a823c5 [project @ 2005-01-21 17:02:21 by jwe]
jwe
parents: 5059
diff changeset
180 }
cda8c0a823c5 [project @ 2005-01-21 17:02:21 by jwe]
jwe
parents: 5059
diff changeset
181 }
cda8c0a823c5 [project @ 2005-01-21 17:02:21 by jwe]
jwe
parents: 5059
diff changeset
182 else if (si_len > 0)
cda8c0a823c5 [project @ 2005-01-21 17:02:21 by jwe]
jwe
parents: 5059
diff changeset
183 {
cda8c0a823c5 [project @ 2005-01-21 17:02:21 by jwe]
jwe
parents: 5059
diff changeset
184 T *py = y.fortran_vec ();
cda8c0a823c5 [project @ 2005-01-21 17:02:21 by jwe]
jwe
parents: 5059
diff changeset
185 T *psi = si.fortran_vec ();
cda8c0a823c5 [project @ 2005-01-21 17:02:21 by jwe]
jwe
parents: 5059
diff changeset
186
cda8c0a823c5 [project @ 2005-01-21 17:02:21 by jwe]
jwe
parents: 5059
diff changeset
187 const T *pb = b.data ();
cda8c0a823c5 [project @ 2005-01-21 17:02:21 by jwe]
jwe
parents: 5059
diff changeset
188 const T *px = x.data ();
cda8c0a823c5 [project @ 2005-01-21 17:02:21 by jwe]
jwe
parents: 5059
diff changeset
189
cda8c0a823c5 [project @ 2005-01-21 17:02:21 by jwe]
jwe
parents: 5059
diff changeset
190 psi += si_offset;
cda8c0a823c5 [project @ 2005-01-21 17:02:21 by jwe]
jwe
parents: 5059
diff changeset
191
5275
23b37da9fd5b [project @ 2005-04-08 16:07:35 by jwe]
jwe
parents: 5114
diff changeset
192 for (octave_idx_type i = 0, idx = x_offset; i < x_len; i++, idx += x_stride)
5114
cda8c0a823c5 [project @ 2005-01-21 17:02:21 by jwe]
jwe
parents: 5059
diff changeset
193 {
cda8c0a823c5 [project @ 2005-01-21 17:02:21 by jwe]
jwe
parents: 5059
diff changeset
194 py[idx] = psi[0] + pb[0] * px[idx];
4844
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
195
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
196 if (si_len > 1)
4153
6b96ce9f5743 [project @ 2002-11-06 20:38:49 by jwe]
jwe
parents: 3419
diff changeset
197 {
5275
23b37da9fd5b [project @ 2005-04-08 16:07:35 by jwe]
jwe
parents: 5114
diff changeset
198 for (octave_idx_type j = 0; j < si_len - 1; j++)
4844
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
199 {
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
200 OCTAVE_QUIT;
4153
6b96ce9f5743 [project @ 2002-11-06 20:38:49 by jwe]
jwe
parents: 3419
diff changeset
201
5114
cda8c0a823c5 [project @ 2005-01-21 17:02:21 by jwe]
jwe
parents: 5059
diff changeset
202 psi[j] = psi[j+1] + pb[j+1] * px[idx];
4844
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
203 }
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
204
5114
cda8c0a823c5 [project @ 2005-01-21 17:02:21 by jwe]
jwe
parents: 5059
diff changeset
205 psi[si_len-1] = pb[si_len] * px[idx];
4153
6b96ce9f5743 [project @ 2002-11-06 20:38:49 by jwe]
jwe
parents: 3419
diff changeset
206 }
4844
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
207 else
5114
cda8c0a823c5 [project @ 2005-01-21 17:02:21 by jwe]
jwe
parents: 5059
diff changeset
208 {
cda8c0a823c5 [project @ 2005-01-21 17:02:21 by jwe]
jwe
parents: 5059
diff changeset
209 OCTAVE_QUIT;
4844
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
210
5114
cda8c0a823c5 [project @ 2005-01-21 17:02:21 by jwe]
jwe
parents: 5059
diff changeset
211 psi[0] = pb[1] * px[idx];
4844
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
212 }
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
213 }
2928
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
214 }
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
215 }
4844
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
216
2928
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
217 return y;
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
218 }
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
219
3109
fcb2f5063e7f [project @ 1997-11-20 08:52:42 by jwe]
jwe
parents: 3107
diff changeset
220 #if !defined (CXX_NEW_FRIEND_TEMPLATE_DECL)
4844
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
221 extern MArrayN<double>
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
222 filter (MArray<double>&, MArray<double>&, MArrayN<double>&,
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
223 MArrayN<double>&, int dim);
2928
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
224
4844
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
225 extern MArrayN<Complex>
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
226 filter (MArray<Complex>&, MArray<Complex>&, MArrayN<Complex>&,
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
227 MArrayN<Complex>&, int dim);
3107
a8dcfbf87ea3 [project @ 1997-11-19 22:34:59 by jwe]
jwe
parents: 3097
diff changeset
228 #endif
2928
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
229
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
230 template <class T>
4844
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
231 MArrayN<T>
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
232 filter (MArray<T>& b, MArray<T>& a, MArrayN<T>& x, int dim = -1)
2928
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
233 {
5059
38dd53c8935e [project @ 2004-11-02 02:42:25 by jwe]
jwe
parents: 5056
diff changeset
234 dim_vector x_dims = x.dims();
4844
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
235
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
236 if (dim < 0)
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
237 {
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
238 // Find first non-singleton dimension
5059
38dd53c8935e [project @ 2004-11-02 02:42:25 by jwe]
jwe
parents: 5056
diff changeset
239 while (dim < x_dims.length () && x_dims(dim) <= 1)
4844
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
240 dim++;
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
241
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
242 // All dimensions singleton, pick first dimension
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
243 if (dim == x_dims.length ())
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
244 dim = 0;
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
245 }
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
246 else
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
247 if (dim < 0 || dim > x_dims.length ())
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
248 {
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
249 error ("filter: filtering over invalid dimension");
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
250 return MArrayN<T> ();
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
251 }
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
252
5275
23b37da9fd5b [project @ 2005-04-08 16:07:35 by jwe]
jwe
parents: 5114
diff changeset
253 octave_idx_type a_len = a.length ();
23b37da9fd5b [project @ 2005-04-08 16:07:35 by jwe]
jwe
parents: 5114
diff changeset
254 octave_idx_type b_len = b.length ();
2928
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
255
5275
23b37da9fd5b [project @ 2005-04-08 16:07:35 by jwe]
jwe
parents: 5114
diff changeset
256 octave_idx_type si_len = (a_len > b_len ? a_len : b_len) - 1;
4844
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
257 dim_vector si_dims = x.dims ();
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
258 for (int i = dim; i > 0; i--)
5059
38dd53c8935e [project @ 2004-11-02 02:42:25 by jwe]
jwe
parents: 5056
diff changeset
259 si_dims(i) = si_dims(i-1);
38dd53c8935e [project @ 2004-11-02 02:42:25 by jwe]
jwe
parents: 5056
diff changeset
260 si_dims(0) = si_len;
4844
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
261
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
262 MArrayN<T> si (si_dims, T (0.0));
2928
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
263
4844
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
264 return filter (b, a, x, si, dim);
2928
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
265 }
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
266
3097
d12c312d1cfb [project @ 1997-11-14 08:03:47 by jwe]
jwe
parents: 3088
diff changeset
267 DEFUN_DLD (filter, args, nargout,
3367
0748b03c3510 [project @ 1999-11-20 14:52:38 by jwe]
jwe
parents: 3109
diff changeset
268 "-*- texinfo -*-\n\
0748b03c3510 [project @ 1999-11-20 14:52:38 by jwe]
jwe
parents: 3109
diff changeset
269 @deftypefn {Loadable Function} {y =} filter (@var{b}, @var{a}, @var{x})\n\
0748b03c3510 [project @ 1999-11-20 14:52:38 by jwe]
jwe
parents: 3109
diff changeset
270 @deftypefnx {Loadable Function} {[@var{y}, @var{sf}] =} filter (@var{b}, @var{a}, @var{x}, @var{si})\n\
4844
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
271 @deftypefnx {Loadable Function} {[@var{y}, @var{sf}] =} filter (@var{b}, @var{a}, @var{x}, [], @var{dim})\n\
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
272 @deftypefnx {Loadable Function} {[@var{y}, @var{sf}] =} filter (@var{b}, @var{a}, @var{x}, @var{si}, @var{dim})\n\
3367
0748b03c3510 [project @ 1999-11-20 14:52:38 by jwe]
jwe
parents: 3109
diff changeset
273 Return the solution to the following linear, time-invariant difference\n\
0748b03c3510 [project @ 1999-11-20 14:52:38 by jwe]
jwe
parents: 3109
diff changeset
274 equation:\n\
0748b03c3510 [project @ 1999-11-20 14:52:38 by jwe]
jwe
parents: 3109
diff changeset
275 @iftex\n\
0748b03c3510 [project @ 1999-11-20 14:52:38 by jwe]
jwe
parents: 3109
diff changeset
276 @tex\n\
0748b03c3510 [project @ 1999-11-20 14:52:38 by jwe]
jwe
parents: 3109
diff changeset
277 $$\n\
0748b03c3510 [project @ 1999-11-20 14:52:38 by jwe]
jwe
parents: 3109
diff changeset
278 \\sum_{k=0}^N a_{k+1} y_{n-k} = \\sum_{k=0}^M b_{k+1} x_{n-k}, \\qquad\n\
0748b03c3510 [project @ 1999-11-20 14:52:38 by jwe]
jwe
parents: 3109
diff changeset
279 1 \\le n \\le P\n\
0748b03c3510 [project @ 1999-11-20 14:52:38 by jwe]
jwe
parents: 3109
diff changeset
280 $$\n\
0748b03c3510 [project @ 1999-11-20 14:52:38 by jwe]
jwe
parents: 3109
diff changeset
281 @end tex\n\
0748b03c3510 [project @ 1999-11-20 14:52:38 by jwe]
jwe
parents: 3109
diff changeset
282 @end iftex\n\
0748b03c3510 [project @ 1999-11-20 14:52:38 by jwe]
jwe
parents: 3109
diff changeset
283 @ifinfo\n\
2928
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
284 \n\
3367
0748b03c3510 [project @ 1999-11-20 14:52:38 by jwe]
jwe
parents: 3109
diff changeset
285 @smallexample\n\
0748b03c3510 [project @ 1999-11-20 14:52:38 by jwe]
jwe
parents: 3109
diff changeset
286 N M\n\
0748b03c3510 [project @ 1999-11-20 14:52:38 by jwe]
jwe
parents: 3109
diff changeset
287 SUM a(k+1) y(n-k) = SUM b(k+1) x(n-k) for 1<=n<=length(x)\n\
0748b03c3510 [project @ 1999-11-20 14:52:38 by jwe]
jwe
parents: 3109
diff changeset
288 k=0 k=0\n\
0748b03c3510 [project @ 1999-11-20 14:52:38 by jwe]
jwe
parents: 3109
diff changeset
289 @end smallexample\n\
0748b03c3510 [project @ 1999-11-20 14:52:38 by jwe]
jwe
parents: 3109
diff changeset
290 @end ifinfo\n\
2928
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
291 \n\
3367
0748b03c3510 [project @ 1999-11-20 14:52:38 by jwe]
jwe
parents: 3109
diff changeset
292 @noindent\n\
0748b03c3510 [project @ 1999-11-20 14:52:38 by jwe]
jwe
parents: 3109
diff changeset
293 where\n\
0748b03c3510 [project @ 1999-11-20 14:52:38 by jwe]
jwe
parents: 3109
diff changeset
294 @ifinfo\n\
0748b03c3510 [project @ 1999-11-20 14:52:38 by jwe]
jwe
parents: 3109
diff changeset
295 N=length(a)-1 and M=length(b)-1.\n\
0748b03c3510 [project @ 1999-11-20 14:52:38 by jwe]
jwe
parents: 3109
diff changeset
296 @end ifinfo\n\
0748b03c3510 [project @ 1999-11-20 14:52:38 by jwe]
jwe
parents: 3109
diff changeset
297 @iftex\n\
0748b03c3510 [project @ 1999-11-20 14:52:38 by jwe]
jwe
parents: 3109
diff changeset
298 @tex\n\
0748b03c3510 [project @ 1999-11-20 14:52:38 by jwe]
jwe
parents: 3109
diff changeset
299 $a \\in \\Re^{N-1}$, $b \\in \\Re^{M-1}$, and $x \\in \\Re^P$.\n\
0748b03c3510 [project @ 1999-11-20 14:52:38 by jwe]
jwe
parents: 3109
diff changeset
300 @end tex\n\
0748b03c3510 [project @ 1999-11-20 14:52:38 by jwe]
jwe
parents: 3109
diff changeset
301 @end iftex\n\
4844
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
302 over the first non-singleton dimension of @var{x} or over @var{dim} if\n\
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
303 supplied. An equivalent form of this equation is:\n\
3367
0748b03c3510 [project @ 1999-11-20 14:52:38 by jwe]
jwe
parents: 3109
diff changeset
304 @iftex\n\
0748b03c3510 [project @ 1999-11-20 14:52:38 by jwe]
jwe
parents: 3109
diff changeset
305 @tex\n\
0748b03c3510 [project @ 1999-11-20 14:52:38 by jwe]
jwe
parents: 3109
diff changeset
306 $$\n\
0748b03c3510 [project @ 1999-11-20 14:52:38 by jwe]
jwe
parents: 3109
diff changeset
307 y_n = -\\sum_{k=1}^N c_{k+1} y_{n-k} + \\sum_{k=0}^M d_{k+1} x_{n-k}, \\qquad\n\
0748b03c3510 [project @ 1999-11-20 14:52:38 by jwe]
jwe
parents: 3109
diff changeset
308 1 \\le n \\le P\n\
0748b03c3510 [project @ 1999-11-20 14:52:38 by jwe]
jwe
parents: 3109
diff changeset
309 $$\n\
0748b03c3510 [project @ 1999-11-20 14:52:38 by jwe]
jwe
parents: 3109
diff changeset
310 @end tex\n\
0748b03c3510 [project @ 1999-11-20 14:52:38 by jwe]
jwe
parents: 3109
diff changeset
311 @end iftex\n\
0748b03c3510 [project @ 1999-11-20 14:52:38 by jwe]
jwe
parents: 3109
diff changeset
312 @ifinfo\n\
2928
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
313 \n\
3367
0748b03c3510 [project @ 1999-11-20 14:52:38 by jwe]
jwe
parents: 3109
diff changeset
314 @smallexample\n\
0748b03c3510 [project @ 1999-11-20 14:52:38 by jwe]
jwe
parents: 3109
diff changeset
315 N M\n\
0748b03c3510 [project @ 1999-11-20 14:52:38 by jwe]
jwe
parents: 3109
diff changeset
316 y(n) = - SUM c(k+1) y(n-k) + SUM d(k+1) x(n-k) for 1<=n<=length(x)\n\
0748b03c3510 [project @ 1999-11-20 14:52:38 by jwe]
jwe
parents: 3109
diff changeset
317 k=1 k=0\n\
0748b03c3510 [project @ 1999-11-20 14:52:38 by jwe]
jwe
parents: 3109
diff changeset
318 @end smallexample\n\
0748b03c3510 [project @ 1999-11-20 14:52:38 by jwe]
jwe
parents: 3109
diff changeset
319 @end ifinfo\n\
0748b03c3510 [project @ 1999-11-20 14:52:38 by jwe]
jwe
parents: 3109
diff changeset
320 \n\
0748b03c3510 [project @ 1999-11-20 14:52:38 by jwe]
jwe
parents: 3109
diff changeset
321 @noindent\n\
0748b03c3510 [project @ 1999-11-20 14:52:38 by jwe]
jwe
parents: 3109
diff changeset
322 where\n\
0748b03c3510 [project @ 1999-11-20 14:52:38 by jwe]
jwe
parents: 3109
diff changeset
323 @ifinfo\n\
0748b03c3510 [project @ 1999-11-20 14:52:38 by jwe]
jwe
parents: 3109
diff changeset
324 c = a/a(1) and d = b/a(1).\n\
0748b03c3510 [project @ 1999-11-20 14:52:38 by jwe]
jwe
parents: 3109
diff changeset
325 @end ifinfo\n\
0748b03c3510 [project @ 1999-11-20 14:52:38 by jwe]
jwe
parents: 3109
diff changeset
326 @iftex\n\
0748b03c3510 [project @ 1999-11-20 14:52:38 by jwe]
jwe
parents: 3109
diff changeset
327 @tex\n\
0748b03c3510 [project @ 1999-11-20 14:52:38 by jwe]
jwe
parents: 3109
diff changeset
328 $c = a/a_1$ and $d = b/a_1$.\n\
0748b03c3510 [project @ 1999-11-20 14:52:38 by jwe]
jwe
parents: 3109
diff changeset
329 @end tex\n\
0748b03c3510 [project @ 1999-11-20 14:52:38 by jwe]
jwe
parents: 3109
diff changeset
330 @end iftex\n\
0748b03c3510 [project @ 1999-11-20 14:52:38 by jwe]
jwe
parents: 3109
diff changeset
331 \n\
0748b03c3510 [project @ 1999-11-20 14:52:38 by jwe]
jwe
parents: 3109
diff changeset
332 If the fourth argument @var{si} is provided, it is taken as the\n\
0748b03c3510 [project @ 1999-11-20 14:52:38 by jwe]
jwe
parents: 3109
diff changeset
333 initial state of the system and the final state is returned as\n\
0748b03c3510 [project @ 1999-11-20 14:52:38 by jwe]
jwe
parents: 3109
diff changeset
334 @var{sf}. The state vector is a column vector whose length is\n\
0748b03c3510 [project @ 1999-11-20 14:52:38 by jwe]
jwe
parents: 3109
diff changeset
335 equal to the length of the longest coefficient vector minus one.\n\
0748b03c3510 [project @ 1999-11-20 14:52:38 by jwe]
jwe
parents: 3109
diff changeset
336 If @var{si} is not supplied, the initial state vector is set to all\n\
2928
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
337 zeros.\n\
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
338 \n\
3367
0748b03c3510 [project @ 1999-11-20 14:52:38 by jwe]
jwe
parents: 3109
diff changeset
339 In terms of the z-transform, y is the result of passing the discrete-\n\
0748b03c3510 [project @ 1999-11-20 14:52:38 by jwe]
jwe
parents: 3109
diff changeset
340 time signal x through a system characterized by the following rational\n\
0748b03c3510 [project @ 1999-11-20 14:52:38 by jwe]
jwe
parents: 3109
diff changeset
341 system function:\n\
0748b03c3510 [project @ 1999-11-20 14:52:38 by jwe]
jwe
parents: 3109
diff changeset
342 @iftex\n\
0748b03c3510 [project @ 1999-11-20 14:52:38 by jwe]
jwe
parents: 3109
diff changeset
343 @tex\n\
0748b03c3510 [project @ 1999-11-20 14:52:38 by jwe]
jwe
parents: 3109
diff changeset
344 $$\n\
0748b03c3510 [project @ 1999-11-20 14:52:38 by jwe]
jwe
parents: 3109
diff changeset
345 H(z) = {\\displaystyle\\sum_{k=0}^M d_{k+1} z^{-k}\n\
0748b03c3510 [project @ 1999-11-20 14:52:38 by jwe]
jwe
parents: 3109
diff changeset
346 \\over 1 + \\displaystyle\\sum_{k+1}^N c_{k+1} z^{-k}}\n\
0748b03c3510 [project @ 1999-11-20 14:52:38 by jwe]
jwe
parents: 3109
diff changeset
347 $$\n\
0748b03c3510 [project @ 1999-11-20 14:52:38 by jwe]
jwe
parents: 3109
diff changeset
348 @end tex\n\
0748b03c3510 [project @ 1999-11-20 14:52:38 by jwe]
jwe
parents: 3109
diff changeset
349 @end iftex\n\
0748b03c3510 [project @ 1999-11-20 14:52:38 by jwe]
jwe
parents: 3109
diff changeset
350 @ifinfo\n\
0748b03c3510 [project @ 1999-11-20 14:52:38 by jwe]
jwe
parents: 3109
diff changeset
351 \n\
0748b03c3510 [project @ 1999-11-20 14:52:38 by jwe]
jwe
parents: 3109
diff changeset
352 @example\n\
0748b03c3510 [project @ 1999-11-20 14:52:38 by jwe]
jwe
parents: 3109
diff changeset
353 M\n\
0748b03c3510 [project @ 1999-11-20 14:52:38 by jwe]
jwe
parents: 3109
diff changeset
354 SUM d(k+1) z^(-k)\n\
0748b03c3510 [project @ 1999-11-20 14:52:38 by jwe]
jwe
parents: 3109
diff changeset
355 k=0\n\
0748b03c3510 [project @ 1999-11-20 14:52:38 by jwe]
jwe
parents: 3109
diff changeset
356 H(z) = ----------------------\n\
0748b03c3510 [project @ 1999-11-20 14:52:38 by jwe]
jwe
parents: 3109
diff changeset
357 N\n\
0748b03c3510 [project @ 1999-11-20 14:52:38 by jwe]
jwe
parents: 3109
diff changeset
358 1 + SUM c(k+1) z(-k)\n\
0748b03c3510 [project @ 1999-11-20 14:52:38 by jwe]
jwe
parents: 3109
diff changeset
359 k=1\n\
0748b03c3510 [project @ 1999-11-20 14:52:38 by jwe]
jwe
parents: 3109
diff changeset
360 @end example\n\
0748b03c3510 [project @ 1999-11-20 14:52:38 by jwe]
jwe
parents: 3109
diff changeset
361 @end ifinfo\n\
0748b03c3510 [project @ 1999-11-20 14:52:38 by jwe]
jwe
parents: 3109
diff changeset
362 @end deftypefn")
2928
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
363 {
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
364 octave_value_list retval;
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
365
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
366 int nargin = args.length ();
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
367
4844
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
368 if (nargin < 3 || nargin > 5)
2928
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
369 {
5823
080c08b192d8 [project @ 2006-05-19 05:32:17 by jwe]
jwe
parents: 5428
diff changeset
370 print_usage ();
2928
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
371 return retval;
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
372 }
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
373
4844
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
374 const char *errmsg = "filter: arguments a and b must be vectors";
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
375
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
376 int dim;
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
377 dim_vector x_dims = args(2).dims ();
2928
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
378
4844
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
379 if (nargin == 5)
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
380 {
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
381 dim = args(4).nint_value() - 1;
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
382 if (dim < 0 || dim >= x_dims.length ())
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
383 {
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
384 error ("filter: filtering over invalid dimension");
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
385 return retval;
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
386 }
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
387 }
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
388 else
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
389 {
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
390 // Find first non-singleton dimension
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
391 dim = 0;
5059
38dd53c8935e [project @ 2004-11-02 02:42:25 by jwe]
jwe
parents: 5056
diff changeset
392 while (dim < x_dims.length () && x_dims(dim) <= 1)
4844
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
393 dim++;
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
394
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
395 // All dimensions singleton, pick first dimension
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
396 if (dim == x_dims.length ())
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
397 dim = 0;
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
398 }
2928
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
399
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
400 if (args(0).is_complex_type ()
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
401 || args(1).is_complex_type ()
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
402 || args(2).is_complex_type ()
4844
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
403 || (nargin >= 4 && args(3).is_complex_type ()))
2928
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
404 {
3419
e71b3d1dd327 [project @ 2000-01-12 05:23:34 by jwe]
jwe
parents: 3418
diff changeset
405 ComplexColumnVector b (args(0).complex_vector_value ());
e71b3d1dd327 [project @ 2000-01-12 05:23:34 by jwe]
jwe
parents: 3418
diff changeset
406 ComplexColumnVector a (args(1).complex_vector_value ());
4844
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
407
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
408 ComplexNDArray x (args(2).complex_array_value ());
2928
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
409
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
410 if (! error_state)
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
411 {
4844
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
412 ComplexNDArray si;
3088
f33738510eb2 [project @ 1997-10-02 06:38:36 by jwe]
jwe
parents: 3014
diff changeset
413
4844
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
414 if (nargin == 3 || args(3).is_empty ())
2928
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
415 {
5275
23b37da9fd5b [project @ 2005-04-08 16:07:35 by jwe]
jwe
parents: 5114
diff changeset
416 octave_idx_type a_len = a.length ();
23b37da9fd5b [project @ 2005-04-08 16:07:35 by jwe]
jwe
parents: 5114
diff changeset
417 octave_idx_type b_len = b.length ();
3088
f33738510eb2 [project @ 1997-10-02 06:38:36 by jwe]
jwe
parents: 3014
diff changeset
418
5275
23b37da9fd5b [project @ 2005-04-08 16:07:35 by jwe]
jwe
parents: 5114
diff changeset
419 octave_idx_type si_len = (a_len > b_len ? a_len : b_len) - 1;
3088
f33738510eb2 [project @ 1997-10-02 06:38:36 by jwe]
jwe
parents: 3014
diff changeset
420
4844
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
421 dim_vector si_dims = x.dims ();
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
422 for (int i = dim; i > 0; i--)
5059
38dd53c8935e [project @ 2004-11-02 02:42:25 by jwe]
jwe
parents: 5056
diff changeset
423 si_dims(i) = si_dims(i-1);
38dd53c8935e [project @ 2004-11-02 02:42:25 by jwe]
jwe
parents: 5056
diff changeset
424 si_dims(0) = si_len;
4844
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
425
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
426 si.resize (si_dims, 0.0);
3088
f33738510eb2 [project @ 1997-10-02 06:38:36 by jwe]
jwe
parents: 3014
diff changeset
427 }
f33738510eb2 [project @ 1997-10-02 06:38:36 by jwe]
jwe
parents: 3014
diff changeset
428 else
4844
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
429 {
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
430 dim_vector si_dims = args (3).dims ();
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
431 bool si_is_vector = true;
5059
38dd53c8935e [project @ 2004-11-02 02:42:25 by jwe]
jwe
parents: 5056
diff changeset
432 for (int i = 0; i < si_dims.length (); i++)
38dd53c8935e [project @ 2004-11-02 02:42:25 by jwe]
jwe
parents: 5056
diff changeset
433 if (si_dims(i) != 1 && si_dims(i) < si_dims.numel ())
4844
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
434 {
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
435 si_is_vector = false;
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
436 break;
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
437 }
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
438
5059
38dd53c8935e [project @ 2004-11-02 02:42:25 by jwe]
jwe
parents: 5056
diff changeset
439 si = args(3).complex_array_value ();
38dd53c8935e [project @ 2004-11-02 02:42:25 by jwe]
jwe
parents: 5056
diff changeset
440
4844
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
441 if (si_is_vector)
5059
38dd53c8935e [project @ 2004-11-02 02:42:25 by jwe]
jwe
parents: 5056
diff changeset
442 si = si.reshape (dim_vector (si.numel (), 1));
4844
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
443 }
3088
f33738510eb2 [project @ 1997-10-02 06:38:36 by jwe]
jwe
parents: 3014
diff changeset
444
f33738510eb2 [project @ 1997-10-02 06:38:36 by jwe]
jwe
parents: 3014
diff changeset
445 if (! error_state)
f33738510eb2 [project @ 1997-10-02 06:38:36 by jwe]
jwe
parents: 3014
diff changeset
446 {
4844
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
447 ComplexNDArray y (filter (b, a, x, si, dim));
3088
f33738510eb2 [project @ 1997-10-02 06:38:36 by jwe]
jwe
parents: 3014
diff changeset
448
f33738510eb2 [project @ 1997-10-02 06:38:36 by jwe]
jwe
parents: 3014
diff changeset
449 if (nargout == 2)
4844
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
450 retval(1) = si;
2928
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
451
4844
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
452 retval(0) = y;
2928
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
453 }
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
454 else
3088
f33738510eb2 [project @ 1997-10-02 06:38:36 by jwe]
jwe
parents: 3014
diff changeset
455 error (errmsg);
2928
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
456 }
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
457 else
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
458 error (errmsg);
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
459 }
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
460 else
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
461 {
3419
e71b3d1dd327 [project @ 2000-01-12 05:23:34 by jwe]
jwe
parents: 3418
diff changeset
462 ColumnVector b (args(0).vector_value ());
e71b3d1dd327 [project @ 2000-01-12 05:23:34 by jwe]
jwe
parents: 3418
diff changeset
463 ColumnVector a (args(1).vector_value ());
4844
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
464
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
465 NDArray x (args(2).array_value ());
2928
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
466
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
467 if (! error_state)
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
468 {
4844
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
469 NDArray si;
3088
f33738510eb2 [project @ 1997-10-02 06:38:36 by jwe]
jwe
parents: 3014
diff changeset
470
4844
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
471 if (nargin == 3 || args(3).is_empty ())
2928
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
472 {
5275
23b37da9fd5b [project @ 2005-04-08 16:07:35 by jwe]
jwe
parents: 5114
diff changeset
473 octave_idx_type a_len = a.length ();
23b37da9fd5b [project @ 2005-04-08 16:07:35 by jwe]
jwe
parents: 5114
diff changeset
474 octave_idx_type b_len = b.length ();
3088
f33738510eb2 [project @ 1997-10-02 06:38:36 by jwe]
jwe
parents: 3014
diff changeset
475
5275
23b37da9fd5b [project @ 2005-04-08 16:07:35 by jwe]
jwe
parents: 5114
diff changeset
476 octave_idx_type si_len = (a_len > b_len ? a_len : b_len) - 1;
3088
f33738510eb2 [project @ 1997-10-02 06:38:36 by jwe]
jwe
parents: 3014
diff changeset
477
4844
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
478 dim_vector si_dims = x.dims ();
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
479 for (int i = dim; i > 0; i--)
5059
38dd53c8935e [project @ 2004-11-02 02:42:25 by jwe]
jwe
parents: 5056
diff changeset
480 si_dims(i) = si_dims(i-1);
38dd53c8935e [project @ 2004-11-02 02:42:25 by jwe]
jwe
parents: 5056
diff changeset
481 si_dims(0) = si_len;
4844
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
482
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
483 si.resize (si_dims, 0.0);
3088
f33738510eb2 [project @ 1997-10-02 06:38:36 by jwe]
jwe
parents: 3014
diff changeset
484 }
f33738510eb2 [project @ 1997-10-02 06:38:36 by jwe]
jwe
parents: 3014
diff changeset
485 else
4844
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
486 {
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
487 dim_vector si_dims = args (3).dims ();
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
488 bool si_is_vector = true;
5059
38dd53c8935e [project @ 2004-11-02 02:42:25 by jwe]
jwe
parents: 5056
diff changeset
489 for (int i = 0; i < si_dims.length (); i++)
38dd53c8935e [project @ 2004-11-02 02:42:25 by jwe]
jwe
parents: 5056
diff changeset
490 if (si_dims(i) != 1 && si_dims(i) < si_dims.numel ())
4844
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
491 {
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
492 si_is_vector = false;
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
493 break;
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
494 }
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
495
5059
38dd53c8935e [project @ 2004-11-02 02:42:25 by jwe]
jwe
parents: 5056
diff changeset
496 si = args(3).array_value ();
38dd53c8935e [project @ 2004-11-02 02:42:25 by jwe]
jwe
parents: 5056
diff changeset
497
4844
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
498 if (si_is_vector)
5059
38dd53c8935e [project @ 2004-11-02 02:42:25 by jwe]
jwe
parents: 5056
diff changeset
499 si = si.reshape (dim_vector (si.numel (), 1));
4844
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
500 }
3088
f33738510eb2 [project @ 1997-10-02 06:38:36 by jwe]
jwe
parents: 3014
diff changeset
501
f33738510eb2 [project @ 1997-10-02 06:38:36 by jwe]
jwe
parents: 3014
diff changeset
502 if (! error_state)
f33738510eb2 [project @ 1997-10-02 06:38:36 by jwe]
jwe
parents: 3014
diff changeset
503 {
4844
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
504 NDArray y (filter (b, a, x, si, dim));
3088
f33738510eb2 [project @ 1997-10-02 06:38:36 by jwe]
jwe
parents: 3014
diff changeset
505
f33738510eb2 [project @ 1997-10-02 06:38:36 by jwe]
jwe
parents: 3014
diff changeset
506 if (nargout == 2)
4844
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
507 retval(1) = si;
2928
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
508
4844
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
509 retval(0) = y;
2928
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
510 }
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
511 else
3088
f33738510eb2 [project @ 1997-10-02 06:38:36 by jwe]
jwe
parents: 3014
diff changeset
512 error (errmsg);
2928
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
513 }
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
514 else
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
515 error (errmsg);
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
516 }
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
517
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
518 return retval;
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
519 }
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
520
4844
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
521 template MArrayN<double>
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
522 filter (MArray<double>&, MArray<double>&, MArrayN<double>&,
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
523 MArrayN<double>&, int dim);
2928
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
524
4844
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
525 template MArrayN<double>
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
526 filter (MArray<double>&, MArray<double>&, MArrayN<double>&, int dim);
2928
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
527
4844
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
528 template MArrayN<Complex>
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
529 filter (MArray<Complex>&, MArray<Complex>&, MArrayN<Complex>&,
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
530 MArrayN<Complex>&, int dim);
2928
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
531
4844
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
532 template MArrayN<Complex>
9f7ef92b50b0 [project @ 2004-04-02 17:26:53 by jwe]
jwe
parents: 4153
diff changeset
533 filter (MArray<Complex>&, MArray<Complex>&, MArrayN<Complex>&, int dim);
2928
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
534
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
535 /*
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
536 ;;; Local Variables: ***
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
537 ;;; mode: C++ ***
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
538 ;;; End: ***
295f037b4b3e [project @ 1997-05-05 05:32:33 by jwe]
jwe
parents:
diff changeset
539 */