function [C,C0] = corrmat(x,lags,method,symmetrize);
%
%CORRMAT compute lagged  correlation matrices
%
%   Use corrmat to compute a set of correlation matrices to be diagonalized
%   with qdiag or any other matrix diagonalization algorithm.
%
%   C = corrmat(x, lags), where x is a TxN matrix and lags is a K-element
%   vector of whole numbers, computes K lagged NxN correlation matrices of
%   the row vectors in x.
%
%   C = corrmat(x, lags, method) uses method to compute the
%   matrices. method can be 'cov' (compute covariance matrices, default)
%   or 'moments' (compute matrices of second moments).
%  
%   C = corrmat(x, lags, method, symmetrize), where symmetrize is a
%   scalar boolean, can be used to explicitely symmetrize the resulting
%   correlation matrices (default: symmetrize==false).
%
%   [C,C0] = corrmat(x, lags, ... )  additionally returns in C0 the matrix
%   C(:,:,k) for which lags(k) = 0 (if there is any).
%  
%   See also qdiag.


% Copyright (C) 2006  TU-Berlin, Roland Vollgraf
% 
% This program is free software; you can redistribute it and/or
% modify it under the terms of the GNU General Public License
% as published by the Free Software Foundation; either version 2
% of the License, or (at your option) any later version.
% 
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
% 
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
% 
% If you use this software, please cite 
% "R. Vollgraf and K. Obermayer, Quadratic Optimization for Simultaneous
% Matrix Diagonalization, IEEE Transaction on Signal Processing, 2006, in
% press"
% in your publications
  
  
  K = length(lags);
  
  
  [T,N] = size(x);
  if T<N
    warning('x may be oriented the wrong way.');
  end
  
  
  C = zeros(N,N,K);
  
  
  if nargin<3
    method = 'cov';
  end

  if nargin<4
    symmetrize = false;
  end
  

  for k=1:K
    
    if lags(k)>=0
      x1 = x(1:end-lags(k),:);
      x2 = x(1+lags(k):end,:);
    else
      x1 = x(1-lags(k):end,:);
      x2 = x(1:end+lags(k),:);
    end  

    T = size(x1,1);

    switch method
      
      case 'cov'
        x1 = x1-repmat(sum(x1,1)/T,T,1);
        x2 = x2-repmat(sum(x2,1)/T,T,1);
        C(:,:,k) = x1'*x2/(T-1);
      
      case 'moments'
        C(:,:,k) = x1'*x2/T;
        
      otherwise
        error(sprintf('%s: unknown method.',char(method)));
        
    end
    
    if symmetrize
      C(:,:,k) = (C(:,:,k)+C(:,:,k)')/2;
  
    if lags(k)==0
      C0 = C(:,:,k);
    end
  end
  
  
  
  
  
  
  
  
