1 | function C = conv(A, B) |
---|
2 | % CONV Convolution and polynomial multiplication. |
---|
3 | % (Quaternion overloading of standard Matlab function.) |
---|
4 | |
---|
5 | % Copyright © 2006 Stephen J. Sangwine and Nicolas Le Bihan. |
---|
6 | % See the file : Copyright.m for further details. |
---|
7 | |
---|
8 | % This function operates in the same way as the standard Matlab function |
---|
9 | % apart from supporting both left and right coefficients. (Since quaternion |
---|
10 | % multiplication is not commutative, the general case requires both left |
---|
11 | % and right multiplication in the convolution product/summation.) To supply |
---|
12 | % left and right coefficients, use the calling profile conv({L,R},v) where |
---|
13 | % L and R are vectors of the same length and orientation. If the first |
---|
14 | % parameter is not a cell array, it is taken to be a left coefficient and |
---|
15 | % the right coefficient array is implicitly ones. |
---|
16 | |
---|
17 | error(nargchk(2, 2, nargin)), error(nargoutchk(0, 1, nargout)) |
---|
18 | |
---|
19 | if iscell(A) |
---|
20 | % The first parameter must be a cell array of two elements, and the |
---|
21 | % elements must have the same size (therefore orientation). |
---|
22 | |
---|
23 | if length(A) ~= 2 |
---|
24 | error('If a cell array, the first parameter must have two elements.') |
---|
25 | end |
---|
26 | |
---|
27 | L = A{1}; R = A{2}; % Separate the cell array into two vectors. |
---|
28 | |
---|
29 | if ~all(size(L) == size(R)) |
---|
30 | error('The elements of the cell array must have the same size.') |
---|
31 | end |
---|
32 | else |
---|
33 | L = A; |
---|
34 | R = ones(size(L)); % Supply the implicit right coefficients. |
---|
35 | end |
---|
36 | |
---|
37 | if ~isvector(L) || ~isvector(B) |
---|
38 | error('Parameters must be vectors.') |
---|
39 | end |
---|
40 | |
---|
41 | m = length(L); |
---|
42 | n = length(B); |
---|
43 | |
---|
44 | % Construct a result vector C matching the orientation of the longer of |
---|
45 | % L or B in order to return a result compatible with the standard Matlab |
---|
46 | % conv function. |
---|
47 | |
---|
48 | if m > n |
---|
49 | % w must have the same orientation as u. |
---|
50 | if size(L, 1) == 1 |
---|
51 | C = quaternion(zeros(1, m + n - 1)); |
---|
52 | else |
---|
53 | C = quaternion(zeros(m + n - 1, 1)); |
---|
54 | end |
---|
55 | end |
---|
56 | if m < n |
---|
57 | % w must have the same orientation as q. |
---|
58 | if size(B, 1) == 1 |
---|
59 | C = quaternion(zeros(1, m + n - 1)); |
---|
60 | else |
---|
61 | C = quaternion(zeros(m + n - 1, 1)); |
---|
62 | end |
---|
63 | end |
---|
64 | if m == n |
---|
65 | % The two parameters are of equal length. In this case Matlab's conv |
---|
66 | % function seems to return a column vector unless both parameters are |
---|
67 | % row vectors. |
---|
68 | if size(L, 1) == 1 && size(B, 1) == 1 |
---|
69 | C = quaternion(zeros(1, m + n - 1)); |
---|
70 | else |
---|
71 | C = quaternion(zeros(m + n - 1, 1)); |
---|
72 | end |
---|
73 | end |
---|
74 | |
---|
75 | % To ensure that we can multiply L .* B .* R we need them to have the same |
---|
76 | % orientation (row, column). If B doesn't match the orientation of L, |
---|
77 | % transpose it. |
---|
78 | |
---|
79 | if (size(L, 1) == 1 && size(B, 1) ~= 1) || (size(L, 2) == 1 && size(B, 2) ~= 1) |
---|
80 | B = B.'; |
---|
81 | end |
---|
82 | |
---|
83 | for k = 1:length(C) % For each element of the output array. |
---|
84 | |
---|
85 | % This code is vectorized: the two arrays are sliced to extract the |
---|
86 | % parts that must be multiplied element by element and summed to give |
---|
87 | % the current element of the output array. This is complicated by the |
---|
88 | % cases where one array overlaps the end of the other. |
---|
89 | |
---|
90 | j1 = max(1, k + 1 - n); % These two indices select a region of the L |
---|
91 | j2 = min(k, m); % and R arrays (the whole of, in the general |
---|
92 | % case). |
---|
93 | |
---|
94 | % We need to use indexing here, but since this is a class method, |
---|
95 | % normal array indexing doesn't work. See the comment in the file |
---|
96 | % svd.m where the same problem and solution occurs. The commented |
---|
97 | % statement below is what we would write if array indexing worked |
---|
98 | % here as per normal. |
---|
99 | % C(k) = sum(L(j1 : j2) .* B(k + 1 - j1 : k + 1 - j2) .* R(j1 : j2)); |
---|
100 | |
---|
101 | s = substruct('()',{k}); |
---|
102 | t = substruct('()',{j1:j2}); |
---|
103 | C = subsasgn(C, s, ... |
---|
104 | sum(subsref(L, t) .* ... |
---|
105 | flipud(fliplr(subsref(B, substruct('()', {k + 1 - j2 : k + 1 - j1})))) .* ... |
---|
106 | subsref(R, t))); |
---|
107 | end |
---|
108 | |
---|
109 | % Implementation note: the use of ones(size(u)) for v in the case where v |
---|
110 | % is omitted is a crude hack and it would be better not to multiply by the |
---|
111 | % right coefficients at all. Also if only two parameters are given it would |
---|
112 | % be better to assign the longer of the two to q, and the shorter to u or v |
---|
113 | % according as to whether it should be multiplied on the left or right. |
---|
114 | |
---|
115 | % The code also needs to be sorted out to tidy up the mess caused by trying |
---|
116 | % to return a result with the same orientation as Matlab's built-in in |
---|
117 | % conv. |
---|