%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Water probe calibration using a falling step
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear;
clf;

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

[m_in,res]=step_rd(strcat(fin,'.csv'),0);

m_in=decimate(m_in,5); %tentative
res=res*5; %tentative

m_ok=step_ideal(m_in);
m_corr=diagonalize(m_in,65536); %tentative
m_ok=diagonalize(m_ok,65536); %tentative
[len,x]=size(m_corr);
lenh=len/2+1;
f=1/res*(0:len/2)/len;
f=f';

%plot measured vs. ideal step
figure(1);
plot(1:len,m_corr,1:len,m_ok);
xlabel('Time');
ylabel('Amplitude');
legend('measured','ideal',0);
title ('CALIBRATION STEP');

Y=fft(m_corr,len);
Yok=fft(m_ok,len);

%plot measured vs. ideal FFT
figure(2);
loglog(f,abs(Y(1:lenh)),f,abs(Yok(1:lenh)));
xlabel('Freq.');
ylabel('FFT');
set(gca,'XLim',[1e4 1e8]);
legend('measured','ideal',0);
title('STEP FFT');

[G,npoints]=calc_fresp(Y,Yok,res);

%plot module
figure(1);
semilogx(f,abs(G(1:lenh)),f,1);
xlabel('Freq.');
ylabel('mod(G)');
set(gca,'XLim',[1e4 1e8]);
title('PROBE FREQUENCY RESPONSE');
%plot phase
figure(2);
semilogx(f,angle(G(1:lenh))*180/3.1416);
xlabel('Freq.');
ylabel('phase(G)');
set(gca,'XLim',[1e4 1e8]);
title('PROBE FREQUENCY RESPONSE');

%save the probe response
save_fresp(strcat(fin,'_G.csv'),f,G,npoints,res);


% Calculates the ideal falling step
% from a sampled one
%
% m = sampled step array

function id=step_ideal(m)

[len,x]=size(m);
f=find_edge(m);
a=mean(m(1:f-200));
b=mean(m(f+2000:len));
id(1:f)=a;
id(f+1:len)=b;
id=id';


% Finds the position of the step falling edge
%
% m = step array

function p=find_edge(m)

[len,a]=size(m);
deltam=0;
p=1;
for i=1:len-1
 delta=abs(m(i)-m(i+1));
 if delta>deltam
  p=i;
  deltam=delta;
 end;
end;


% Removes a diagonal from a sampled step
%
% m = sampled step
% s = desired size (if > 0)

function r=diagonalize(m,s)

[len,x]=size(m);
f=find_edge(m);
a=mean(m(1:f-200));
b=mean(m(f+2000:len));

%make the diagonal and subtract it
k=(b-a)/(len-1);
i=1:len;
diag(i)=a+k*(i-1);

r=m(:)-diag(:);

%pad with "zeroes"
if s>0
 r(len+1:s)=0;
end;


% Calculates the probe freq. response
%
% meas = measured FFT
% ideal = ideal FFT
% res = sampling resolutions (timebase)
% resp =freq. response
% np = n. of valid points (<>1)

function [resp,np]=calc_fresp(meas,ideal,res)

[len,x]=size(meas);
resp=zeros(1,len);
for i=1:len
 if (ideal(i)==0) & (meas(i)==0)
  resp(i)=resp(i-1);
 elseif (ideal(i)==0)
  resp(i)=0;
 else
  resp(i)=meas(i)/ideal(i);
 end;
end;
resp=resp.';

%lowpass filter with cut @ 12 MHz
np=round(12e6*res*len);
resp(np+1:len-np)=1;

return;


% Saves a frequency response to a csv file
%
% fout = output file name
% fvals = frequency values (sized L/2+1)
% G = fresp matrix (sized L)
% np = n. of points to save
% srate = sampling rate

function save_fresp(fout,fvals,G,np,srate)

[len,x]=size(G);
e=zeros(np*2+1,3);
%save sampling rate and original size
e(1,1)=srate;
e(1,2)=len;
%save the first half
i=2:np+1;
j=1:np;
e(i,1)=fvals(j);
e(i,2)=real(G(j));
e(i,3)=imag(G(j));
%save the (mirrored) second half
i=np+2:np*2+1;
j=len-np+1:len;
e(i,2)=real(G(j));
e(i,3)=imag(G(j));

csvwrite(fout,e(:,:));