% BCMP_LP_WTT - create range-lag ambiguity figure for Barker-coded multipulse code lag profile
% function hs = bcmp_lp_wtt(pattern [, len [, baud [, lags [, ax [, plotfunc]]]]])
%  - pattern is the intervals (gaps) sequence for the multipulse code
%	[Farley, /Radio Science/ 7(6), 661, 1972], table 1
%  - len is the Barker code or the length of the Barker code
%  - baud is the length of a transmitted baud in basic time units
%	(default: 1)
%  - 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, MULTIPULSE, MP_LP_WTT, BARKERCODE, BC_WTT, BCAC_LP_WTT.

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

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


function hs = bcmp_lp_wtt(pattern, len, baud, lags, ax, plotfunc, varargin)

if nargin < 1, help bcmp_lp_wtt, return, end

if nargin < 2 | isempty(len),   len = 13;             end
if nargin < 3 | isempty(baud), baud = 3;              end
if nargin < 4 | isempty(lags), lags = 0:sum(pattern); end
if nargin < 5                    ax = gca;            end
if nargin < 6 | isempty(plotfunc),
    plotfunc = 'imagesc';
else
    if isempty(ax), ax = gca; end
end

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

cbaud = (length(code))*baud;

h   = baudexpand(code, baud);
env = fliplr(baudexpand(multipulse(pattern, code), baud));

% wta(tau, r) = h(t-tau) * env(tau-r)
wta = zeros(length(h), length(env)+length(h));

for tau = 1:length(h)
    wta(tau,tau:tau+length(env)-1) = (h(tau) * env);
end

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

for lgind = lags
    if lgind == 0,
	tmp = convolve(wta', wta')';
	wtt(1:2*cbaud-1,:) = tmp;
    else 
	lag = (cbaud)*lgind;

	% wtt contains the maximum non-zero entry
	tmp  = convolve(wta(:,lag+1:end)', wta(:,1:end-lag)')';
	wtmp = wtt(lag+[1:2*cbaud-1],lag+1:end);
	ii   = find(wtmp == 0 | (tmp ~= 0 & tmp > wtmp));
	wtmp(ii) = tmp(ii);

	wtt(lag+[1:2*cbaud-1],lag+1:end) = wtmp;
    end
end


% ii = find(wtt == 0); wtt(ii) = NaN*ones(size(ii));

% 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), wtt', varargin{:});
axis xy;

if nargout > 0, hs = h; end