view idcti.m @ 17:0da5eb32e49d tip

minor changes

function x = idcti(c,type)
% IDCTI inverse discrete cosine transform of types I-IV
%
% Example 
% X = IDCTI(C,TYPE) returns the idct of the signal c of the specified type
%
% Input
% - c: vector or matrix containing the dct coefficients. If matrix, the
% function acts columnwise.
% - type: (optional, default = 'II') string containing the idct type. Valid
% types are 'I','II','III' and 'IV'
%
% OUTPUT
% - x: vector or matrix contaning the inverse dct coefficients
%
% References
% 1) S. Mallat, A wavelet tour of signal processing
%
% See also DCTI, 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 idct
[N M] = size(c);
switch upper(type)
    case 'I'
        % Explicit calculation of IDCT-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
        x = dctMat*c;
    case 'II' % Fast calculation of DCT II. Uses method employed in idct.M
        % Compute wieghts
        w = sqrt(2*N) * exp(1i*(0:N-1)*pi/(2*N)).';
        if rem(N,2) || ~isreal(c), % odd case
            % Form intermediate even-symmetric matrix.
            w(1) = w(1) * sqrt(2);
            W = w(:,ones(1,M));
            yy = zeros(2*N,M);
            yy(1:N,:) = W.*c;
            yy(N+2:2*N,:) = -1i*W(2:N,:).*flipud(c(2:N,:));
            
            y = ifft(yy);
            
            % Extract inverse DCT
            x = y(1:N,:);
            
        else % even case
            w(1) = w(1)/sqrt(2);
            W = w(:,ones(1,M));
            yy = W.*c;
            
            % Compute x tilde using equation (5.93) in Jain
            y = ifft(yy);
            
            % Re-order elements of each column according to equations (5.93) and
            % (5.94) in Jain
            x = zeros(N,M);
            x(1:2:N,:) = y(1:N/2,:);
            x(2:2:N,:) = y(N:-1:N/2+1,:);
        end
    case 'III' % Explicit calculation of IDCT 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);
        x = dctMat*c;
    case 'IV' % Fast calculation of IDCT IV.
        y = zeros(8*N,M);
        y(2:2:2*N,:) = c;
        yy = fft(y,[],1);
        x = sqrt(2/N)*real(yy(2:2:2*N,:));
    otherwise
        error('Unsupported idct type');
end

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

%EOF