changeset 1:856ac6a53070

Added dct function

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/dcti.m	Tue Jun 14 15:27:00 2011 +0100
@@ -0,0 +1,80 @@
+function c = dcti(x,type)
+%DCTI discrete cosine transform of types I-IV
+%
+% USAGE
+% c = dcti(x,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)
+%
+% SEE ALSO
+% DCT, FFT, IDCT, IFFT
+
+%   Author(s): Daniele Barchiesi
+%   Copyright QMUL
+%   $Revision: 1
+
+
+%% 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
+        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);
+        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
+        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
\ No newline at end of file