% AC_LP_WTT - create range-lag ambiguity plot for alternating code lag profiles
% function hs = ac_lp_wtt(codeset [, baud [, 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
%  - baud is the length of a transmitted baud in basic time units
%	(default: 20)
%  - lags is the lags that will be computed, in units of 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, AC_GATE_WTT, AC_LPH_WTT, ALTCODE, BCAC_LP_WTT, MP_LP_WTT.

%   $Id: ac_lp_wtt.m,v 1.4 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 = ac_lp_wtt(codeset, baud, lags, ax, plotfunc, varargin)

if nargin < 1, help ac_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(baud), baud = 20;        end
if nargin < 3 | isempty(lags), lags = 0:nbaud-1; end
if nargin < 4                  ax   = gca;       end
if nargin < 5 | isempty(plotfunc),
    plotfunc = 'imagesc';
else
    if isempty(ax), ax = gca; end
end

h = baudexpand(1, baud)/baud;

for ii = 1:scancount,
    env      = fliplr(baudexpand(codeset(:,ii), baud));
    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,
    ptt = zeros(2*baud-1, 2*baud*(nbaud-lgind));
    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
    % fprintf('done lag %d\n', lgind)
    % pause
    
    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
	    % lwtt{jj} = zeros(size(cwtt{1}));
	    %for ii = 1:scancount,
	    %    lwtt{jj} = lwtt{jj} + decode(ii)*cwtt{ii};
	    %end
	    % subplot(2, nplot, nplot+jj); imagesc(wtmp'); axis xy
	    %ptt(:, baud*(jj-1)+[1:baud*(nbaud-lgind+1)]) = ...
	    %    ptt(:, baud*(jj-1)+[1:baud*(nbaud-lgind+1)]) ...
	    %    + lwtt{jj};

	    % wtt will contain the maximum non-zero value
	    % wtmp = wtmp/((nbaud-lgind)*scancount);
	    % rind = baud*(lgind-1)+[1:size(ptt,2)];
	    tmp  = wtt(lag+[1:2*baud-1], lag+[1:size(wtmp,2)]);
	    ii   = find(tmp == 0 | (wtmp ~= 0 & wtmp > tmp));
	    tmp(ii) = wtmp(ii);
	    wtt(lag+[1:2*baud-1], lag+[1:size(wtmp,2)]) = tmp;
	end
	% pause
    end
end

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

wtt = wtt / max(wtt(:));

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

if nargout > 0, hs = h; end