
% Convolve a source timeseries with the channel impulse response to
% get the received waveform.
%
% This is a delay-and-sum process in which echoes of the source waveform
% are combined based on their amplitude and arrival time
%
% mbp 8/96, Universidade do Algarve
% 4/09 addition of Doppler effects

v  = [ 0, 0 ];   % Receiver motion vector (vr, vz ) in m/s
c0 = 1500;       % reference speed to convert v to proportional Doppler

c = 1520.0;  % reduction velocity (should exceed fastest possible arrival)
T = 10.0;     % time window to capture
sample_rate = 1000;

Narrmx = 100;

fileroot = 'MunkB';

isd = 1;   % source depth index

% source timeseries
% STSFIL = 'Ricker.wav'
% STSFIL = 'sample_band-A_pkt-2_2frames.wav'
% [ stsTemp, sample_rate ] = audioread( STSFIL );
% sts = stsTemp( 1 : 2 : 5 * 96000 );  % sub-sample down to 48000/s
% sample_rate = 48000;

% generated by 'cans'
% sample_rate = 20000;
% TT          = linspace( -0.2, 0.8, 20000 )';
% sample_rate = 20000;
% omega       = 2 * pi * 4000
% Pulse       = 'T'
% [ sts, PulseTitle ] = cans( TT, omega, Pulse );

% read from Sandipa's mat file
% load rate2_48kHz
% sts         = y_pb48k;
% sample_rate = 48000;
% deltat      = 1 / sample_rate;
% nsts         = length( sts );
% TT          = linspace( 0, ( nts - 1 ) * deltat, nts );

% generated by 'cans'
TStart = -1;
TEnd   =  2;
t_sts       = linspace( TStart, TEnd, ( TEnd - TStart ) * sample_rate )';
omega       = 2 * pi * 10;
Pulse       = 'P';
[ sts, PulseTitle ] = cans( t_sts, omega, Pulse );

% normalize the source time series so that it has 0 dB source level
% (measured based on the peak power in the timeseries)
sts = -sts / max( abs( sts ) );   % the - is to match other AT models
nsts = length( sts );   % number of samples in the source time series

% figure
% plot( t_sts, sts )
% xlabel( 'Time (s)' )
% ylabel( 's(t)' )

%%

% optional Doppler effects
% pre-calculate Dopplerized versions of the source time series

% following can be further optimized if we know that ray-angle limits
% further restrict possible Doppler factors

if ( norm( v ) > 0 )
   disp( 'Setting up Dopplerized waveforms' )
   %v_vec = linspace( min( v ), max( v ), 10 );   % vector of Doppler velocity bins
   v_vec       = linspace( 1.9, 2, 51 );   % vector of Doppler velocity bins
   alpha_vec   = 1 - v_vec / c0;                     % Doppler factors
   nalpha      = length( alpha_vec );
   sts_hilbert = zeros( nsts, nalpha );
   
   % loop over Doppler factors (should be further vectorized ...)
   for ialpha = 1 : length( alpha_vec )
      disp( ialpha )
      sts_hilbert( :, ialpha ) = hilbert( arbitrary_alpha( sts', alpha_vec( ialpha ) )' ); % Dopplerize
   end
else
   sts_hilbert = hilbert( sts );   % need the Hilbert transform for phase changes (imag. part of amp)
   % sts_hilbert = sts;
end

figure
plot( t_sts, imag( sts_hilbert ) )
xlabel( 'Time (s)' )
ylabel( 's^+(t)' )

%%
% Main routine:
% read the arrivals file and convolve with the source waveform

ARRFIL   = [ fileroot '.arr' ];

[ Arr, Pos ] = read_arrivals_asc( ARRFIL );  % read the arrivals file
disp( 'Done reading arrivals' )
omega = 2 * pi * Pos.freq;

% select which source/receiver index to use
nsd = length( Pos.s.z );
nrr = length( Pos.r.r );
nrd = length( Pos.r.z );

deltat = 1 / sample_rate;
nrts   = round( T / deltat );	  % number of samples in the receiver time series
rtsmat = zeros( nrts, nrd, nrr );

for ird = 1 : nrd   % loop over receiver depths
   disp( ird )
   
   for irr = 1 : nrr   % loop over receiver ranges
      tstart = Pos.r.r( irr ) / c - 0.1;   % min( delay( ir, :, ird ) )

      % compute channel transfer function
      
      rts = zeros( nrts, 1 );	% initialize the time series
      
      Narr = Arr( irr, ird, isd ).Narr;
        
      for iarr = 1 : Narr   % loop over arrivals

         % delay is a complex value with the imaginary part holding loss due to volume attenuation
         % Thus, e( -i omega * delay ) has a phase delay and an amplitude decay.
         % Need to check that omega is a central frequency and test the VolumeLoss

         Tarr = real( Arr( irr, ird, isd ).delay( iarr ) ) - tstart;   % arrival time relative to tstart
         VolumeLoss = exp( omega * imag( Arr( irr, ird, isd ).delay( iarr ) ) );
         
         irts1 = round( Tarr / deltat + 1 );            % starting time index in rts for that delay
         irts2 = irts1 + nsts - 1;                      % ending   time index in rts
         ists1 = 1;                                      % starting time index in sts
         ists2 = nsts;                                   % ending   time index in sts
         
         % clip to make sure [ it1, it2 ] is inside the limits of rts
         if ( irts1 < 1 )
            ists1 = ists1 + ( 1 - irts1 );  % shift right by 1 - it1 samples
            irts1  = 1;
         end
         
         if ( irts2 > nrts )
            ists2 = ists2 - ( irts2 - nrts );  % shift left by it2 - nt samples
            irts2  = nrts;
         end
         
         if ( norm( v ) > 0 )
            % identify the Doppler bin
            theta_ray = Arr( irr, ird, isd ).RcvrDeclAngle( iarr ) * pi / 180;   % convert ray angle to radians
            tan_ray   = [ cos( theta_ray ) sin( theta_ray ) ];        % form unit tangent
            alpha     = 1 - dot( v / c0, tan_ray );                   % project Doppler vector onto ray tangent
            ialpha    = 1 + round( ( alpha - alpha_vec( 1 ) ) / ( alpha_vec( end ) - alpha_vec( 1 ) ) * ( length( alpha_vec ) - 1 ) );
            
            % check alpha index within bounds?
            if ( ialpha < 1 || ialpha > nalpha )
               disp( 'Doppler exceeds pre-tabulated values' )
               ialpha = max( ialpha, 1 );
               ialpha = min( ialpha, nalpha );
            end
            
            % load the weighted and Dopplerized waveform into the received time series
            rts( irts1 : irts2 ) = rts( irts1 : irts2 ) + VolumeLoss * ...
               real( Arr( irr, ird, isd ).A( iarr ) * sts_hilbert( ists1 : ists2, ialpha ) );
         else
            
            % rts( it1 : it2 ) = rts( it1 : it2 ) + real( Arr( irr, ird, isd ).A( iarr ) ) * real( sts_hilbert( its1 : its2, 1 ) ) ...
            %                                     - imag( Arr( irr, ird, isd ).A( iarr ) ) * imag( sts_hilbert( its1 : its2, 1 ) );
            
            % following is math-equivalent to above, but runs faster in
            % Matlab even though it looks like more work ...
            rts( irts1 : irts2 ) = rts( irts1 : irts2 ) + VolumeLoss * ...
               real( Arr( irr, ird, isd ).A( iarr ) * sts_hilbert( ists1 : ists2, 1 ) );
         end
      end   % next arrival, iarr
      
      %       % write to file
      %       if ( ird == 1 && ir == 1 )
      %          fid_out = fopen( pcmfile, 'w', 'ieee-le' );
      %       else
      %          fid_out = fopen( pcmfile, 'a', 'ieee-le' );
      %       end
      %       if ( fid_out == -1 )
      %          error('Can''t open PCM file for output.');
      %       end
      %
      %       fwrite( fid_out, 2^15 * 0.95 * rts / max( abs( rts ) ), 'int16' );
      %       fclose( fid_out );
      
      %wavwrite( 0.95 * rts / max( abs( rts ) ) , sample_rate, [ 'rts_Rd_' num2str( ird ) '_Rr_' num2str( ir ) '.wav' ] );
      %rtsmat( :, ir ) = 0.95 * rts / max( abs( rts ) );
      rtsmat( :, ird, irr ) = rts;
   end
   %eval( [ ' save ' fileroot '_Rd_' num2str( ird ) ' rtsmat Pos sample_rate' ] );
end

%%
% plot snapshots

% figure
% for it = 1 : nt
%    imagesc( squeeze( rtsmat( it, :, : ) ) )
%    colorbar
%    drawnow
% end

%%
% save to a Matlab timeseries file

irr = 26;   % select range
%irr = 6;   % select range

% calculate time vector that goes with that range
tstart = Pos.r.r( irr ) / c - 0.1;   % min( delay( ir, :, ird ) )
tend   = tstart + T - deltat;
tout   = tstart : deltat : tend;

PlotTitle = fileroot;

RTS = -squeeze( rtsmat( :, :, irr ) );
save( [ fileroot '.rts.mat' ], 'PlotTitle', 'Pos', 'tout', 'RTS' )

% save autec
%
% figure; pcolor( Pos.r.r, linspace( 0, T, nt ), 20 * log10 ( abs( hilbert( rtsmat ) )' ) );
% shading flat
% colorbar
% xlabel( 'Time (s)' )
% ylabel( 'Range (m)' )
