% BCAC_LP_WTT - create range-lag ambiguity plot for Barker-coded alternating code lag profiles
% function hs = bcac_lp_wtt(codeset [, bclen [, bcbaud [, lags [, ax [, plotfunc]]]]])
%  - codeset is the sign sequence for the alternating code,
%	one scan per row, or
%     + a number of bits (4, 8, 16) for a strong codeset from
%       [Nygrén, 1996, table 6.2], or 
%     + the string 'tau0' for the codeset used in that experiment
%  - bclen is the length of the Barker code in bauds (default: 13)
%  - bcbaud is the length of a Barker code baud in basic time units
%	(default: 3)
%  - lags is the lags that will be computed, in units of bclen*baud
%	(default: all integer lags)
%  - ax is the axis in which to make the plot (default: current axis)
%       if empty and plotfunc is not given, no plot is made and wtt is returned
%  - plotfunc is the function to use for the ambiguity figure
%       (default: imagesc)
%
%   See also WTT, BCAC_GATE_WTT, ALTCODE, AC_LP_WTT, BCMP_LP_WTT,
%            BARKERCODE, BC_WTT.

%   $Id: bcac_lp_wtt.m,v 1.3 2008-11-18 04:57:14 tom Exp $

%   Copyright © 2003-2008 Tom Grydeland <Tom.Grydeland@phys.uit.no>
%   This is free software, licensed under GNU GPL version 2 or later

function hs = bcac_lp_wtt(codeset, bclen, bcbaud, lags, ax, plotfunc, varargin)

if nargin < 1, help bcac_lp_wtt, return, end

if length(codeset) == 1,
    codeset = altcode(codeset);
elseif strncmp(codeset, 'tau0', 4)
    codeset = tau0code;
end

codeset = codeset';

[nbaud, scancount] = size(codeset);

if nargin < 2 | isempty(bclen),  bclen  = 13;        end
if nargin < 3 | isempty(bcbaud), bcbaud =  3;        end
if nargin < 4 | isempty(lags),     lags = 0:nbaud-1; end
if nargin < 5                        ax = gca;       end
if nargin < 6 | isempty(plotfunc),
    plotfunc = 'imagesc';
else
    if isempty(ax), ax = gca; end
end

if length(bclen) > 1
    code = bclen;
else
    code = barkercode(bclen);
end

baud = (length(code))*bcbaud;
h    = baudexpand(code, bcbaud);

for ii = 1:scancount,
    env      = fliplr(baudexpand(codeset(:,ii), baudexpand(code, bcbaud)));
    cwta{ii} = zeros(length(h), length(env)+length(h));
    
    for tau = 1:length(h)
	cwta{ii}(tau,tau:tau+length(env)-1) = (h(tau) * env);
    end
    % imagesc(cwta{ii}'); axis xy; pause
end

wtt = zeros(baud*(max(lags)+2), length(env)+length(h));

nplot = max(size(codeset));
for lgind = lags,
    for ii = 1:scancount,
	if lgind == 0,
	    cwtt{ii} = convolve(cwta{ii}', cwta{ii}')';
	    wtt(1:2*baud-1, :) = wtt(1:2*baud-1, :) + cwtt{ii};
	else
	    lag = baud*lgind;
	    cwtt{ii} = ...
		convolve(cwta{ii}(:,lag+1:end)', cwta{ii}(:,1:end-lag)')';
	end
	% subplot(2, nplot, ii); imagesc(cwtt{ii}'); axis xy
    end
    
    if lgind ~= 0,
	% decode alt. codes
	for jj = 1:nbaud-lgind,
	    decode = codeset(jj, :) .* codeset(jj+lgind, :);
	    wtmp = zeros(size(cwtt{1}));
	    for ii = 1:scancount,
		wtmp = wtmp + decode(ii)*cwtt{ii};
	    end
	    % subplot(2, nplot, nplot+jj); imagesc(wtmp'); axis xy

	    % wtt contains the maximum non-zero value
	    tmp = wtt(lag+[1:2*baud-1], lag+1:end);
	    ii  = find(tmp == 0 | (wtmp ~= 0 & wtmp > tmp));
	    tmp(ii) = wtmp(ii);
	    wtt(lag+[1:2*baud-1], lag+1:end) = tmp;
	end
	% pause
    end
end

% If ax is empty, return wtt'
if isempty(ax), hs = wtt'; return; end

% Otherwise, plot wtt' and return handle to figure
axes(ax);
h = feval(plotfunc, ([1:size(wtt,1)]-baud)/(baud), ...
	[1:size(wtt,2)]/(baud), wtt', varargin{:});
axis xy;

if nargout > 0, hs = h; end