view dcti.m @ 17:0da5eb32e49d tip

minor changes

function c = dcti(x,type)
% DCTI discrete cosine transform of types I-IV
%
% Example 
% C = DCTI(X,TYPE) returns the dct of the signal x of the specified type
%
% Input
% - x: vector or matrix containing the signal to be analized. If matrix, the
% function acts columnwise.
% - type: (optional, default = 'II') string containing the dct type. Valid
% types are 'I','II','III' and 'IV'
%
% Output
% - c: vector or matrix contaning the dct coefficients
%
% References
% 1) S. Mallat, A wavelet tour of signal processing
%
% See also IDCTI, DCT, FFT, IDCT, IFFT

% Author(s): Daniele Barchiesi
% Copyright 2011-2011

%% Check inputs & defaults
error(nargchk(1,2,nargin,'struct'));
if nargin==1, type = 'II';  end

%% Compute dct
[N M] = size(x);
switch upper(type)
    case 'I'
        % Explicit calculation of DCT-I. Replace with fast algorithm (todo)
        dctMat = zeros(N);
        dctMat(:,1) = 1/sqrt(N)*ones(N,1);
        for k=1:N-1
            dctMat(:,k+1) = sqrt(2/N)*cos(k*pi/N*(0.5+(0:N-1)));
        end
        c = dctMat'*x;
    case 'II' % Fast calculation of DCT II. Uses method employed in dct.m
        w = (exp(-1i*(0:N-1)*pi/(2*N))/sqrt(2*N)).';
        w(1) = w(1)/sqrt(2);
        if rem(N,2) || ~isreal(x), % odd case
            % Form intermediate even-symmetric matrix
            y = zeros(2*N,M);
            y(1:N,:) = x;
            y(N+1:2*N,:) = flipud(x);
            
            % Compute the FFT and keep the appropriate portion:
            yy = fft(y);
            yy = yy(1:N,:);
        else % even case
            % Re-order the elements of the columns of x
            y = [x(1:2:N,:); x(N:-2:2,:)];
            yy = fft(y);
            w = 2*w;  % Double the weights for even-length case
        end
        % Multiply FFT by weights:
        c = w(:,ones(1,M)).*yy;
    case 'III' 
        % Explicit calculation of DCT III. Replace with fast algorithm (todo)
        dctMat = zeros(N);
        for k=0:N-1
            dctMat(:,k+1) = [sqrt(2)/2; cos(pi/N*(k+1/2)*(1:N-1))'];
        end
        dctMat = normcol(dctMat);
        c = dctMat'*x;
    case 'IV' % Fast calculation of DCT IV.
        y = zeros(8*N,M);
        y(2:2:2*N,:) = x;
        yy = fft(y,[],1);
        c = sqrt(2/N)*real(yy(2:2:2*N,:));
    otherwise
        error('Unsupported dct type');
end

%% Postprocess
if isreal(x), c = real(c); end

%EOF