%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Water probe measurement reconstruction
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear;
clf;

fcal=input('Calibration file name [.csv]? ','s');
fin=input('Measure file name [.csv]? ','s');

[G,f,srate,np]=load_fresp(strcat(fcal,'.csv'));
disp(sprintf('Calibration file resolution = %d ns',srate/1e-9));

[m,x,rrate]=wave_restore(strcat(fin,'.csv'),G,srate);

%lowpass filter at 12 MHz
b=fir1(256,[12e6*2*rrate],'low');
j=filtfilt(b,1,x);

[lstep,x]=size(j);

figure(5);
t=0:rrate:rrate*(lstep-1);
plot(t,m(1:lstep),t,j);
xlabel('Time');
ylabel('Amplitude');
legend('measured','reconstructed',0);
title(sprintf('FILTERING %s WITH %s',fin,fcal));


% Loads a frequency response from a csv file
%
% fin = input file name
% fvals = frequency values (sized L/2+1)
% G = fresp matrix (sized L)
% srate = sampling rate
% np = number of useful points (<>1)

function [G,fvals,srate,np]=load_fresp(fin)

e=csvread(fin);
[sz,x]=size(e);
np=(sz-1)/2; %n. of useful points

srate=e(1,1);
len=e(1,2); %original length
G=zeros(len,1);

fvals=1/srate*(0:len/2)/len;
fvals=fvals';

G(1:np)=complex(e(2:np+1,2),e(2:np+1,3));
G(np+1:len-np)=1;
G(len-np+1:len)=complex(e(np+2:sz,2),e(np+2:sz,3));


% Reconstructs the original waveform from the measured step
% in the frequency plane using overlap-add method
%
% fname = measured wave filename
% fresp = probe freq. response
% frate = freq. response sample rate
% wi = measured wave
% r = reconstruct wave
% wrate = reconstructed wave rate

function [wi,r,wrate]=wave_restore(fname,fresp,frate)

segsize=1024; %size of sliced segments

wi=csvread(fname);
wrate=abs(wi(1,1)-wi(2,1));
disp(sprintf('Wave file resolution = %d ns',wrate/1e-9));
wi=wi(:,2);

%correct for odd sizes
[sz,a]=size(wi);
if (mod(sz,2)>0)
 wi=wi(1:sz-1);
end;

%interpolate/decimate if needed
ifact=wrate/frate;
if(ifact>1)
 ifact=round(ifact);
end;
if(ifact>1)
 m=interp(wi,round(ifact));
else
 if (ifact<1)
  m=decimate(wi,round(1/ifact));
 else
  m=wi;
 end;
end;
[sz,a]=size(m);

%pad with zeros
[fsize,x]=size(fresp);
wsize=(floor(sz/segsize+1)*segsize);
m(sz:wsize)=0;

%slice in segments and add together
w=zeros(fsize,1);
R=zeros(fsize,1);
r=zeros(fsize,1);
q=zeros(wsize,1);
h=hann(segsize); %Hanning windowing -> segments must overlap 50%

for i=1:wsize/segsize*2-1
 %cut a slice
 w(1:segsize)=m(1+(i-1)*segsize/2:(i+1)*segsize/2);
 %apply a Hanning window
 s=1:segsize;
 w(s)=w(s).*h(s);
 %take the FFT of the slice
 W=fft(w);
 %filter the slice FFT
 s=1:fsize;
 R(s)=W(s)./fresp(s);
 %calculate the IFFT of the filtered slice
 r=ifft(R);
 %sum to the resulting step
 sa=1+(i-1)*segsize/2; %q interval start
 ua=1; %r interval start
 slen=min(wsize-sa+1,fsize); %intervals length
 s=sa:sa+slen-1;
 u=ua:ua+slen-1;
 q(s)=q(s)+r(u);
end;

%fix the beginning (half slice)
w=zeros(fsize,1);
w(1:segsize/2)=m(1:segsize/2);
s=1:segsize/2;
w(s)=w(s).*h(s+segsize/2);
W=fft(w);
s=1:fsize;
R(s)=W(s)./fresp(s);
r=ifft(R);
s=1:min(fsize,wsize);
q(s)=q(s)+r(s);

%fix the end (half slice)
w=zeros(fsize,1);
w(1:segsize/2)=m(wsize-segsize/2+1:wsize);
s=1:segsize/2;
w(s)=w(s).*h(s);
W=fft(w);
s=1:fsize;
R(s)=W(s)./fresp(s);
r=ifft(R);
s=wsize-segsize/2+1:wsize;
q(s)=q(s)+r(1:segsize/2);

%cut to the original length
r=q(1:sz);
r=real(r);

%decimate/interpolate if needeed
if(ifact>1)
 r=decimate(r,round(ifact));
else
 if (ifact<1)
  r=interp(r,round(1/ifact));
 end;
end;