clear all;
close all;
clc;
%% Script of the time constant estimation of signals used in the sponge experiments
%
%% Description
%
% There are several steps to estimate the time constant of signals used in
% the sponge experiments. Firstly, the part of signal representing a step
% response is chosen. Then a maximum and minimum of this part is found out 
% and the function linspace is used to determine a line. The value corresponding 
% to 37% of that line is found (for highpass filters this value have to corresponds to 37%,
% for lowpass filters this value have to be 63%). This script is sensitive to 
% manual setting of borders of the step response area. The time constants
% are estimated for each output separately.
%
%% Loading of signals used in the sponge experiments
%
infoorigsignal = h5info('sponge_C1_bip_01_2015.06.15_13.17.11.hdf5');               
origin = h5read('sponge_C1_bip_01_2015.06.15_13.17.11.hdf5', '/RawData/Samples');
%
%% Definition of the sampling frequency
fs = 4800;
%
%% Plotting of signal to determine the step response area
%
figure
plot(origin(1,:))
hold on
plot(origin(4,:),'g')
%
%% Definition of borders of the step response part of signal
%
response = origin(4,1973983:(1973983+2500));
%
%% Definition of these borders for the simplest signals used in the sponge experiments
% AP response = origin(2, 1897040:1899040);
% LP response = origin(2, 1891000:1900000);
% HP response = origin(2, 1879110:1881110);
%
%% Plotting of the step response to verify 
% that the shape of that response corresponds to highpass/lowpass filter
%
figure
plot(response,'g')
%
%% Estimation of the maximal and minimal value of the step response
%
maximum = max(response);
minimum = min(response);
%
%% Using linspace function to create line with lower limit 
% corresponding to minimum value and upper limit corresponding to maximum
% value
%
y = linspace(minimum,maximum);
%
%% The value corresponding to 37% of step response is chosen
%
y_37 = y(37);
ind11 = find(response>y_37);
ind1 = ind11(end);
ind2 = response(ind1);
%
%% Estimation of the time constant
%
time_constant = ind1*(1/fs);
%
%% Plotting the time constant in figure shown the step response of the system
%
x = [ind1+100 ind1+100];
y = [58400 ind2];
%
x2 = [0 ind1+100];
y2 = [ind2 ind2];
%
figure
plot(origin(4,1973883:(1973883+2500)), 'g')
xlabel('Number of samples')
ylabel('Output')
ylim([58400 maximum])
hold on
line(x,y, 'Color', 'r', 'LineStyle', '--')
hold on
line(x2,y2, 'Color', 'r', 'LineStyle', '--')
legend('output response on the step function','x value of a time constant', 'y value of a time constant')
