[37] | 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. |
---|