1 | function Y = qfft2(X, A, L) |
---|
2 | % QFFT2 Discrete 2D quaternion Fourier transform. |
---|
3 | % |
---|
4 | % This function calculates the fast quaternion 2D 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: IQFFT2, FFT2, IFFT2. |
---|
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 = fft2(complex(scalar(X), x(X))); |
---|
74 | C2 = fft2(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 = fft2(complex(real(scalar(X)), real(x(X)))); |
---|
85 | C2 = fft2(complex(imag(scalar(X)), imag(x(X)))); |
---|
86 | C3 = fft2(complex(real( y(X)), S .* real(z(X)))); |
---|
87 | C4 = fft2(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. |
---|