clc;
clf;
clear all;
warning('OFF', 'MATLAB:table:ModifiedAndSavedVarnames')

set(figure(1), 'Renderer', 'painters');

myfile = uigetfile('*.csv','MultiSelect','on');
M1 = readtable(myfile);
A1 = table2array(M1);
A1(1,:) = [];

Time = A1(:,1);
Input = A1(:,2);
Output = A1(:,3);

Timeplot = Time*1e6;

filterbw = [0.6e6 3e6];
filterlength = 2000;

%d= distance between wedges + distance to beam (centre of wedge)
df = finddelay(Input, Output);
Fs = 1/((Time(2,1))-(Time(1,1)));
St = 1/Fs;
dftimeToF = df*St;

%low pass filter
bpinput = bandpass(Input, filterbw, Fs,'ImpulseResponse','iir');

%low pass filter
bpoutput = bandpass(Output, filterbw, Fs,'ImpulseResponse','iir');

%normalise data to -1 - 1
factorin = max([max(bpinput) abs(min(bpinput))]);
normin = bpinput/factorin;
factorout = max([max(bpoutput) abs(min(bpoutput))]);
normout = bpoutput/factorout;

%create envelope of Input
[upin,loin] = envelope(normin,filterlength,'analytic');

%find peak of envelope Input
[pkkin,locsin] = findpeaks(upin,Fs,'MinPeakHeight',0.5);
peaklocin = locsin(1,1)+min(Time);

%create envelope of Output
[upout,loout] = envelope(normout,filterlength,'analytic');

%find peak of envelope Output
[pkkout,locsout] = findpeaks(upout,Fs,'MinPeakHeight',0.5);
peaklocout = locsout+min(Time);

%calc ToF
ToF = (peaklocout(1,1)-peaklocin)

%fft
Fs = 1/((Time(2,1))-(Time(1,1)));
St = 1/Fs;     % Sampling period
                   
L = length(bpoutput);   % Length of signal
t = (0:L-1)*St;        % Time vector

Y1 = fft(bpinput);
Y2 = fft(bpoutput);

PI2 = abs(Y1/L);
PI1 = PI2(1:L/2+1);
PI1(2:end-1) = 2*PI1(2:end-1);

PO2 = abs(Y2/L);
PO1 = PO2(1:L/2+1);
PO1(2:end-1) = 2*PO1(2:end-1);

f = Fs*(0:(L/2))/L;

figure(1)
hold on
grid on
%title('Simulated Signals')
xlim([0 max(Timeplot)])

yyaxis left
plot(Timeplot,normin)
plot(Timeplot,upin,'-','LineWidth',4.0)
ylim([-1.2 1.2])
yticks(-1.2:0.2:1.2)
ylabel('Normalised Amplitude')

plot([peaklocin*1e6 peaklocin*1e6],[0 max(normin)],'b:','LineWidth',4.0)
plot([peaklocout(1,1)*1e6 peaklocout(1,1)*1e6],[0 max(normout)],'b:','LineWidth',4.0)

yyaxis right
hold on
plot(Timeplot,normout)
plot(Timeplot,upout,'-','LineWidth',4.0)
ylabel('Normalised Amplitude')
set(gca,'fontname','Libertinus Sans','fontsize', 20)
ylim([-1.2 1.2])
yticks(-1.2:0.2:1.2)
xlabel('Time (μs)')
hold off

x0=10;
y0=10;
width=1600;
height=1000;
set(gcf,'position',[x0,y0,width,height])
set(gca,'fontname','Libertinus Sans','fontsize', 20)

% %plot FFT
% figure(2)
% 
% tiledlayout(2,1)
% ax1=nexttile;
% 
% plot(f,PI1) 
% title('Single-Sided Amplitude Spectrum of Excitation Pulse')
% xlabel('f (Hz)')
% xlim([0 2e6])
% ylabel('|P1(f)|')
% grid on
%  
% ax2=nexttile;
% plot(f,PO1) 
% title('Single-Sided Amplitude Spectrum of Output')
% xlabel('f (Hz)')
% xlim([0 2e6])
% ylabel('|P1(f)|')
% grid on
% 
% x0=10;
% y0=10;
% width=1600;
% height=1000;
% set(gcf,'position',[x0,y0,width,height])

set(gcf, 'Color', 'w');
%export_fig simpulse20csensors.eps -nocrop -transparent