
% plot replica time series
% mbp 8/96, Universidade do Algarve

clear all
ARRFIL = 'ARRFIL';

T = 1.0;		% time window over which impulse response is calculated
c = 1530.0;	% optional reduction velocity

% set up source spectrum

fmin = 380;
fmax = 420;
T = 17.5;
sample_rate = 1200;

[ s, time ] = lfm( fmin, fmax, T, sample_rate );

% hanning weighting to reduce ripple
%s = s .* ( hanning( length( s ) ).^(1/4) )';
wavwrite( s, sample_rate, 'lfm_in' )

shat = fft( s );
shat = shat .* conj( shat );	% propagate the autocorrelation instead

%shat = [ zeros( 1, 300 ) shat zeros( 1, 3199 ) ]';   % in [ 0, 1000 ] Hz
nfreq = size( shat, 2 );

LoadArrfil

% compute the sum

deltat = 1 / sample_rate;
nt = round( T / deltat );	% number of time points


for isd = 1 % :nsd
  H = zeros( nt, nr, nrd );	% clear the transfer function
  for ird = 1 % :nrd
    for ir = 1:nr
      [ isd ird ir ]
      last = narrmat( ir, ird, isd );

      tstart = rr( ir ) / c - 0.1;   % min( delay( ir, 1:last, ird ) )
      tend   = tstart + T - deltat;
      time   = tstart:deltat:tend;

      % compute channel transfer function

      for iarr = 1 : narrmat( ir, ird, isd )

        % compute time index for delay
        Tarr = delay( ir, iarr, ird, isd ) - tstart;
        it = round( Tarr / deltat + 1 );
        % s = Tarr / deltat + 1 - it;		% proportional interval
        if ( it >= 1 & it <= nt )
           H( it,   ir, ird ) = H( it,   ir, ird ) + amp( ir, iarr, ird );
        end

      end   % next arrival, iarr.

    end   % next range, ir.

  end   % next rd

  % form r( t ) as convolution of impulse response with s( t )

  rmod = zeros( nt, nr );

  for ird = 1
    for ir = 1:nr
      gr = real( H( :, ir, ird ) );
      gi = imag( H( :, ir, ird ) );
      grhat = fft( gr );
      giquad = imag( hilbert( gi ) );
      ghat = grhat + fft( giquad ); 

      % note: ifft below has small im. parts due to round-off
      % - makes sign consistent with KRAKEN/SCOOTER
      rmod( :, ir ) = real( ifft( -ghat .* shat.' ) );
    end
  end

  % take log and do some normalization

  rmod = abs( hilbert( rmod ) );

  % normalize

  for ird = 1
     for ir = 1:nr
      scalef = max( rmod( :, ir ) );
        if scalef == 0
         scalef = 1;
      end

      temp = 10 * log10( rmod( :, ir ) / scalef ) + 10;
      I = find( temp < 0 );
      temp( I ) = zeros( size( I ) );
      scalef = norm( temp );
      if scalef == 0
         scalef = 1;
      end
      rmod( :, ir ) = temp / scalef;
    end
  end

  % rmod = rmod( 1:4:4000, :, : );
  nrr = nr;
  fname = [ 'model' int2str( isd ) ]
  eval( [ 'save ' fname ' nrr rr rmod' ] );
end;	% next source depth

tl = 10 * log10( rmod' );
figure
imagesc( time, rr , tl ); zoom on
xlabel( 'Time (s)')
ylabel( 'Range (m)' )
peak = max( max( tl ) );
caxis( [ peak - 30, peak ] );
colorbar
title( 'BELLHOP impulse response' )    

