danieleb@1: function c = dcti(x,type)
danieleb@14: % DCTI discrete cosine transform of types I-IV
danieleb@1: %
danieleb@14: % Example 
danieleb@14: % C = DCTI(X,TYPE) returns the dct of the signal x of the specified type
danieleb@1: %
danieleb@14: % Input
danieleb@14: % - x: vector or matrix containing the signal to be analized. If matrix, the
danieleb@1: % function acts columnwise.
danieleb@14: % - type: (optional, default = 'II') string containing the dct type. Valid
danieleb@1: % types are 'I','II','III' and 'IV'
danieleb@1: %
danieleb@14: % Output
danieleb@14: % - c: vector or matrix contaning the dct coefficients
danieleb@1: %
danieleb@14: % References
danieleb@14: % 1) S. Mallat, A wavelet tour of signal processing
danieleb@1: %
danieleb@14: % See also IDCTI, DCT, FFT, IDCT, IFFT
danieleb@1: 
danieleb@4: % Author(s): Daniele Barchiesi
danieleb@14: % Copyright 2011-2011
danieleb@1: 
danieleb@1: %% Check inputs & defaults
danieleb@1: error(nargchk(1,2,nargin,'struct'));
danieleb@1: if nargin==1, type = 'II';  end
danieleb@1: 
danieleb@1: %% Compute dct
danieleb@1: [N M] = size(x);
danieleb@1: switch upper(type)
danieleb@9:     case 'I'
danieleb@13:         % Explicit calculation of DCT-I. Replace with fast algorithm (todo)
danieleb@1:         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@1:         end
danieleb@1:         c = dctMat'*x;
danieleb@1:     case 'II' % Fast calculation of DCT II. Uses method employed in dct.m
danieleb@1:         w = (exp(-1i*(0:N-1)*pi/(2*N))/sqrt(2*N)).';
danieleb@1:         w(1) = w(1)/sqrt(2);
danieleb@1:         if rem(N,2) || ~isreal(x), % odd case
danieleb@1:             % Form intermediate even-symmetric matrix
danieleb@1:             y = zeros(2*N,M);
danieleb@1:             y(1:N,:) = x;
danieleb@1:             y(N+1:2*N,:) = flipud(x);
danieleb@1:             
danieleb@1:             % Compute the FFT and keep the appropriate portion:
danieleb@1:             yy = fft(y);
danieleb@1:             yy = yy(1:N,:);
danieleb@1:         else % even case
danieleb@1:             % Re-order the elements of the columns of x
danieleb@1:             y = [x(1:2:N,:); x(N:-2:2,:)];
danieleb@1:             yy = fft(y);
danieleb@1:             w = 2*w;  % Double the weights for even-length case
danieleb@1:         end
danieleb@1:         % Multiply FFT by weights:
danieleb@1:         c = w(:,ones(1,M)).*yy;
danieleb@13:     case 'III' 
danieleb@13:         % Explicit calculation of DCT III. Replace with fast algorithm (todo)
danieleb@1:         dctMat = zeros(N);
danieleb@1:         for k=0:N-1
danieleb@1:             dctMat(:,k+1) = [sqrt(2)/2; cos(pi/N*(k+1/2)*(1:N-1))'];
danieleb@1:         end
danieleb@1:         dctMat = normcol(dctMat);
danieleb@1:         c = dctMat'*x;
danieleb@1:     case 'IV' % Fast calculation of DCT IV.
danieleb@1:         y = zeros(8*N,M);
danieleb@1:         y(2:2:2*N,:) = x;
danieleb@1:         yy = fft(y,[],1);
danieleb@1:         c = sqrt(2/N)*real(yy(2:2:2*N,:));
danieleb@1:     otherwise
danieleb@1:         error('Unsupported dct type');
danieleb@1: end
danieleb@1: 
danieleb@1: %% Postprocess
danieleb@1: if isreal(x), c = real(c); end
danieleb@1: 
danieleb@1: %EOF