Mercurial > hg > octave-lyh

comparison scripts/polynomial/conv.m @ 11082:4558aad4c41d

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
conv.m: handle "same" shape argument
author John W. Eaton <jwe@octave.org>
date Thu, 07 Oct 2010 03:09:07 -0400
parents f6e0404421f4
children 0f6c5efce96e
comparison
equal deleted inserted replaced
11081:f9284142a060 11082:4558aad4c41d
16 ## You should have received a copy of the GNU General Public License 16 ## You should have received a copy of the GNU General Public License
17 ## along with Octave; see the file COPYING. If not, see 17 ## along with Octave; see the file COPYING. If not, see
18 ## <http://www.gnu.org/licenses/>. 18 ## <http://www.gnu.org/licenses/>.
19 19
20 ## -*- texinfo -*- 20 ## -*- texinfo -*-
21 ## @deftypefn {Function File} {} conv (@var{a}, @var{b}) 21 ## @deftypefn {Function File} {} conv (@var{a}, @var{b}, @var{shape})
22 ## Convolve two vectors. 22 ## Convolve two vectors.
23 ## 23 ##
24 ## @code{y = conv (a, b)} returns a vector of length equal to 24 ## @code{y = conv (a, b)} returns a vector of length equal to
25 ## @code{length (a) + length (b) - 1}. 25 ## @code{length (a) + length (b) - 1}.
26 ## If @var{a} and @var{b} are polynomial coefficient vectors, @code{conv} 26 ## If @var{a} and @var{b} are polynomial coefficient vectors, @code{conv}
27 ## returns the coefficients of the product polynomial. 27 ## returns the coefficients of the product polynomial.
28 ##
29 ## The optional @var{shape} parameter may be
30 ##
31 ## @table @asis\n\
32 ## @item @var{shape} = "full"
33 ## Return the full convolution.
34 ## \n\
35 ## @item @var{shape} = "same"
36 ## Return central part of the convolution with the same size as @var{a}.
37 ## @end table
38 ##
39 ## @noindent
40 ## By default, @var{shape} is @samp{"full"}.
41 ##
28 ## @seealso{deconv, poly, roots, residue, polyval, polyderiv, polyint} 42 ## @seealso{deconv, poly, roots, residue, polyval, polyderiv, polyint}
29 ## @end deftypefn 43 ## @end deftypefn
30 44
31 ## Author: Tony Richardson <arichard@stark.cc.oh.us> 45 ## Author: Tony Richardson <arichard@stark.cc.oh.us>
32 ## Created: June 1994 46 ## Created: June 1994
33 ## Adapted-By: jwe 47 ## Adapted-By: jwe
34 48
35 function y = conv (a, b) 49 function y = conv (a, b, shape = "full")
36 50
37 if (nargin != 2) 51 if (nargin < 2 || nargin > 3)
38 print_usage (); 52 print_usage ();
39 endif 53 endif
40 54
41 if (! (isvector (a) && isvector (b))) 55 if (! (isvector (a) && isvector (b)))
42 error ("conv: both arguments must be vectors"); 56 error ("conv: both arguments must be vectors");
45 la = length (a); 59 la = length (a);
46 lb = length (b); 60 lb = length (b);
47 61
48 ly = la + lb - 1; 62 ly = la + lb - 1;
49 63
50 ## Use the shortest vector as the coefficent vector to filter. 64 if (ly == 0)
51 ## Preserve the row/column orientation of the longer input. 65 y = zeros (1, 0);
52 if (la <= lb) 66 else
53 if (ly > lb) 67 ## Use the shortest vector as the coefficent vector to filter.
54 if (size (b, 1) <= size (b, 2)) 68 ## Preserve the row/column orientation of the longer input.
55 x = [b, (zeros (1, ly - lb))]; 69 if (la <= lb)
70 if (ly > lb)
71 if (size (b, 1) <= size (b, 2))
72 x = [b, (zeros (1, ly - lb))];
73 else
74 x = [b; (zeros (ly - lb, 1))];
75 endif
56 else 76 else
57 x = [b; (zeros (ly - lb, 1))]; 77 x = b;
58 endif 78 endif
79 y = filter (a, 1, x);
59 else 80 else
60 x = b; 81 if (ly > la)
82 if (size (a, 1) <= size (a, 2))
83 x = [a, (zeros (1, ly - la))];
84 else
85 x = [a; (zeros (ly - la, 1))];
86 endif
87 else
88 x = a;
89 endif
90 y = filter (b, 1, x);
61 endif 91 endif
62 y = filter (a, 1, x); 92 if (strcmp (shape, "same"))
63 else 93 idx = ceil ((ly - la) / 2);
64 if (ly > la) 94 y = y(idx+1:idx+la);
65 if (size (a, 1) <= size (a, 2))
66 x = [a, (zeros (1, ly - la))];
67 else
68 x = [a; (zeros (ly - la, 1))];
69 endif
70 else
71 x = a;
72 endif 95 endif
73 y = filter (b, 1, x);
74 endif 96 endif
75 97
76 endfunction 98 endfunction
77 99
78 %!test 100 %!test
97 %! b = 1:3; 119 %! b = 1:3;
98 %! assert (size(conv(a,b)), [1, numel(a)+numel(b)-1]) 120 %! assert (size(conv(a,b)), [1, numel(a)+numel(b)-1])
99 %! assert (size(conv(b,a)), [1, numel(a)+numel(b)-1]) 121 %! assert (size(conv(b,a)), [1, numel(a)+numel(b)-1])
100 122
101 %!test 123 %!test
124 %! a = 1:10;
125 %! b = 1:3;
126 %! assert (size(conv(a,b)), [1, numel(a)+numel(b)-1])
127 %! assert (size(conv(b,a)), [1, numel(a)+numel(b)-1])
128
102 %! a = (1:10).'; 129 %! a = (1:10).';
103 %! b = 1:3; 130 %! b = 1:3;
104 %! assert (size(conv(a,b)), [numel(a)+numel(b)-1, 1]) 131 %! assert (size(conv(a,b)), [numel(a)+numel(b)-1, 1])
105 %! assert (size(conv(b,a)), [numel(a)+numel(b)-1, 1]) 132 %! assert (size(conv(b,a)), [numel(a)+numel(b)-1, 1])
106 133
107 %!test 134 %!test
108 %! a = 1:10; 135 %! a = 1:10;
136 %! b = 1:3;
137 %! assert (conv(a,b,'same'), [4, 10, 16, 22, 28, 34, 40, 46, 52, 47])
138 %! assert (conv(b,a,'same'), [28, 34, 40])
139
140 %!test
141 %! a = 1:10;
109 %! b = (1:3).'; 142 %! b = (1:3).';
110 %! assert (size(conv(a,b)), [1, numel(a)+numel(b)-1]) 143 %! assert (size(conv(a,b)), [1, numel(a)+numel(b)-1])
111 %! assert (size(conv(b,a)), [1, numel(a)+numel(b)-1]) 144 %! assert (size(conv(b,a)), [1, numel(a)+numel(b)-1])