danieleb@2: function x = idcti(c,type)
danieleb@14: % IDCTI inverse discrete cosine transform of types I-IV
danieleb@2: %
danieleb@14: % Example 
danieleb@14: % X = IDCTI(C,TYPE) returns the idct of the signal c of the specified type
danieleb@2: %
danieleb@14: % Input
danieleb@14: % - c: vector or matrix containing the dct coefficients. If matrix, the
danieleb@2: % function acts columnwise.
danieleb@14: % - type: (optional, default = 'II') string containing the idct type. Valid
danieleb@2: % types are 'I','II','III' and 'IV'
danieleb@2: %
danieleb@2: % OUTPUT
danieleb@14: % - x: vector or matrix contaning the inverse dct coefficients
danieleb@2: %
danieleb@14: % References
danieleb@14: % 1) S. Mallat, A wavelet tour of signal processing
danieleb@2: %
danieleb@14: % See also DCTI, DCT, FFT, IDCT, IFFT
danieleb@2: 
danieleb@4: % Author(s): Daniele Barchiesi
danieleb@14: % Copyright 2011-2011
danieleb@2: 
danieleb@2: %% Check inputs & defaults
danieleb@2: error(nargchk(1,2,nargin,'struct'));
danieleb@2: if nargin==1, type = 'II';  end
danieleb@2: 
danieleb@2: %% Compute idct
danieleb@2: [N M] = size(c);
danieleb@2: switch upper(type)
danieleb@9:     case 'I'
danieleb@13:         % Explicit calculation of IDCT-I. Replace with fast algorithm (todo)
danieleb@2:         dctMat = zeros(N);
danieleb@9:         dctMat(:,1) = 1/sqrt(N)*ones(N,1);
danieleb@9:         for k=1:N-1
danieleb@9:             dctMat(:,k+1) = sqrt(2/N)*cos(k*pi/N*(0.5+(0:N-1)));
danieleb@2:         end
danieleb@2:         x = dctMat*c;
danieleb@2:     case 'II' % Fast calculation of DCT II. Uses method employed in idct.M
danieleb@2:         % Compute wieghts
danieleb@2:         w = sqrt(2*N) * exp(1i*(0:N-1)*pi/(2*N)).';
danieleb@2:         if rem(N,2) || ~isreal(c), % odd case
danieleb@2:             % Form intermediate even-symmetric matrix.
danieleb@2:             w(1) = w(1) * sqrt(2);
danieleb@2:             W = w(:,ones(1,M));
danieleb@2:             yy = zeros(2*N,M);
danieleb@2:             yy(1:N,:) = W.*c;
danieleb@2:             yy(N+2:2*N,:) = -1i*W(2:N,:).*flipud(c(2:N,:));
danieleb@2:             
danieleb@2:             y = ifft(yy);
danieleb@2:             
danieleb@2:             % Extract inverse DCT
danieleb@2:             x = y(1:N,:);
danieleb@2:             
danieleb@2:         else % even case
danieleb@2:             w(1) = w(1)/sqrt(2);
danieleb@2:             W = w(:,ones(1,M));
danieleb@2:             yy = W.*c;
danieleb@2:             
danieleb@2:             % Compute x tilde using equation (5.93) in Jain
danieleb@2:             y = ifft(yy);
danieleb@2:             
danieleb@2:             % Re-order elements of each column according to equations (5.93) and
danieleb@2:             % (5.94) in Jain
danieleb@2:             x = zeros(N,M);
danieleb@2:             x(1:2:N,:) = y(1:N/2,:);
danieleb@2:             x(2:2:N,:) = y(N:-1:N/2+1,:);
danieleb@2:         end
danieleb@13:     case 'III' % Explicit calculation of IDCT III. Replace with fast algorithm (todo)
danieleb@2:         dctMat = zeros(N);
danieleb@2:         for k=0:N-1
danieleb@2:             dctMat(:,k+1) = [sqrt(2)/2; cos(pi/N*(k+1/2)*(1:N-1))'];
danieleb@2:         end
danieleb@2:         dctMat = normcol(dctMat);
danieleb@2:         x = dctMat*c;
danieleb@2:     case 'IV' % Fast calculation of IDCT IV.
danieleb@2:         y = zeros(8*N,M);
danieleb@2:         y(2:2:2*N,:) = c;
danieleb@2:         yy = fft(y,[],1);
danieleb@2:         x = sqrt(2/N)*real(yy(2:2:2*N,:));
danieleb@2:     otherwise
danieleb@2:         error('Unsupported idct type');
danieleb@2: end
danieleb@2: 
danieleb@2: %% Postprocess
danieleb@2: if isreal(c), x = real(x); end
danieleb@2: 
danieleb@2: %EOF