function ExtractFeatures(dir_name, file_struct, f_subsample, threshold_amp, threshold_dist, width_th, verbose)
%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This function acts on the data extracted using the ExtractAbfData 
% function. It filters the subsampled data with a high pass filter in order
% to further use in the detection and characterisation of depolarisation 
% events.
%
% dir_name = a string defining the path of the folder where the input data 
% is and the output data should be saved.
% file_struct = file structure with the datasets from the desired session:
%
% file_struct(file_number, feature_number) = [session_number recording_number
% baseline/response control/target dpf errors]
% 
% Note that:
% session_number = number of the session - should be present in the
% folder name;
% recording_number = number of the recording - should be present in the
% last digits of the recording name;
% baseline/response = 0 if the trace is recorded only during spontaneous
% activity; 1 if trace contains response to PTZ 
% control/target = 0 if the trace is recorded from a Mismatch control fish;
% 1 if trace is recorded from a DEPDC5 animal
%
% f_subsample [Hz] = the sampling frequency of the signal
% threshold_amp [mV] = the amplitude threshold for event detection
% threshold_dist [s] = threshold for inter-event distance in order to be 
% considered separate depolarization events
% width_th [s] = threshold for event length in order to be considered a 
% depolarization event
% verbose = flag to select wether to plot intermediate computations (1) or 
% not (0)
%
%
% output: saves a events and features in peak_matrix in with the following 
% format:
% peak_matrix (events, features) = [event_start, event_end, event_height,
% event_width, data_length] 
%
% Note that:
% event_start [samples] = timestamp where the event starts;
% event_end [samples] = timestamp where the event ends;
% event_height [mV] = amplitude of event;
% event_width [samples] = duration of event;
% data_length [samples] = length of dataset the event originates in;
%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% convert threshold_dist and width_th from seconds to samples
threshold_dist_idx = threshold_dist*f_subsample;
width_th_idx = width_th*f_subsample;

% find a list of files in the source folder
list = dir(dir_name);
if length(list)>2
    filename_base = list(3,1).name;
    if file_struct < 0
        file_struct = 0:length(list)-3;
    end
    clear list;
    
    % iterate through datasets
    for i_datasets=1:length(file_struct)
        
        % construct filename
        if file_struct(i_datasets)<10
            filename = [filename_base(1:7) num2str(file_struct(i_datasets))];
        else
            filename = [filename_base(1:6) num2str(file_struct(i_datasets))];
        end
      
        if exist([dir_name '\' filename '_filt.mat' ], 'file')
            % load file containing filtered traces: ANALYSIS and DETECT signals
            load([dir_name '\' filename '_filt' ]);

            % find indices where amplitude is higher then the threshold
            index = find(abs(V_filt_extract)>threshold_amp);
            
            if (index)
                % define event start and end as the begining and end of
                % threshold crossings in the DETECTION signal. 
                dif_index = diff(index);
                hist = find(dif_index>threshold_dist_idx);
                ev_start = [max(index(1)-1, 1); index(hist+1)-1];
                ev_end = [index(hist)+1; min(length(V_filt_extract),index(end)+1)];
                
                % compute event features: amplitude, width;
                for i_events=1:length(ev_start)
                    ev_max = max(V_filt_analyze(ev_start(i_events):ev_end(i_events)));
                    ev_min = min(V_filt_analyze(ev_start(i_events):ev_end(i_events))); 
                    ev_avg = mean(V_filt_analyze(ev_start(i_events)-10:ev_start(i_events)-1));
                    if abs(ev_avg-ev_max)>abs(ev_avg-ev_min)
                        ev_sgn = -1;
                    else
                        ev_sgn = 1;
                    end
                    ev_height(i_events) = ev_sgn*(ev_max-ev_min);  
                end
                ev_width = ev_end-ev_start;
                
                % store events and features in a matrix
                data_length = ones(length(ev_width),1)*length(V_filt_extract);
                peak_matrix = [ev_start, ev_end, ev_height', ev_width, data_length];
                peak_matrix(find(peak_matrix(:,4)<=width_th_idx),:) = [];
                   
                % plotting (if verbose active)
                if verbose
                    t=0:1/f_subsample:(length(V_filt_extract)-1)/f_subsample;

                    figure;
                    subplot(2,1,1);
                    plot(t, V_filt_extract);
                    hold on;
                    plot((peak_matrix(:,1)-1)/f_subsample, V_filt_extract(peak_matrix(:,1)), '*g');
                    plot((peak_matrix(:,2)-1)/f_subsample, V_filt_extract(peak_matrix(:,2)), '*k');

                    subplot(2,1,2);
                    plot(t, V_filt_analyze);
                    hold on;
                    plot((peak_matrix(:,1)-1)/f_subsample, V_filt_analyze(peak_matrix(:,1)), '*g');
                    plot((peak_matrix(:,2)-1)/f_subsample, V_filt_analyze(peak_matrix(:,2)), '*k');
                    title(filename);
                end
            else
                peak_matrix = [];
            end
            save([dir_name '\' filename '_peaks' ], 'peak_matrix');
            clear peak_matrix ev_start ev_end ev_height ev_width index;
        end
    end
end


