clear

num_channels = 2; % ADC input info from header file
fileinfo = dir('analogin.dat');
num_samples = fileinfo.bytes/(num_channels * 2); % uint16 = 2 bytes
fid = fopen('analogin.dat', 'r');

v = fread(fid, [num_channels, num_samples], 'uint16');
fclose(fid);
v = v * 0.000050354; % convert to volts
v(2,1:8e4) = 0;
%v(2,1:96e5) = 0;

%uncomment following for saline measurement
%v = v(2:-1:1,:);
%%
fs = 2e4; % Hz
stim_dur = 10; % s
trial_dur = 30; % s
stim_on = find(v(2,:)>2);
start_ind = stim_on(1) - 4*fs;
Nlength = trial_dur *fs;
Ntrial = 30;
vv = zeros(1,Nlength);
time_plot = linspace(0,trial_dur, trial_dur*fs);
vv1 = vv;
for i = 1:Ntrial
    vv = vv + v(1,start_ind+(i-1)*Nlength:start_ind+i*Nlength-1)./(Ntrial);
    vv1 = vv1 +v(2,start_ind+(i-1)*Nlength:start_ind+i*Nlength-1)./(Ntrial);
end

fc=128;% cut off frequency
order = 2;
[b14 a14]=butter(order,(fc/fs),'low');
vv_low=filtfilt(b14,a14,vv);
%plot(vv)
figure(90)
plot(time_plot,vv1)
hold on
plot(time_plot,vv_low);


%
Vmean = mean(vv_low(1*fs:3.5*fs));

Vs = 3.3; %V
dV = vv_low - Vmean;
Beta = 3454;%K ^2
T0 = 310; % K

dT = NTC_temperature(dV,Vmean, Vs, Beta,T0);

figure(99)
plot(time_plot,dT);
hold on


function dT = NTC_temperature(dV,Vmean, Vs, Beta,T0)

VV = (Vs/(Vmean*(Vs-Vmean)) *dV);
dT = - T0^2/Beta * log((1 + VV));

end

