[37] | 1 | function Y = qfft(X, A, L) |
---|
| 2 | % QFFT Discrete quaternion Fourier transform. |
---|
| 3 | % |
---|
| 4 | % This function calculates the fast quaternion Fourier transform of |
---|
| 5 | % (columns of) X. A specifies the transform axis (that is the direction |
---|
| 6 | % in 3-space of the vector part of the hypercomplex exponential). It must |
---|
| 7 | % be a pure quaternion (real or complex) but it need not have unit modulus. |
---|
| 8 | % L specifies whether the quaternion exponential is on the left ('L') or |
---|
| 9 | % right ('R'). |
---|
| 10 | % |
---|
| 11 | % See also: IQFFT, FFT, IFFT. |
---|
| 12 | |
---|
| 13 | % Copyright © 2005 Stephen J. Sangwine and Nicolas Le Bihan. |
---|
| 14 | % See the file : Copyright.m for further details. |
---|
| 15 | |
---|
| 16 | % Reference: |
---|
| 17 | % |
---|
| 18 | % Salem Said, Nicolas Le Bihan and Stephen J. Sangwine, |
---|
| 19 | % Fast complexified quaternion Fourier transform, |
---|
| 20 | % arXiv:math.NA/0603578, 24 March 2006. Available at http://www.arxiv.org/. |
---|
| 21 | |
---|
| 22 | error(nargchk(3, 3, nargin)), error(nargoutchk(0, 1, nargout)) |
---|
| 23 | |
---|
| 24 | if size(A) ~= [1, 1] |
---|
| 25 | error('The transform axis cannot be a matrix or vector.'); |
---|
| 26 | end |
---|
| 27 | |
---|
| 28 | if ~isa(A, 'quaternion') | ~ispure(A) |
---|
| 29 | error('The transform axis must be a pure quaternion.') |
---|
| 30 | end |
---|
| 31 | |
---|
| 32 | if L ~= 'L' & L ~= 'R' |
---|
| 33 | error('L must have the value ''L'' or ''R''.'); |
---|
| 34 | end |
---|
| 35 | |
---|
| 36 | S = 1; |
---|
| 37 | if L == 'R' |
---|
| 38 | S = -1; % S is a sign bit used (in effect) to conjugate one of the complex |
---|
| 39 | % components below when the exponential is on the right. In fact, |
---|
| 40 | % instead of conjugating the exponential (which would require an |
---|
| 41 | % inverse fft (ifft), we conjugate the complex component before and |
---|
| 42 | % after the transformation. This achieves the same effect because |
---|
| 43 | % the inverse transform may always be computed by taking the |
---|
| 44 | % conjugate before and after the transformation (this is a |
---|
| 45 | % standard DFT trick). |
---|
| 46 | end |
---|
| 47 | |
---|
| 48 | A = unit(A); % Ensure that A is a unit (pure) quaternion. |
---|
| 49 | B = orthonormal_basis(A); |
---|
| 50 | |
---|
| 51 | % Compute the transform. This is done by changing the basis of all the |
---|
| 52 | % elements of the matrix or vector X to the transform axis (which may be |
---|
| 53 | % real or complex). The quaternion elements can then be regrouped as |
---|
| 54 | % complex values and their FFTs computed with the standard Matlab function. |
---|
| 55 | % We then have to regroup the complex results into quaternion form and |
---|
| 56 | % invert the change of basis. |
---|
| 57 | |
---|
| 58 | % Note that we compute here a real quaternion FFT (two complex FFTs) if |
---|
| 59 | % both X and the axis are real, otherwise we compute a complex quaternion |
---|
| 60 | % FFT (four complex FFTs). The decision is made after the change of basis, |
---|
| 61 | % because by this point, if X is still real, a real quaternion FFT is |
---|
| 62 | % needed. It is theoretically possible to use the complex code for both |
---|
| 63 | % cases, but in the real case this would result in two unnecessary FFTs |
---|
| 64 | % being computed, which is a heavy cost. |
---|
| 65 | |
---|
| 66 | X = change_basis(X, B); |
---|
| 67 | |
---|
| 68 | if isreal(X) |
---|
| 69 | |
---|
| 70 | % Compute the two complex FFTs using the standard Matlab complex FFT |
---|
| 71 | % function. |
---|
| 72 | |
---|
| 73 | C1 = fft(complex(scalar(X), x(X))); |
---|
| 74 | C2 = fft(complex( y(X), S .* z(X))); |
---|
| 75 | |
---|
| 76 | % Compose a real quaternion result from the two complex results. |
---|
| 77 | |
---|
| 78 | Y = quaternion(real(C1), imag(C1), real(C2), S .* imag(C2)); |
---|
| 79 | else |
---|
| 80 | |
---|
| 81 | % Compute the four complex FFTs using the standard Matlab complex FFT |
---|
| 82 | % function. |
---|
| 83 | |
---|
| 84 | C1 = fft(complex(real(scalar(X)), real(x(X)))); |
---|
| 85 | C2 = fft(complex(imag(scalar(X)), imag(x(X)))); |
---|
| 86 | C3 = fft(complex(real( y(X)), S .* real(z(X)))); |
---|
| 87 | C4 = fft(complex(imag( y(X)), S .* imag(z(X)))); |
---|
| 88 | |
---|
| 89 | % Notice here that conjugation of the C3 and C4 components |
---|
| 90 | % (multiplication by S) is implemented by negating the complex number |
---|
| 91 | % formed from their imaginary parts. |
---|
| 92 | |
---|
| 93 | Y = quaternion(complex(real(C1), real(C2)), ... |
---|
| 94 | complex(imag(C1), imag(C2)), ... |
---|
| 95 | complex(real(C3), real(C4)), ... |
---|
| 96 | S .* complex(imag(C3), imag(C4))); |
---|
| 97 | end |
---|
| 98 | |
---|
| 99 | Y = change_basis(Y, B.'); % Change back to the original basis. |
---|