% AC_GATE_WTT - create range-lag ambiguity plot for an alternating code gate
% function hs = ac_gate_wtt(codeset [, baud [, lags [, weight [, 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 non-zero integer lags)
%  - weight is a scalar or a vector with as many elements as there are
%       rows in 'codeset', giving the weight to use for successive pulses
%       (e.g. to simulate time-varying scatterers)  (default: 1)
%  - 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_LP_WTT, ALTCODE, BCAC_GATE_WTT.

%   $Id: ac_wgate_wtt.m,v 1.1 2008-11-18 04:57:57 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_gate_wtt(codeset, baud, lags, weight, ax, plotfunc, varargin)

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

if nargin < 5 | isempty(plotfunc),
    plotfunc = 'imagesc';
else
    if isempty(ax), ax = gca; end
end
if nargin < 6 | isempty(weight),
    weight = ones(1, scancount);
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
end

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

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

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

% ii = find(wtt <= 0); wtt(ii) = 1e-20*ones(size(ii));

% imagesc(10*log10(wtt')); set(gca, 'clim', [-20 0]);

wtt = wtt / max(wtt(:));
% 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*(nbaud-1))/(baud), wtt', varargin{:});
axis xy;

if nargout > 0, hs = h; end