www.gusucode.com > 图像FFT_DCT源码程序 > 图像FFT_DCT源码程序/spectral_saliency_matlab/qtfm/@quaternion/qdct2.m

    function Y = qdct2(X, A, L)
  % QDCT2 calcultes the Quaternion DCT (see qfft2)
  %
  % @author: B. Schauerte
  % @date: 2011
  
  if nargin < 2, A = dft_axis(isreal(X)); end
  if nargin < 3, L = 'L'; end
  
  assert(isreal(X)); % don't want to implement the stuff for complex X
  
  %error(nargchk(3, 3, nargin)), error(nargoutchk(0, 1, nargout))

  if ~isscalar(A)
      error('The transform axis cannot be a matrix or vector.');
  end

  if ~isa(A, 'quaternion') || ~isempty(A.w)
      error('The transform axis must be a pure quaternion.')
  end

  if L ~= 'L' && L ~= 'R'
      error('L must have the value ''L'' or ''R''.');
  end

  S = 1;
  if L == 'R'
      S = -1;  % S is a sign bit used (in effect) to conjugate one of the complex
               % components below when the exponential is on the right.  In fact,
               % instead of conjugating the exponential (which would require an
               % inverse fft (ifft), we conjugate the complex component before and 
               % after the transformation. This achieves the same effect because
               % the inverse transform may always be computed by taking the
               % conjugate before and after the transformation (this is a
               % standard DFT trick).
  end

  A = unit(A); % Ensure that A is a unit (pure) quaternion.
  B = orthonormal_basis(A);

  X = change_basis(X, B);

  if isreal(X)
    % Compute the two complex FFTs using the standard Matlab complex FFT
    % function.

    C1 = dct2(complex(scalar(X),      x(X)));
    C2 = dct2(complex(     y(X), S .* z(X)));

    % Compose a real quaternion result from the two complex results.

    Y = quaternion(real(C1), imag(C1), real(C2), S .* imag(C2));
  else
    error('qdct2 is not implemented for complex X')
  end
  
  Y = change_basis(Y, B.'); % Change back to the original basis.