spm_dctmtx

 Creates basis functions for Discrete Cosine Transform.
 FORMAT C = spm_dctmtx(N,K,n)
     OR C = spm_dctmtx(N,K)
     OR D = spm_dctmtx(N,K,n,'diff')
     OR D = spm_dctmtx(N,K,'diff')
 N - dimension
 K - order
 n - optional points to sample
____________________________________________________________________________
 spm_dctmtx creates a matrix for the first few basis functions of a one
 dimensional discrete cosine transform.
 With the 'diff' argument, spm_dctmtx produces the derivatives of the
 DCT.

 See:    Fundamentals of Digital Image Processing (p 150-154).
         Anil K. Jain 1989.
____________________________________________________________________________
 @(#)spm_dctmtx.m    2.1 John Ashburner MRCCU/FIL 01/08/28

0001 function C = spm_dctmtx(N,K,n,f)
0002 % Creates basis functions for Discrete Cosine Transform.
0003 % FORMAT C = spm_dctmtx(N,K,n)
0004 %     OR C = spm_dctmtx(N,K)
0005 %     OR D = spm_dctmtx(N,K,n,'diff')
0006 %     OR D = spm_dctmtx(N,K,'diff')
0007 % N - dimension
0008 % K - order
0009 % n - optional points to sample
0010 %____________________________________________________________________________
0011 % spm_dctmtx creates a matrix for the first few basis functions of a one
0012 % dimensional discrete cosine transform.
0013 % With the 'diff' argument, spm_dctmtx produces the derivatives of the
0014 % DCT.
0015 %
0016 % See:    Fundamentals of Digital Image Processing (p 150-154).
0017 %         Anil K. Jain 1989.
0018 %____________________________________________________________________________
0019 % @(#)spm_dctmtx.m    2.1 John Ashburner MRCCU/FIL 01/08/28
0020 
0021 d = 0;
0022 
0023 if nargin == 1, K = N; end;
0024 
0025 if any(nargin == [1 2]),
0026     n = (0:(N-1))';
0027 elseif nargin == 3,
0028     if strcmp(n,'diff'),
0029         d = 1;
0030         n = (0:(N-1))';
0031     elseif strcmp(n,'diff2'),
0032         d = 2;
0033         n = (0:(N-1))';
0034     else
0035         n = n(:);
0036     end
0037 elseif nargin == 4,
0038     n = n(:);
0039     if strcmp(f,'diff'),
0040         d = 1;
0041     elseif strcmp(n,'diff2'),
0042         d = 2;
0043     else
0044         error('Incorrect Usage');
0045     end
0046 else
0047     error('Incorrect Usage');
0048 end
0049 
0050 C = zeros(size(n,1),K);
0051 
0052 if d == 0,
0053     C(:,1)=ones(size(n,1),1)/sqrt(N);
0054     for k=2:K
0055         C(:,k) = sqrt(2/N)*cos(pi*(2*n+1)*(k-1)/(2*N));
0056     end
0057 elseif d == 1,
0058     for k=2:K
0059         C(:,k) = -2^(1/2)*(1/N)^(1/2)*sin(1/2*pi*(2*n*k-2*n+k-1)/N)*pi*(k-1)/N;
0060     end
0061 elseif d == 2,
0062     for k=2:K,
0063         C(:,k) = -2^(1/2)*(1/N)^(1/2)*cos(1/2*pi*(2*n+1)*(k-1)/N)*pi^2*(k-1)^2/N^2;
0064     end;
0065 else,
0066     error('Can''t do this');
0067 end
0068 
0069