# HG changeset patch
# User danieleb@code.soundsoftware.ac.uk
# Date 1308063167 -3600
# Node ID e124001d2437f76523e4f0271985650f5afcc5f6
# Parent  856ac6a530703d34a82d5f9db74d33adf7527deb
inverse dct

diff -r 856ac6a53070 -r e124001d2437 idcti.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/idcti.m	Tue Jun 14 15:52:47 2011 +0100
@@ -0,0 +1,90 @@
+function x = idcti(c,type)
+%IDCTI inverse discrete cosine transform of types I-IV
+%
+% USAGE
+% x = idcti(c,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
+% c: vector or matrix contaning the inverse dct coefficients
+%
+% REFERENCES
+%   1)
+%
+% SEE ALSO
+% DCTI, DCT, IDCT
+
+%   Author(s): Daniele Barchiesi
+%   Copyright QMUL
+%   $Revision: 1
+
+
+%% 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
+        dctMat = zeros(N);
+        for k=0:N-1
+            dctMat(:,k+1) = [sqrt(2)/2; cos(pi/(N-1)*k*(1:N-2))'; (-1)^k*sqrt(2)/2];
+        end
+        dctMat = normcol(dctMat);
+        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
+        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
\ No newline at end of file