drta_read_Intan_RHD2000_header.m

function [draq_p,draq_d]=drta_read_Intan_RHD2000_header(filename, which_protocol, handles)

% read_Intan_RHD2000_file
%
% Version 2.0, 20 October 2016
%
% Reads Intan Technologies RHD2000 data file generated by evaluation board
% GUI or Intan Recording Controller.  Data are parsed and placed into
% variables that appear in the base MATLAB workspace.  Therefore, it is
% recommended to execute a 'clear' command before running this program tort
% clear all other variables from the base workspace.
%
% Example:
% >> clear
% >> read_Intan_RHD2000_file
% >> whos
% >> amplifier_channels(1)
% >> plot(t_amplifier, amplifier_data(1,:))

% [file, path, filterindex] = ...
%     uigetfile('*.rhd', 'Select an RHD2000 Data File', 'MultiSelect', 'off')
% 
% if (file == 0)
%     return;
% end

% Read most recent file automatically.
% path = 'C:\Users\Reid\Documents\RHD2132\testing\';
% d = dir([path '*.rhd']);
% file = d(end).name;

tic;
fid = fopen(filename, 'r');

s = dir(filename);
filesize = s.bytes;

% Check 'magic number' at beginning of file to make sure this is an Intan
% Technologies RHD2000 data file.
magic_number = fread(fid, 1, 'uint32');
if magic_number ~= hex2dec('c6912702')
    error('Unrecognized file type.');
end

% Read version number.
data_file_main_version_number = fread(fid, 1, 'int16');
data_file_secondary_version_number = fread(fid, 1, 'int16');

fprintf(1, '\n');
fprintf(1, 'Reading Intan Technologies RHD2000 Data File, Version %d.%d\n', ...
    data_file_main_version_number, data_file_secondary_version_number);
fprintf(1, '\n');

if (data_file_main_version_number == 1)
    num_samples_per_data_block = 60;
else
    num_samples_per_data_block = 128;
end

% Read information of sampling rate and amplifier frequency settings.
sample_rate = fread(fid, 1, 'single');
dsp_enabled = fread(fid, 1, 'int16');
actual_dsp_cutoff_frequency = fread(fid, 1, 'single');
actual_lower_bandwidth = fread(fid, 1, 'single')
actual_upper_bandwidth = fread(fid, 1, 'single')

desired_dsp_cutoff_frequency = fread(fid, 1, 'single');
desired_lower_bandwidth = fread(fid, 1, 'single');
desired_upper_bandwidth = fread(fid, 1, 'single');

fprintf(1, 'Sample rate = %d\nLower bandwidth = %d\nUpper_bandwidth = %d\n', ...
    sample_rate, actual_lower_bandwidth,actual_upper_bandwidth);

% This tells us if a software 50/60 Hz notch filter was enabled during
% the data acquisition.
notch_filter_mode = fread(fid, 1, 'int16');
notch_filter_frequency = 0;
if (notch_filter_mode == 1)
    notch_filter_frequency = 50;
elseif (notch_filter_mode == 2)
    notch_filter_frequency = 60;
end

desired_impedance_test_frequency = fread(fid, 1, 'single');
actual_impedance_test_frequency = fread(fid, 1, 'single');

% Place notes in data strucure
notes = struct( ...
    'note1', fread_QString(fid), ...
    'note2', fread_QString(fid), ...
    'note3', fread_QString(fid) );
    
% If data file is from GUI v1.1 or later, see if temperature sensor data
% was saved.
num_temp_sensor_channels = 0;
if ((data_file_main_version_number == 1 && data_file_secondary_version_number >= 1) ...
    || (data_file_main_version_number > 1))
    num_temp_sensor_channels = fread(fid, 1, 'int16');
end

% If data file is from GUI v1.3 or later, load eval board mode.
eval_board_mode = 0;
if ((data_file_main_version_number == 1 && data_file_secondary_version_number >= 3) ...
    || (data_file_main_version_number > 1))
    eval_board_mode = fread(fid, 1, 'int16');
end

% If data file is from v2.0 or later (Intan Recording Controller),
% load name of digital reference channel.
if (data_file_main_version_number > 1)
    reference_channel = fread_QString(fid);
end

% Place frequency-related information in data structure.
frequency_parameters = struct( ...
    'amplifier_sample_rate', sample_rate, ...
    'aux_input_sample_rate', sample_rate / 4, ...
    'supply_voltage_sample_rate', sample_rate / num_samples_per_data_block, ...
    'board_adc_sample_rate', sample_rate, ...
    'board_dig_in_sample_rate', sample_rate, ...
    'desired_dsp_cutoff_frequency', desired_dsp_cutoff_frequency, ...
    'actual_dsp_cutoff_frequency', actual_dsp_cutoff_frequency, ...
    'dsp_enabled', dsp_enabled, ...
    'desired_lower_bandwidth', desired_lower_bandwidth, ...
    'actual_lower_bandwidth', actual_lower_bandwidth, ...
    'desired_upper_bandwidth', desired_upper_bandwidth, ...
    'actual_upper_bandwidth', actual_upper_bandwidth, ...
    'notch_filter_frequency', notch_filter_frequency, ...
    'desired_impedance_test_frequency', desired_impedance_test_frequency, ...
    'actual_impedance_test_frequency', actual_impedance_test_frequency );

% Define data structure for spike trigger settings.
spike_trigger_struct = struct( ...
    'voltage_trigger_mode', {}, ...
    'voltage_threshold', {}, ...
    'digital_trigger_channel', {}, ...
    'digital_edge_polarity', {} );

new_trigger_channel = struct(spike_trigger_struct);
spike_triggers = struct(spike_trigger_struct);

% Define data structure for data channels.
channel_struct = struct( ...
    'native_channel_name', {}, ...
    'custom_channel_name', {}, ...
    'native_order', {}, ...
    'custom_order', {}, ...
    'board_stream', {}, ...
    'chip_channel', {}, ...
    'port_name', {}, ...
    'port_prefix', {}, ...
    'port_number', {}, ...
    'electrode_impedance_magnitude', {}, ...
    'electrode_impedance_phase', {} );

new_channel = struct(channel_struct);

% Create structure arrays for each type of data channel.
amplifier_channels = struct(channel_struct);
aux_input_channels = struct(channel_struct);
supply_voltage_channels = struct(channel_struct);
board_adc_channels = struct(channel_struct);
board_dig_in_channels = struct(channel_struct);
board_dig_out_channels = struct(channel_struct);

amplifier_index = 1;
aux_input_index = 1;
supply_voltage_index = 1;
board_adc_index = 1;
board_dig_in_index = 1;
board_dig_out_index = 1;

% Read signal summary from data file header.

number_of_signal_groups = fread(fid, 1, 'int16');

for signal_group = 1:number_of_signal_groups
    signal_group_name = fread_QString(fid);
    signal_group_prefix = fread_QString(fid);
    signal_group_enabled = fread(fid, 1, 'int16');
    signal_group_num_channels = fread(fid, 1, 'int16');
    signal_group_num_amp_channels = fread(fid, 1, 'int16');
  
    if (signal_group_num_channels > 0 && signal_group_enabled > 0)
        new_channel(1).port_name = signal_group_name;
        new_channel(1).port_prefix = signal_group_prefix;
        new_channel(1).port_number = signal_group;
        for signal_channel = 1:signal_group_num_channels
            new_channel(1).native_channel_name = fread_QString(fid);
            new_channel(1).custom_channel_name = fread_QString(fid);
            new_channel(1).native_order = fread(fid, 1, 'int16');
            new_channel(1).custom_order = fread(fid, 1, 'int16');
            signal_type = fread(fid, 1, 'int16');
            channel_enabled = fread(fid, 1, 'int16');
            new_channel(1).chip_channel = fread(fid, 1, 'int16');
            new_channel(1).board_stream = fread(fid, 1, 'int16');
            new_trigger_channel(1).voltage_trigger_mode = fread(fid, 1, 'int16');
            new_trigger_channel(1).voltage_threshold = fread(fid, 1, 'int16');
            new_trigger_channel(1).digital_trigger_channel = fread(fid, 1, 'int16');
            new_trigger_channel(1).digital_edge_polarity = fread(fid, 1, 'int16');
            new_channel(1).electrode_impedance_magnitude = fread(fid, 1, 'single');
            new_channel(1).electrode_impedance_phase = fread(fid, 1, 'single');
             
            if (channel_enabled)
                switch (signal_type)
                    case 0
                        amplifier_channels(amplifier_index) = new_channel;
                        spike_triggers(amplifier_index) = new_trigger_channel;
                        amplifier_index = amplifier_index + 1;
                    case 1
                        aux_input_channels(aux_input_index) = new_channel;
                        aux_input_index = aux_input_index + 1;
                    case 2
                        supply_voltage_channels(supply_voltage_index) = new_channel;
                        supply_voltage_index = supply_voltage_index + 1;
                    case 3
                        board_adc_channels(board_adc_index) = new_channel;
                        board_adc_index = board_adc_index + 1;
                    case 4
                        board_dig_in_channels(board_dig_in_index) = new_channel;
                        board_dig_in_index = board_dig_in_index + 1;
                    case 5
                        board_dig_out_channels(board_dig_out_index) = new_channel;
                        board_dig_out_index = board_dig_out_index + 1;
                    otherwise
                        error('Unknown channel type');
                end
            end
            
        end
    end
end

% Summarize contents of data file.
num_amplifier_channels = amplifier_index - 1;
num_aux_input_channels = aux_input_index - 1;
num_supply_voltage_channels = supply_voltage_index - 1;
num_board_adc_channels = board_adc_index - 1;
num_board_dig_in_channels = board_dig_in_index - 1;
num_board_dig_out_channels = board_dig_out_index - 1;
% 
% fprintf(1, 'Found %d amplifier channel%s.\n', ...
%     num_amplifier_channels, plural(num_amplifier_channels));
% fprintf(1, 'Found %d auxiliary input channel%s.\n', ...
%     num_aux_input_channels, plural(num_aux_input_channels));
% fprintf(1, 'Found %d supply voltage channel%s.\n', ...
%     num_supply_voltage_channels, plural(num_supply_voltage_channels));
% fprintf(1, 'Found %d board ADC channel%s.\n', ...
%     num_board_adc_channels, plural(num_board_adc_channels));
% fprintf(1, 'Found %d board digital input channel%s.\n', ...
%     num_board_dig_in_channels, plural(num_board_dig_in_channels));
% fprintf(1, 'Found %d board digital output channel%s.\n', ...
%     num_board_dig_out_channels, plural(num_board_dig_out_channels));
% fprintf(1, 'Found %d temperature sensors channel%s.\n', ...
%     num_temp_sensor_channels, plural(num_temp_sensor_channels));
% fprintf(1, '\n');

% Determine how many samples the data file contains.

% Each data block contains num_samples_per_data_block amplifier samples.
bytes_per_block = num_samples_per_data_block * 4;  % timestamp data
bytes_per_block = bytes_per_block + num_samples_per_data_block * 2 * num_amplifier_channels;
% Auxiliary inputs are sampled 4x slower than amplifiers
bytes_per_block = bytes_per_block + (num_samples_per_data_block / 4) * 2 * num_aux_input_channels;
% Supply voltage is sampled once per data block
bytes_per_block = bytes_per_block + 1 * 2 * num_supply_voltage_channels;
% Board analog inputs are sampled at same rate as amplifiers
bytes_per_block = bytes_per_block + num_samples_per_data_block * 2 * num_board_adc_channels;
% Board digital inputs are sampled at same rate as amplifiers
if (num_board_dig_in_channels > 0)
    bytes_per_block = bytes_per_block + num_samples_per_data_block * 2;
end
% Board digital outputs are sampled at same rate as amplifiers
if (num_board_dig_out_channels > 0)
    bytes_per_block = bytes_per_block + num_samples_per_data_block * 2;
end
% Temp sensor is sampled once per data block
if (num_temp_sensor_channels > 0)
   bytes_per_block = bytes_per_block + 1 * 2 * num_temp_sensor_channels; 
end

% How many data blocks remain in this file?
data_present = 0;
bytes_remaining = filesize - ftell(fid);
if (bytes_remaining > 0)
    data_present = 1;
end

num_data_blocks = bytes_remaining / bytes_per_block;

num_amplifier_samples = num_samples_per_data_block * num_data_blocks;
num_aux_input_samples = (num_samples_per_data_block / 4) * num_data_blocks;
num_supply_voltage_samples = 1 * num_data_blocks;
num_board_adc_samples = num_samples_per_data_block * num_data_blocks;
num_board_dig_in_samples = num_samples_per_data_block * num_data_blocks;
num_board_dig_out_samples = num_samples_per_data_block * num_data_blocks;

record_time = num_amplifier_samples / sample_rate;

% if (data_present)
%     fprintf(1, 'File contains %0.3f seconds of data.  Amplifiers were sampled at %0.2f kS/s.\n', ...
%         record_time, sample_rate / 1000);
%     fprintf(1, '\n');
% else
%     fprintf(1, 'Header file contains no data.  Amplifiers were sampled at %0.2f kS/s.\n', ...
%         sample_rate / 1000);
%     fprintf(1, '\n');
% end

if (data_present)
    
    % Pre-allocate memory for data.
%     fprintf(1, 'Allocating memory for data...\n');

    %The following are read, but are not stored
%     t_amplifier = zeros(1, num_amplifier_samples);
%     amplifier_data = zeros(num_amplifier_channels, num_amplifier_samples);
%     aux_input_data = zeros(num_aux_input_channels, num_aux_input_samples);
%     supply_voltage_data = zeros(num_supply_voltage_channels, num_supply_voltage_samples);
%     temp_sensor_data = zeros(num_temp_sensor_channels, num_supply_voltage_samples);
%     board_adc_data = zeros(num_board_adc_channels, num_board_adc_samples);
%     board_dig_out_data = zeros(num_board_dig_out_channels, num_board_dig_out_samples);
%     board_dig_out_raw = zeros(1, num_board_dig_out_samples);
    
    amplifier_data = zeros(num_amplifier_channels, num_samples_per_data_block);
    aux_input_data = zeros(num_aux_input_channels, num_samples_per_data_block);
    supply_voltage_data = zeros(num_supply_voltage_channels, 1);
    temp_sensor_data = zeros(num_temp_sensor_channels, 1); 
    board_adc_data = zeros(num_board_adc_channels, num_samples_per_data_block);
    board_dig_out_data = zeros(num_board_dig_out_channels, int64(num_board_dig_out_samples));
    board_dig_out_raw = zeros(1, num_samples_per_data_block);
    
    %We do read in all the t_amplifier and dig data to allow determining the time gaps and the trials
    t_amplifier = zeros(1, int64(num_amplifier_samples));
     
    switch which_protocol
        case {1,3,5,6,7,8}
        board_dig_in_raw = zeros(1, int64(num_board_dig_in_samples));
    end
    board_dig_in_data = zeros(num_board_dig_in_channels, int64(num_board_dig_in_samples));

    % Read sampled data from file.
%     fprintf(1, 'Reading data from file...\n');

    amplifier_index = 1;
    aux_input_index = 1;
    supply_voltage_index = 1;
    board_adc_index = 1;
    board_dig_in_index = 1;
    board_dig_out_index = 1;

    print_increment = 10;
    percent_done = print_increment;
    
    for i=1:num_data_blocks
        % In version 1.2, we moved from saving timestamps as unsigned
        % integeters to signed integers to accomidate negative (adjusted)
        % timestamps for pretrigger data.
        
        if ((data_file_main_version_number == 1 && data_file_secondary_version_number >= 2) ...
                || (data_file_main_version_number > 1))
            t_amplifier(amplifier_index:(amplifier_index + num_samples_per_data_block - 1)) = fread(fid, num_samples_per_data_block, 'int32');
        else
            t_amplifier(amplifier_index:(amplifier_index + num_samples_per_data_block - 1)) = fread(fid, num_samples_per_data_block, 'uint32');
        end
        draq_d.offset_start_ch(i)=ftell(fid);
        if (num_amplifier_channels > 0)
            %amplifier_data(:, amplifier_index:(amplifier_index + num_samples_per_data_block - 1)) = fread(fid, [num_samples_per_data_block, num_amplifier_channels], 'uint16')';
            amplifier_data = fread(fid, [num_samples_per_data_block, num_amplifier_channels], 'uint16')';
        end
        draq_d.offset_end_ch(i)=ftell(fid);
        if (num_aux_input_channels > 0)
            aux_input_data = fread(fid, [(num_samples_per_data_block / 4), num_aux_input_channels], 'uint16')';
            %aux_input_data(:, aux_input_index:(aux_input_index + (num_samples_per_data_block / 4) - 1)) = fread(fid, [(num_samples_per_data_block / 4), num_aux_input_channels], 'uint16')';
        end
        if (num_supply_voltage_channels > 0)
            supply_voltage_data = fread(fid, [1, num_supply_voltage_channels], 'uint16')';
        end
        if (num_temp_sensor_channels > 0)
            temp_sensor_data = fread(fid, [1, num_temp_sensor_channels], 'int16')';
        end
        draq_d.offset_start_adc(i)=ftell(fid);
        if (num_board_adc_channels > 0)
            switch which_protocol
                case {1,3,5,6,7,8}
                    board_adc_data = fread(fid, [num_samples_per_data_block, num_board_adc_channels], 'uint16')';
                case {2,4}
                    board_adc_data(:, board_adc_index:(board_adc_index + num_samples_per_data_block - 1)) = fread(fid, [num_samples_per_data_block, num_board_adc_channels], 'uint16')';
            end
        end
        draq_d.offset_end_adc(i)=ftell(fid);
        draq_d.offset_start_dig(i)=ftell(fid);
        if (num_board_dig_in_channels > 0)
            switch which_protocol
                case {1,3,5,6,7,8}
                board_dig_in_raw(board_dig_in_index:(board_dig_in_index + num_samples_per_data_block - 1)) = fread(fid, num_samples_per_data_block, 'uint16');
                case {2,4}
                board_dig_in_raw = fread(fid, num_samples_per_data_block, 'uint16');
            end
        end
        draq_d.offset_end_dig(i)=ftell(fid);
        if (num_board_dig_out_channels > 0)
            board_dig_out_raw = fread(fid, num_samples_per_data_block, 'uint16')';
        end

        amplifier_index = amplifier_index + num_samples_per_data_block;
        aux_input_index = aux_input_index + (num_samples_per_data_block / 4);
        supply_voltage_index = supply_voltage_index + 1;
        board_adc_index = board_adc_index + num_samples_per_data_block;
        board_dig_in_index = board_dig_in_index + num_samples_per_data_block;
        board_dig_out_index = board_dig_out_index + num_samples_per_data_block;

        fraction_done = 100 * (i / num_data_blocks);
        if (fraction_done >= percent_done)
            fprintf(1, '%d%% done...\n', percent_done);
            percent_done = percent_done + print_increment;
        end
    end

    % Make sure we have read exactly the right amount of data.
    bytes_remaining = filesize - ftell(fid);
    if (bytes_remaining ~= 0)
        %error('Error: End of file not reached.');
    end

end

% Close data file.
fclose(fid);

if (data_present)
    
    fprintf(1, 'Parsing data...\n');
    
    switch which_protocol
        case {1,3,5,6,7,8}
            % Extract digital input channels to separate variables.
            for i=1:num_board_dig_in_channels
                mask = 2^(board_dig_in_channels(i).native_order) * ones(size(board_dig_in_raw));
                board_dig_in_data(i, :) = (bitand(board_dig_in_raw, mask) > 0);
            end
            try
                for i=1:num_board_dig_out_channels
                    mask = 2^(board_dig_out_channels(i).native_order) * ones(size(board_dig_out_raw));
                    board_dig_out_data(i, :) = (bitand(board_dig_out_raw, mask) > 0);
                end
            catch
            end
    end

    % Scale voltage levels appropriately.
    amplifier_data = 0.195 * (amplifier_data - 32768); % units = microvolts
    aux_input_data = 37.4e-6 * aux_input_data; % units = volts
    supply_voltage_data = 74.8e-6 * supply_voltage_data; % units = volts
    
    if (eval_board_mode == 1)
        board_adc_data = 152.59e-6 * (board_adc_data - 32768); % units = volts
    elseif (eval_board_mode == 13) % Intan Recording Controller
        board_adc_data = 312.5e-6 * (board_adc_data - 32768); % units = volts    
    else
        board_adc_data = 50.354e-6 * board_adc_data; % units = volts
    end
    temp_sensor_data = temp_sensor_data / 100; % units = deg C

    % Check for gaps in timestamps.
    num_gaps = sum(diff(t_amplifier) ~= 1);
    if (num_gaps == 0)
        fprintf(1, 'No missing timestamps in data.\n');
    else
        fprintf(1, 'Warning: %d gaps in timestamp data found.  Time scale will not be uniform!\n', ...
            num_gaps);
    end

    % Scale time steps (units = seconds).
    t_amplifier = t_amplifier / sample_rate;
    t_aux_input = t_amplifier(1:4:end);
    t_supply_voltage = t_amplifier(1:num_samples_per_data_block:end);
    t_board_adc = t_amplifier;
    t_dig = t_amplifier;
    t_temp_sensor = t_supply_voltage;



end

fprintf(1, 'Finding full trials...\n');

draq_p.dgordra=3;  %3 is rhd
draq_p.show_plot=0;
draq_p.ActualRate=frequency_parameters.board_adc_sample_rate;
draq_p.srate=frequency_parameters.board_adc_sample_rate;
switch which_protocol
    case {1,5,6,7}
        %dropcspm and dropc_conc
%         draq_p.sec_before_trigger=6;
%         draq_p.sec_per_trigger=9;  
        draq_p.sec_before_trigger=handles.pre_dt;
        draq_p.sec_per_trigger=handles.trial_duration;
    case 2
        %laser(Ming)
        draq_p.sec_before_trigger=2;
        draq_p.sec_per_trigger=9;
    case 3
        %dropcnsampler
        draq_p.sec_before_trigger=handles.pre_dt;
        draq_p.sec_per_trigger=handles.trial_duration;
    case 4
        %laser(Merouann)
        draq_p.sec_before_trigger=10;
        draq_p.sec_per_trigger=21;
    case 8
        %working memory
        draq_p.sec_before_trigger=2;
        draq_p.sec_per_trigger=13.5;
end

draq_p.pre_gain=0;
draq_p.scaling=1;
draq_p.offset=0;

draq_d.noTrials=0;
draq_d.trig=1;
draq_d.num_amplifier_channels=num_amplifier_channels;
draq_d.num_samples_per_data_block=num_samples_per_data_block;
draq_d.num_board_adc_channels=num_board_adc_channels;
draq_d.num_board_dig_in_channels=num_board_dig_in_channels;
draq_d.eval_board_mode=eval_board_mode;
draq_d.board_dig_in_channels=board_dig_in_channels;

 

if ~isempty(board_dig_in_data)
    
    switch which_protocol
        case {1,3,6}
            digital_input=board_dig_in_data(1,:)+2*board_dig_in_data(2,:)+4*board_dig_in_data(3,:)...
                +8*board_dig_in_data(4,:)+16*board_dig_in_data(5,:)+32*board_dig_in_data(6,:)...
                +64*board_dig_in_data(7,:);
            digital_input_ns=board_dig_in_data(1,:)+2*board_dig_in_data(2,:)+4*board_dig_in_data(3,:)...
                +8*board_dig_in_data(4,:)+16*board_dig_in_data(5,:);
            digital_input_no2=2*board_dig_in_data(2,:)+4*board_dig_in_data(3,:)...
                +8*board_dig_in_data(4,:)+16*board_dig_in_data(5,:)+32*board_dig_in_data(6,:);
            
        otherwise
            for ii=1:num_board_dig_in_channels
                if ii==1
                    digital_input=board_dig_in_data(1,:);
                else
                    digital_input=digital_input+(2^(ii-1))*board_dig_in_data(ii,:);
                end
            end
            
            if num_board_dig_in_channels>=2
                for ii=2:num_board_dig_in_channels
                    if ii==2
                        digital_input_no2=2*board_dig_in_data(2,:);
                    else
                        digital_input=digital_input+(2^(ii-1))*board_dig_in_data(ii,:);
                    end
                end
                %         digital_input_no2=2*board_dig_in_data(2,:)+4*board_dig_in_data(3,:)...
                %             +8*board_dig_in_data(4,:)+16*board_dig_in_data(5,:)+32*board_dig_in_data(6,:);
            end
    end
else
    digital_input=[];
end
 
%Find the trials
at_end=0;
ii=1;
trials_to_sort=[];
full_trial_start=[];
full_trial_end=[];
 
    
switch which_protocol
    case {1,5,6,8}
        %Find the full trials (excluding short trials)
        while at_end==0
            
            %Find trial (skip changes in the signal that are not actually bonified
            %trials
            found_bonified=0;
            switch which_protocol
                case 1
                    %The trial starts with an output of 6
                    delta_ii_first=find(digital_input_no2(ii:end)==6,1,'first');
                    if ~isempty(delta_ii_first)
                        %The final valve on step has to last at least 0.9 sec
                        ii_step=find(digital_input_no2(ii+delta_ii_first-1:end)>6,1,'first');
                        if ~isempty(ii_step)
                            dt_step1=ii_step/draq_p.ActualRate;
                            
                            if (dt_step1>0.9)&(digital_input_no2(ii+delta_ii_first-1+ii_step-1)==18)
                                %Is there an odor on for >2.4 sec and <4 sec
                                ii_odor_on=find(digital_input_no2(ii+delta_ii_first-1+ii_step-1:end)<18,1,'first');
                                if ~isempty(ii_odor_on)
                                    if ~(((ii_odor_on/draq_p.ActualRate)<2.4)||((ii_odor_on/draq_p.ActualRate)>4))
                                        found_bonified=1;
                                    end
                                end
                            end
                        end
                    end
                     
                case 5
                    %The trial starts with an output of 1
                    delta_ii_first=find(digital_input(ii:end)==1,1,'first');
                    if ~isempty(delta_ii_first)
                        %The final valve on step has to last at least 0.9 sec
                        ii_step=find(digital_input(ii+delta_ii_first-1:end)>1,1,'first');
                        if ~isempty(ii_step)
                            dt_step1=ii_step/draq_p.ActualRate;
                            
                            if (dt_step1>0.9)&(digital_input(ii+delta_ii_first-1+ii_step-1)>=2)&(digital_input(ii+delta_ii_first-1+ii_step-1)<=7)
                                %Is there an odor on for >2.4 sec and <4 sec
                                ii_odor_on=find(digital_input(ii+delta_ii_first-1+ii_step-1:end)~=digital_input(ii+delta_ii_first-1+ii_step-1),1,'first');
                                if ~isempty(ii_odor_on)
                                    if ~(((ii_odor_on/draq_p.ActualRate)<2.4)||((ii_odor_on/draq_p.ActualRate)>4))
                                        found_bonified=1;
                                    end
                                end
                            end
                        end
                    end
                     
                case 6
                    %dropcspm_hf 
                     %The trial starts with an output of 6
                    delta_ii_first=find(digital_input_no2(ii:end)==6,1,'first');
                    if ~isempty(delta_ii_first)
                        %The final valve on step has to last at least 0.5 sec
                        ii_step=find(digital_input_no2(ii+delta_ii_first-1:end)>6,1,'first');
                        if ~isempty(ii_step)
                            dt_step1=ii_step/draq_p.ActualRate;
                            
                            if (dt_step1>0.5)&(digital_input_no2(ii+delta_ii_first-1+ii_step-1)==18)
                                %Is there an odor on for >2.4 sec and <5.5 sec
                                ii_odor_on=find(digital_input_no2(ii+delta_ii_first-1+ii_step-1:end)<18,1,'first');
                                if ~isempty(ii_odor_on)
                                    if (((ii_odor_on/draq_p.ActualRate)>1)||((ii_odor_on/draq_p.ActualRate)<5.5))
                                        found_bonified=1;
                                    end
                                end
                            end
                        end
                    end
                    
                case 8 
                    %Working memory
                    %The trial starts with an output of 1
                    shift_digital_input=bitand(digital_input,1+2+4+8+16);
                    delta_ii_first=find(shift_digital_input(ii:end)==1,1,'first');
                    if ~isempty(delta_ii_first)
                        %The final valve on step has to last at least 0.9 sec
                        ii_step=find(shift_digital_input(ii+delta_ii_first-1:end)>1,1,'first');
                        if ~isempty(ii_step)
                            dt_step1=ii_step/draq_p.ActualRate;
                            
                            if (dt_step1>0.05)&(shift_digital_input(ii+delta_ii_first-1+ii_step-1)>=2)&(shift_digital_input(ii+delta_ii_first-1+ii_step-1)<=7)
                                %Is there an odor on for >2.4 sec and <4 sec
                                ii_odor_on=find(shift_digital_input(ii+delta_ii_first-1+ii_step-1:end)~=shift_digital_input(ii+delta_ii_first-1+ii_step-1),1,'first');
                                if ~isempty(ii_odor_on)
                                    if (ii_odor_on/draq_p.ActualRate)>0.5
                                        found_bonified=1;
                                    end
                                end
                            end
                        end
                    end
            end
            
            
             
            if found_bonified==1
                %Found a bonified digital signal
                %Find the interval
                switch which_protocol
                    case 8
                        ii=ii+delta_ii_first-1;
                        delta_ii_last=(draq_p.sec_per_trigger-draq_p.sec_before_trigger)*draq_p.ActualRate;
                       
                        if ii+delta_ii_last<length(shift_digital_input)
 
                            %Full trial
                            draq_d.noTrials=draq_d.noTrials+1;
                            %Trial start time goes in column 1
                            trials_to_sort(draq_d.noTrials,1)=((ii-draq_p.sec_before_trigger*draq_p.ActualRate)/draq_p.ActualRate);
                            %Block no at start of trial goes in column 2
                            trials_to_sort(draq_d.noTrials,2)=ceil((trials_to_sort(draq_d.noTrials,1)*draq_p.ActualRate)/num_samples_per_data_block);
                            %Block no at end of trial goes in column 3
                            trials_to_sort(draq_d.noTrials,3)=trials_to_sort(draq_d.noTrials,2)+ceil((draq_p.sec_per_trigger*draq_p.ActualRate)/num_samples_per_data_block);
                            full_trial_start(draq_d.noTrials)=ii-draq_p.sec_before_trigger*draq_p.ActualRate;
                            full_trial_end(draq_d.noTrials)=full_trial_start(draq_d.noTrials)+draq_p.sec_per_trigger*draq_p.ActualRate;
                            
%                             figure(11)
%                             plot([1:length(shift_digital_input(floor(trials_to_sort(draq_d.noTrials,1)*draq_p.ActualRate):ceil((trials_to_sort(draq_d.noTrials,1)+draq_p.sec_per_trigger)*draq_p.ActualRate)))]/draq_p.ActualRate...
%                                 ,shift_digital_input(floor(trials_to_sort(draq_d.noTrials,1)*draq_p.ActualRate):ceil((trials_to_sort(draq_d.noTrials,1)+draq_p.sec_per_trigger)*draq_p.ActualRate)))
%                             pffft=1;
                            
                            
                            ii=ii+delta_ii_last;
                        else
                            at_end=1;
                        end
                    otherwise
                        ii=ii+delta_ii_first;
                        delta_ii_last=find(digital_input(ii:end)<1,1,'first');
                        if ~isempty(delta_ii_last)
                            ii=ii+delta_ii_last;
                            if (delta_ii_last/draq_p.ActualRate)>1
                                if found_bonified==1
                                    %Full trial
                                    draq_d.noTrials=draq_d.noTrials+1;
                                    %Trial start time goes in column 1
                                    trials_to_sort(draq_d.noTrials,1)=(ii/draq_p.ActualRate)-(draq_p.sec_before_trigger+0.5);
                                    %Block no at start of trial goes in column 2
                                    trials_to_sort(draq_d.noTrials,2)=ceil((trials_to_sort(draq_d.noTrials,1)*draq_p.ActualRate)/num_samples_per_data_block);
                                    %Block no at end of trial goes in column 3
                                    trials_to_sort(draq_d.noTrials,3)=trials_to_sort(draq_d.noTrials,2)+ceil((draq_p.sec_per_trigger*draq_p.ActualRate)/num_samples_per_data_block);
                                    full_trial_start(draq_d.noTrials)=ii-draq_p.ActualRate*(draq_p.sec_before_trigger+0.5);
                                    full_trial_end(draq_d.noTrials)=full_trial_start(draq_d.noTrials)+draq_p.sec_per_trigger*draq_p.ActualRate;
                                    
                                    %                             figure(11)
                                    %                             plot([1:length(digital_input(floor(trials_to_sort(draq_d.noTrials,1)*draq_p.ActualRate):ceil((trials_to_sort(draq_d.noTrials,1)+draq_p.sec_per_trigger)*draq_p.ActualRate)))]/draq_p.ActualRate...
                                    %                                 ,digital_input(floor(trials_to_sort(draq_d.noTrials,1)*draq_p.ActualRate):ceil((trials_to_sort(draq_d.noTrials,1)+draq_p.sec_per_trigger)*draq_p.ActualRate)))
                                    %                             pffft=1;
                                end
                            end
                        else
                            at_end=1;
                        end
                end
            else
                if isempty(delta_ii_first)
                    at_end=1;
                else
                    ii=ii+delta_ii_first;
                    delta_ii_last=find(digital_input(ii:end)<1,1,'first');
                    if ~isempty(delta_ii_last)
                        ii=ii+delta_ii_last;
                    else
                        at_end=1;
                    end
                end
            end
            
            if isempty(delta_ii_first)
                at_end=1;
            end
        end
        
        full_trials=draq_d.noTrials;
        
        fprintf(1, 'Found %d full trials...\n',full_trials);
        
%         fprintf(1, 'Finding short trials...\n');
%         
%         
%         at_end=0;
%         ii=1;
%         
%         while at_end==0
%             delta_ii_first=find(digital_input(ii:end)>=1,1,'first');
%             %Found a digital signal
%             if ~isempty(delta_ii_first)
%                 %Find the interval
%                 ii=ii+delta_ii_first;
%                 delta_ii_last=find(digital_input(ii:end)<1,1,'first');
%                 if ~isempty(delta_ii_last)
%                     isthere=find((ii-full_trial_start>0)&(ii-full_trial_end<0));
%                     ii=ii+delta_ii_last;
%                     if isempty(isthere)
%                         if (delta_ii_last/draq_p.ActualRate)>1
%                             %Short trial
%                             draq_d.noTrials=draq_d.noTrials+1;
%                             %Trial start time goes in column 1
%                             trials_to_sort(draq_d.noTrials,1)=(ii/draq_p.ActualRate)-(draq_p.sec_before_trigger-1.5);
%                             %Block no at start of trial go in column 2
%                             trials_to_sort(draq_d.noTrials,2)=ceil((trials_to_sort(draq_d.noTrials,1)*draq_p.ActualRate)/num_samples_per_data_block);
%                             %Block no at end of trial go in column 3
%                             trials_to_sort(draq_d.noTrials,3)=trials_to_sort(draq_d.noTrials,2)+ceil((draq_p.sec_per_trigger*draq_p.ActualRate)/num_samples_per_data_block);
%                         end
%                     end
%                 else
%                     at_end=1;
%                 end
%                 
%             else
%                 at_end=1;
%             end
%             figure(10)
%             plot(digital_input(ii-(draq_p.sec_before_trigger+1)*draq_p.ActualRate: ii-(draq_p.sec_before_trigger+1)*draq_p.ActualRate+draq_p.sec_per_trigger*draq_p.ActualRate))
%             pffft=1;
%         end
%         
%         short_trials=draq_d.noTrials-full_trials;
%         fprintf(1, 'Found %d short trials ...\n',short_trials);
        
%         %Find empty trials
%         last_trial=draq_d.noTrials;
%         
%         trials_to_sort=sortrows(trials_to_sort);
%         
%         fprintf(1, 'Finding inter trials...\n');
%         
%         %Empty trial before the first trial
%         if trials_to_sort(1,1)>draq_p.sec_per_trigger+draq_p.sec_before_trigger
%             draq_d.noTrials=draq_d.noTrials+1;
%             trials_to_sort(draq_d.noTrials,1)= trials_to_sort(1,1)-draq_p.sec_per_trigger;
%             trials_to_sort(draq_d.noTrials,2)=ceil((trials_to_sort(draq_d.noTrials,1)*draq_p.ActualRate)/num_samples_per_data_block);
%             trials_to_sort(draq_d.noTrials,3)=trials_to_sort(draq_d.noTrials,2)+ceil((draq_p.sec_per_trigger*draq_p.ActualRate)/num_samples_per_data_block);
%             
% %             figure(10)
% %             %                         plot(digital_input(ii-(draq_p.sec_before_trigger+1)*draq_p.ActualRate: ii-(draq_p.sec_before_trigger+1)*draq_p.ActualRate+draq_p.sec_per_trigger*draq_p.ActualRate))
% %             plot(digital_input(floor(trials_to_sort(draq_d.noTrials,1)*draq_p.ActualRate):ceil((trials_to_sort(draq_d.noTrials,1)+draq_p.sec_per_trigger)*draq_p.ActualRate)))
% %             pffft=1;
%             
%         end
%         
%         %Empty trials between tirals
%         for ii=2:last_trial
%             if (trials_to_sort(ii,1)-trials_to_sort(ii-1,1)-draq_p.sec_per_trigger)>2*draq_p.sec_per_trigger
%                 
%                 %Trial after this trial
%                 if ii~=last_trial
%                     draq_d.noTrials=draq_d.noTrials+1;
%                     trials_to_sort(draq_d.noTrials,1)= trials_to_sort(ii-1,1)+1.2*draq_p.sec_per_trigger;
%                     trials_to_sort(draq_d.noTrials,2)=ceil((trials_to_sort(draq_d.noTrials,1)*draq_p.ActualRate)/num_samples_per_data_block);
%                     trials_to_sort(draq_d.noTrials,3)=trials_to_sort(draq_d.noTrials,2)+ceil((draq_p.sec_per_trigger*draq_p.ActualRate)/num_samples_per_data_block);
%                     
%                     figure(10)
%                     %                         plot(digital_input(ii-(draq_p.sec_before_trigger+1)*draq_p.ActualRate: ii-(draq_p.sec_before_trigger+1)*draq_p.ActualRate+draq_p.sec_per_trigger*draq_p.ActualRate))
%                     plot(digital_input(floor(trials_to_sort(draq_d.noTrials,1)*draq_p.ActualRate):ceil((trials_to_sort(draq_d.noTrials,1)+draq_p.sec_per_trigger)*draq_p.ActualRate)))
%                     pffft=1;
%                         
%                 end
%                 
% %                 %Trial before this trial
% %                 draq_d.noTrials=draq_d.noTrials+1;
% %                 trials_to_sort(draq_d.noTrials,1)= trials_to_sort(ii,1)-draq_p.sec_per_trigger;
% %                 trials_to_sort(draq_d.noTrials,2)=ceil((trials_to_sort(draq_d.noTrials,1)*draq_p.ActualRate)/num_samples_per_data_block);
% %                 trials_to_sort(draq_d.noTrials,3)=trials_to_sort(draq_d.noTrials,2)+ceil((draq_p.sec_per_trigger*draq_p.ActualRate)/num_samples_per_data_block);
% %                 
%             end
%         end
%         
%         fprintf(1, 'Found %d inter trials ...\n',draq_d.noTrials-full_trials);
         
    case 2
        draq_d.max_laser=max(board_adc_data(4,:));
        draq_d.min_laser=min(board_adc_data(4,:));
        
        %Find the full trials (excluding short trials)
        while at_end==0
            delta_ii_first=find(board_adc_data(4,ii:end)>=0.5*(draq_d.max_laser-draq_d.min_laser)+draq_d.min_laser,1,'first');
            %Found a digital signal
            if ~isempty(delta_ii_first)
                %Find the interval
                ii=ii+delta_ii_first;
                
                if ii+(draq_p.sec_per_trigger-draq_p.sec_before_trigger)*draq_p.ActualRate<length(board_adc_data(4,:))
                    
                    %Full trial
                    draq_d.noTrials=draq_d.noTrials+1;
                    %Trial start time goes in column 1
                    trials_to_sort(draq_d.noTrials,1)=(ii/draq_p.ActualRate)-(draq_p.sec_before_trigger+1);
                    %Block no at start of trial go in column 2
                    trials_to_sort(draq_d.noTrials,2)=ceil((trials_to_sort(draq_d.noTrials,1)*draq_p.ActualRate)/num_samples_per_data_block);
                    %Block no at end of trial go in column 3
                    trials_to_sort(draq_d.noTrials,3)=trials_to_sort(draq_d.noTrials,2)+ceil((draq_p.sec_per_trigger*draq_p.ActualRate)/num_samples_per_data_block);
                    full_trial_start(draq_d.noTrials)=ii-draq_p.ActualRate*(draq_p.sec_before_trigger+1);
                    full_trial_end(draq_d.noTrials)=full_trial_start(draq_d.noTrials)+draq_p.sec_per_trigger*draq_p.ActualRate;
                    ii=ii+ceil((draq_p.sec_per_trigger-draq_p.sec_before_trigger)*draq_p.ActualRate);
                else
                    at_end=1;
                end
            else
                at_end=1;
            end
        end
        
        full_trials=draq_d.noTrials;
         
        fprintf(1, 'Found %d full trials...\n',full_trials);
        
        %Find empty trials
        last_trial=draq_d.noTrials;
        
        trials_to_sort=sortrows(trials_to_sort);
        
        fprintf(1, 'Finding inter trials...\n');
        
        %Empty trial before the first trial
        if trials_to_sort(1,1)>draq_p.sec_per_trigger+draq_p.sec_before_trigger
            draq_d.noTrials=draq_d.noTrials+1;
            trials_to_sort(draq_d.noTrials,1)= trials_to_sort(1,1)-draq_p.sec_per_trigger;
            trials_to_sort(draq_d.noTrials,2)=ceil((trials_to_sort(draq_d.noTrials,1)*draq_p.ActualRate)/num_samples_per_data_block);
            trials_to_sort(draq_d.noTrials,3)=trials_to_sort(draq_d.noTrials,2)+ceil((draq_p.sec_per_trigger*draq_p.ActualRate)/num_samples_per_data_block);
        end
        
        %Empty trials between tirals
        for ii=2:last_trial
            if (trials_to_sort(ii,1)-trials_to_sort(ii-1,1)-draq_p.sec_per_trigger)>draq_p.sec_per_trigger
                
                %Trial after the last one
                if ii~=last_trial
                    draq_d.noTrials=draq_d.noTrials+1;
                    trials_to_sort(draq_d.noTrials,1)= trials_to_sort(ii-1,1)+draq_p.sec_per_trigger;
                    trials_to_sort(draq_d.noTrials,2)=ceil((trials_to_sort(draq_d.noTrials,1)*draq_p.ActualRate)/num_samples_per_data_block);
                    trials_to_sort(draq_d.noTrials,3)=trials_to_sort(draq_d.noTrials,2)+ceil((draq_p.sec_per_trigger*draq_p.ActualRate)/num_samples_per_data_block);
                end
                
                %Trial before the next one
                draq_d.noTrials=draq_d.noTrials+1;
                trials_to_sort(draq_d.noTrials,1)= trials_to_sort(ii,1)-draq_p.sec_per_trigger;
                trials_to_sort(draq_d.noTrials,2)=ceil((trials_to_sort(draq_d.noTrials,1)*draq_p.ActualRate)/num_samples_per_data_block);
                trials_to_sort(draq_d.noTrials,3)=trials_to_sort(draq_d.noTrials,2)+ceil((draq_p.sec_per_trigger*draq_p.ActualRate)/num_samples_per_data_block);
            end
        end
        
        fprintf(1, 'Found %d inter trials ...\n',draq_d.noTrials-full_trials);
         
    case 3
        %dropcnsampler
        %Find the full trials (excluding short trials)
        while at_end==0
            delta_ii_first=find(digital_input_ns(ii:end)==16,1,'first');
            %Found the end of an odor on?
            if ~isempty(delta_ii_first)
                %Find the end of odor on
                ii=ii+delta_ii_first;
                delta_ii_last=find(digital_input_ns(ii:end)~=16,1,'first');
                %Was it found?
                if ~isempty(delta_ii_last)
                    ii=ii+delta_ii_last;
                    %Is it too close to the beggining?
                    if (ii-(draq_p.sec_before_trigger*draq_p.ActualRate))>0
                        %Is is too close to the end?
                        if (ii+((draq_p.sec_per_trigger-draq_p.sec_before_trigger)*draq_p.ActualRate))<length(digital_input_ns)
                            %Full trial
                            draq_d.noTrials=draq_d.noTrials+1;
                            %Trial start time goes in column 1
                            trials_to_sort(draq_d.noTrials,1)=(ii/draq_p.ActualRate)-(draq_p.sec_before_trigger+1);
                            %Block no at start of trial go in column 2
                            trials_to_sort(draq_d.noTrials,2)=ceil((trials_to_sort(draq_d.noTrials,1)*draq_p.ActualRate)/num_samples_per_data_block);
                            %Block no at end of trial go in column 3
                            trials_to_sort(draq_d.noTrials,3)=trials_to_sort(draq_d.noTrials,2)+ceil((draq_p.sec_per_trigger*draq_p.ActualRate)/num_samples_per_data_block);
                            full_trial_start(draq_d.noTrials)=ii-draq_p.ActualRate*(draq_p.sec_before_trigger+1);
                            full_trial_end(draq_d.noTrials)=full_trial_start(draq_d.noTrials)+draq_p.sec_per_trigger*draq_p.ActualRate;
%                             figure(1)
%                             plot(digital_input_ns(full_trial_start(draq_d.noTrials):full_trial_end(draq_d.noTrials)))
                        else
                            at_end=1;
                        end
                    end
                else
                    at_end=1;
                end
            else
                at_end=1;
            end
        end
        
        full_trials=draq_d.noTrials;
        
        fprintf(1, 'Found %d full trials...\n',full_trials);
        
    case 4
        draq_d.max_laser=max(board_adc_data(1,:));
        draq_d.min_laser=min(board_adc_data(1,:));
        
        %Find the full trials (excluding short trials)
        while at_end==0
            delta_ii_first=find(board_adc_data(1,ii:end)>=0.5*(draq_d.max_laser-draq_d.min_laser)+draq_d.min_laser,1,'first');
            %Found a digital signal
            if ~isempty(delta_ii_first)
                %Find the interval
                ii=ii+delta_ii_first;
                
                if ii+(draq_p.sec_per_trigger-draq_p.sec_before_trigger)*draq_p.ActualRate<length(board_adc_data(1,:))
                    
                    %Full trial
                    draq_d.noTrials=draq_d.noTrials+1;
                    %Trial start time goes in column 1
                    trials_to_sort(draq_d.noTrials,1)=(ii/draq_p.ActualRate)-(draq_p.sec_before_trigger+1);
                    %Block no at start of trial go in column 2
                    trials_to_sort(draq_d.noTrials,2)=ceil((trials_to_sort(draq_d.noTrials,1)*draq_p.ActualRate)/num_samples_per_data_block);
                    %Block no at end of trial go in column 3
                    trials_to_sort(draq_d.noTrials,3)=trials_to_sort(draq_d.noTrials,2)+ceil((draq_p.sec_per_trigger*draq_p.ActualRate)/num_samples_per_data_block);
                    full_trial_start(draq_d.noTrials)=ii-draq_p.ActualRate*(draq_p.sec_before_trigger+1);
                    full_trial_end(draq_d.noTrials)=full_trial_start(draq_d.noTrials)+draq_p.sec_per_trigger*draq_p.ActualRate;
                    ii=ii+ceil((draq_p.sec_per_trigger-draq_p.sec_before_trigger)*draq_p.ActualRate);
                else
                    at_end=1;
                end
            else
                at_end=1;
            end
        end
        
        full_trials=draq_d.noTrials;
        
        fprintf(1, 'Found %d full trials...\n',full_trials);
        
        %Find empty trials
        last_trial=draq_d.noTrials;
        
        trials_to_sort=sortrows(trials_to_sort);
        
        fprintf(1, 'Finding inter trials...\n');
        
        %Empty trial before the first trial
        if trials_to_sort(1,1)>draq_p.sec_per_trigger+draq_p.sec_before_trigger
            draq_d.noTrials=draq_d.noTrials+1;
            trials_to_sort(draq_d.noTrials,1)= trials_to_sort(1,1)-draq_p.sec_per_trigger;
            trials_to_sort(draq_d.noTrials,2)=ceil((trials_to_sort(draq_d.noTrials,1)*draq_p.ActualRate)/num_samples_per_data_block);
            trials_to_sort(draq_d.noTrials,3)=trials_to_sort(draq_d.noTrials,2)+ceil((draq_p.sec_per_trigger*draq_p.ActualRate)/num_samples_per_data_block);
        end
        
        %Empty trials between tirals
        for ii=2:last_trial
            if (trials_to_sort(ii,1)-trials_to_sort(ii-1,1)-draq_p.sec_per_trigger)>draq_p.sec_per_trigger
                
                %Trial after the last one
                if ii~=last_trial
                    draq_d.noTrials=draq_d.noTrials+1;
                    trials_to_sort(draq_d.noTrials,1)= trials_to_sort(ii-1,1)+draq_p.sec_per_trigger;
                    trials_to_sort(draq_d.noTrials,2)=ceil((trials_to_sort(draq_d.noTrials,1)*draq_p.ActualRate)/num_samples_per_data_block);
                    trials_to_sort(draq_d.noTrials,3)=trials_to_sort(draq_d.noTrials,2)+ceil((draq_p.sec_per_trigger*draq_p.ActualRate)/num_samples_per_data_block);
                end
                
                %Trial before the next one
                draq_d.noTrials=draq_d.noTrials+1;
                trials_to_sort(draq_d.noTrials,1)= trials_to_sort(ii,1)-draq_p.sec_per_trigger;
                trials_to_sort(draq_d.noTrials,2)=ceil((trials_to_sort(draq_d.noTrials,1)*draq_p.ActualRate)/num_samples_per_data_block);
                trials_to_sort(draq_d.noTrials,3)=trials_to_sort(draq_d.noTrials,2)+ceil((draq_p.sec_per_trigger*draq_p.ActualRate)/num_samples_per_data_block);
            end
        end
        
        fprintf(1, 'Found %d inter trials ...\n',draq_d.noTrials-full_trials);
    
    case 7
        %continuous read
        %Find the full trials (excluding short trials)
        ii=1;
        if ~isempty(digital_input)
            total_length=length(digital_input);
        else
            total_length=board_adc_index;
        end
        while at_end==0
               
            %Is is too close to the end?
%             if (ii+((draq_p.sec_per_trigger-draq_p.sec_before_trigger)*draq_p.ActualRate)+0.1)<total_length
                if (ii+(draq_p.sec_per_trigger*draq_p.ActualRate)+0.1)<total_length
                
                %Note: For filtering I add a 0.2 sec time_pad to avoid filter artifacts.
                %Make sure to use the same time_pad that is used for filtering in drgMaster
                time_pad=0.2;
                
                %Full trial
                draq_d.noTrials=draq_d.noTrials+1;
               
                %Shift ii as if a digital trial had been detected
%                 ii=ii+draq_p.ActualRate*(draq_p.sec_before_trigger+0.5);
                
                %Trial start time goes in column 1
                trials_to_sort(draq_d.noTrials,1)=(ii/draq_p.ActualRate);
                %Block no at start of trial goes in column 2
                trials_to_sort(draq_d.noTrials,2)=ceil((trials_to_sort(draq_d.noTrials,1)*draq_p.ActualRate)/num_samples_per_data_block);
                %Block no at end of trial goes in column 3
                trials_to_sort(draq_d.noTrials,3)=trials_to_sort(draq_d.noTrials,2)+ceil(((draq_p.sec_per_trigger)*draq_p.ActualRate)/num_samples_per_data_block)-1;
                
                full_trial_start(draq_d.noTrials)=ii;
                full_trial_end(draq_d.noTrials)=full_trial_start(draq_d.noTrials)+(draq_p.sec_per_trigger+2*time_pad)*draq_p.ActualRate;
                
                 %Shift ii to the end of the trial minus time_pad
                ii=ii+draq_p.ActualRate*(draq_p.sec_per_trigger);
                pffft=1;
                
                %Old code
%                 %Full trial
%                 draq_d.noTrials=draq_d.noTrials+1;
%                
%                 %Shift ii as if a digital trial had been detected
%                 ii=ii+draq_p.ActualRate*(draq_p.sec_before_trigger+0.5);
%                 
%                 %Trial start time goes in column 1
%                 trials_to_sort(draq_d.noTrials,1)=(ii/draq_p.ActualRate)-(draq_p.sec_before_trigger+0.5);
%                 %Block no at start of trial goes in column 2
%                 trials_to_sort(draq_d.noTrials,2)=ceil((trials_to_sort(draq_d.noTrials,1)*draq_p.ActualRate)/num_samples_per_data_block);
%                 %Block no at end of trial goes in column 3
%                 trials_to_sort(draq_d.noTrials,3)=trials_to_sort(draq_d.noTrials,2)+ceil((draq_p.sec_per_trigger*draq_p.ActualRate)/num_samples_per_data_block);
%                 
%                 full_trial_start(draq_d.noTrials)=ii-draq_p.ActualRate*(draq_p.sec_before_trigger+0.5);
%                 full_trial_end(draq_d.noTrials)=full_trial_start(draq_d.noTrials)+draq_p.sec_per_trigger*draq_p.ActualRate;
%                 
%                  %Shift ii to the end of the trial plus 0.1 sec
%                 ii=ii+draq_p.ActualRate*(draq_p.sec_per_trigger-draq_p.sec_before_trigger+0.1);
                
            else
                at_end=1;
            end
            
            
        end
        
        full_trials=draq_d.noTrials;
        
        fprintf(1, 'Found %d full trials...\n',full_trials);
end
%Sort the trials and assign them to draq_d

if ~isempty(digital_input)
    total_length=length(digital_input);
else
    total_length=board_adc_index;
end
 
fprintf(1, 'Sorting trials...\n');
sorted_trials=sortrows(trials_to_sort);
  
if sorted_trials(1,2)<1
    draq_d.t_trial=sorted_trials(sorted_trials(:,2)>1,1)';
    draq_d.start_blockNo=sorted_trials(sorted_trials(:,2)>1,2);
    draq_d.end_blockNo=sorted_trials(sorted_trials(:,2)>1,3);
    draq_d.noTrials=sum(sorted_trials(:,2)>1);
else
    draq_d.t_trial=sorted_trials(:,1)';
    draq_d.start_blockNo=sorted_trials(:,2);
    draq_d.end_blockNo=sorted_trials(:,3);
end

draq_d.t_end=draq_d.t_trial(end)+draq_p.sec_per_trigger;
draq_d.total_length=total_length;
fprintf(1, 'Done reading rhd header...\n');

return


function a = fread_QString(fid)

% a = read_QString(fid)
%
% Read Qt style QString.  The first 32-bit unsigned number indicates
% the length of the string (in bytes).  If this number equals 0xFFFFFFFF,
% the string is null.

a = '';
length = fread(fid, 1, 'uint32');
if length == hex2num('ffffffff')
    return;
end
% convert length from bytes to 16-bit Unicode words
length = length / 2;

for i=1:length
    a(i) = fread(fid, 1, 'uint16');
end

return


function s = plural(n)

% s = plural(n)
% 
% Utility function to optionally plurailze words based on the value
% of n.

if (n == 1)
    s = '';
else
    s = 's';
end

return


function out = notch_filter(in, fSample, fNotch, Bandwidth)

% out = notch_filter(in, fSample, fNotch, Bandwidth)
%
% Implements a notch filter (e.g., for 50 or 60 Hz) on vector 'in'.
% fSample = sample rate of data (in Hz or Samples/sec)
% fNotch = filter notch frequency (in Hz)
% Bandwidth = notch 3-dB bandwidth (in Hz).  A bandwidth of 10 Hz is
%   recommended for 50 or 60 Hz notch filters; narrower bandwidths lead to
%   poor time-domain properties with an extended ringing response to
%   transient disturbances.
%
% Example:  If neural data was sampled at 30 kSamples/sec
% and you wish to implement a 60 Hz notch filter:
%
% out = notch_filter(in, 30000, 60, 10);

tstep = 1/fSample;
Fc = fNotch*tstep;

L = length(in);

% Calculate IIR filter parameters
d = exp(-2*pi*(Bandwidth/2)*tstep);
b = (1 + d*d)*cos(2*pi*Fc);
a0 = 1;
a1 = -b;
a2 = d*d;
a = (1 + d*d)/2;
b0 = 1;
b1 = -2*cos(2*pi*Fc);
b2 = 1;

out = zeros(size(in));
out(1) = in(1);  
out(2) = in(2);
% (If filtering a continuous data stream, change out(1) and out(2) to the
%  previous final two values of out.)

% Run filter
for i=3:L
    out(i) = (a*b2*in(i-2) + a*b1*in(i-1) + a*b0*in(i) - a2*out(i-2) - a1*out(i-1))/a0;
end

return


function move_to_base_workspace(variable)

% move_to_base_workspace(variable)
%
% Move variable from function workspace to base MATLAB workspace so
% user will have access to it after the program ends.

variable_name = inputname(1);
assignin('base', variable_name, variable);

return;