function r = chromxcorr(Q,T,L)
% r = chromxcorr(Q,T,L)
%   Cross-correlate two chroma ftr vecs in both time and
%   transposition
%   Both Q and T can be long, result is full convolution
%   (length(Q) + length(T) - 1 columns, in T order).
%   L is the maximum lag to search to - default 100.
%   Optimized version.
% 2006-07-14 dpwe@ee.columbia.edu

%   Copyright (c) 2006 Columbia University.
% 
%   This file is part of LabROSA-coversongID
% 
%   LabROSA-coversongID is free software; you can redistribute it and/or modify
%   it under the terms of the GNU General Public License version 2 as
%   published by the Free Software Foundation.
% 
%   LabROSA-coversongID 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 LabROSA-coversongID; if not, write to the Free Software
%   Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA
%   02110-1301 USA
% 
%   See the file "COPYING" for the text of the license.

if nargin < 3;  L = 100; end

[nchr,nbts1] = size(Q);
[nchr2,nbts2] = size(T);

if nchr ~= nchr2
  error('chroma sizes dont match');
end

% % If test item is longer than query, truncate it to 
% % avoid excessively large xcorrs from very long test items
% if nbts2 > nbts1
%   % take the middle part of T
%   T = T(:,floor((nbts2-nbts1)/2)+[1:nbts1]);
%   nbts2 = nbts1;
% end

r = zeros(nchr, 2*L+1);

%for i = 1:nchr
%  rr = 0;
%  for j = 1:nchr
%    rr = rr + xcorr(T(1+rem(j+i-2,nchr),:),Q(j,:),L);
%  end
%  r(i,:) = rr;
%end

% super speed up by Jesper hoevang jensen jhj@es.aau.dk 2007-05-04
% execution time of 80x80 comparison on macbook pro goes from 3700s 
% to 336s!!!  (only 891s to 418s on hog)

t=max(length(Q), length(T))+2*L+1;
t2=ifft2(fft2(T, nchr, t).*conj(fft2(Q, nchr, t)));
r = [t2(:,end-L+1:end) t2(:,1:L+1)];


% Normalize by shorter vector so max poss val is 1
%r = r/min(nbts1,nbts2);

% with normalization: cover80 = 46.25% (42.5)
% without normalization: cover80 = 57.5% (51.1)