% BCAC_GATE_WTT - create range-lag ambiguity plot for a Barker-coded alternating code gate
% function hs = bcac_gate_wtt(codeset [, bclen [, 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
%  - bclen is the length of the Barker code in bauds (default: 13)
%  - baud is the length of a transmitted baud in basic time units
%	(default: 3)
%  - lags is the lags that will be computed, in units of baud
%	(default: all non-zero 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_LP_WTT, ALTCODE, AC_GATE_WTT, BARKERCODE, BC_WTT.

%   $Id: bcac_gate_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_gate_wtt(codeset, bclen, baud, lags, ax, plotfunc, varargin)

if nargin < 1, help bcac_gate_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(baud),   baud = 3;         end
if nargin < 4 | isempty(lags),   lags = 1: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

cbaud = (length(code))*baud;

h   = baudexpand(code, baud);

for ii = 1:scancount,
    env      = fliplr(baudexpand(codeset(:,ii), baudexpand(code, 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
end

% The Wtt' range-lag ambiguity functions for a complete range gate
wtt = zeros(cbaud*(max(lags)+2), 2*cbaud*(nbaud-1));

for lgind = lags,
    ptt = zeros(2*cbaud-1, 2*cbaud*(nbaud-lgind));
    for ii = 1:scancount,
	if lgind == 0,
	    cwtt{ii} = convolve(cwta{ii}', cwta{ii}')';
	else
	    lag = cbaud*lgind;
	    cwtt{ii} = ...
		convolve(cwta{ii}(:,lag+1:end)', cwta{ii}(:,1:end-lag)')';
	end
	% subplot(2, scancount, 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, scancount, scancount+jj); imagesc(wtmp'); axis xy

	    % Range indices slightly tricky
	    rind = cbaud*(jj-1)+[1:cbaud*(nbaud-lgind+1)];
	    ptt(:, rind) = ptt(:, rind) + wtmp;
	    % imagesc(wtmp'); axis xy; pause
	end
	wtmp = ptt/(nbaud-lgind);
	% imagesc(ptt'); axis xy; pause

	% wtt contains the maximum non-zero value, range index tricky again
	rind = cbaud*(lgind-1)+[1:size(wtmp,2)];
	tmp = wtt(lag+[1:2*cbaud-1], rind);
	ii  = find(tmp == 0 | (wtmp ~= 0 & wtmp > tmp));
	tmp(ii) = wtmp(ii);
	wtt(lag+[1:2*cbaud-1], rind) = tmp;
    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)]-cbaud)/(cbaud), ...
	([1:size(wtt,2)]-cbaud*(nbaud-1))/(cbaud), wtt', varargin{:});
axis xy;

if nargout > 0, hs = h; end