clc
clear

load W3_SaltWater.csv
RawData = W3_SaltWater;
clear W3_SaltWater

%load W2_PureWater.csv
%RawData = W2_PureWater;
%clear W2_PureWater

Rx = RawData(:,1);
Tx = RawData(:,2);
clear RawData

SPR = 64.0E+6;
TotalPnts = length(Tx);
dt = 1.0E+6/SPR;
tt = [1:TotalPnts]*dt;

TxMean = mean(Tx)
RxMean = mean(Rx)

Tx = Tx - TxMean;
Rx = RxMean - Rx;

scaling = 0.2;
TxPower = mean(abs(Tx))*scaling;
RxPower = mean(abs(Rx))*scaling;



bPeakFound = 0;  PeakCount = 0;
PeakVal = 0;  PeakIndex = 0;  IndexCount = 0;
for kk = 1:TotalPnts
  if ( bPeakFound == 1 )  bPeakFound = 0; PeakVal = 0;  PeakIndex = 0;  IndexCount = 0;  end
  Value =  Tx(kk); 
  if ( Value > TxPower ) 
      if ( PeakVal < Value ) 
          PeakVal = Value;  PeakIndex = kk;  IndexCount = IndexCount + 1;
      end    
  else
      if (IndexCount > 1) bPeakFound = 1; PeakCount = PeakCount+1; Peaks(PeakCount) = PeakIndex; end
  end
end  

TxTrigger = Peaks(4)
TxPeak = zeros(TotalPnts,1);
for jj = 1:PeakCount
    TxPeak(Peaks(jj)) = 200;
end  

bPeakFound = 0;  PeakCount = 0;
PeakVal = 0;  PeakIndex = 0;  IndexCount = 0;
for kk = 1:TotalPnts
  if ( bPeakFound == 1 )  bPeakFound = 0; PeakVal = 0;  PeakIndex = 0;  IndexCount = 0;  end
  Value =  Rx(kk); 
  if ( Value > TxPower ) 
      if ( PeakVal < Value ) 
          PeakVal = Value;  PeakIndex = kk;  IndexCount = IndexCount + 1;
      end    
  else
      if (IndexCount > 1) bPeakFound = 1;  PeakCount = PeakCount+1;  Peaks(PeakCount) = PeakIndex; end
  end
end  

RxTrigger = Peaks(4)
RxPeak = zeros(TotalPnts,1);
for jj = 1:PeakCount
    RxPeak(Peaks(jj)) = 200;
end  


ax = [8 15 -256 256 ];
figure(1)
subplot(211), plot(tt,Tx,'-c',tt,TxPeak,'-b'), axis(ax), title('Transmit Waveform')
subplot(212), plot(tt,Rx,'-c',tt,RxPeak,'-r'), axis(ax), title('Received Waveform')

%TxTrigger*dt
%RxTrigger*dt
BigShift = RxTrigger - TxTrigger;
BigShift = BigShift - mod(BigShift,16)   % cycle


%start = 964;
start = RxTrigger;
span = 28*16;
stop = start + span - 1;
CX1 = Tx(start:stop).*Rx(start:stop);
CX2 = Tx(start:stop).*Rx(start+4:stop+4);
M1 = sum(CX1)
M2 = sum(CX2)
AtanV = M2/M1
Theta = atan(AtanV);                    % in radian
if ( (M2 > 0) & (M1 > 0) ) Phase = Theta;           end
if ( (M2 > 0) & (M1 < 0) ) Phase = Theta + pi;      end
if ( (M2 < 0) & (M1 < 0) ) Phase = Theta + pi;      end
if ( (M2 < 0) & (M1 > 0) ) Phase = Theta + 2*pi;    end

PhaseShift = 8*Phase/pi                 % in sampling points
PhaseShiftDeg = Phase/pi*180            % in degree
TotalShift = BigShift + PhaseShift      % in sampling points

BigTime = BigShift*dt
SmallTime = PhaseShift*dt
TimeOfFlight = BigTime + SmallTime


%figure(3)
%subplot(211), plot(CX1,'-r'), title('CX1')
%subplot(212), plot(CX2,'-b'), title('CX2')

