function [header, groups, data] = c3dread(file_path)

fid = fopen(file_path,'rb');
if fid < 0
    error(['Cannot open file: ' file_path])
end
%Read the header
header_data = int16(fread(fid,256,'int16'));

%Chech if C3D file
if header_data(1) ~= hex2dec('5002') 
    error('The input file is not of C3D format.')
end

header.parameters_start_block = bitand(header_data(1),hex2dec('00FF'));
header.number_trajectories = header_data(2);

if header.number_trajectories == 0
    warning('C3D file does not contain any trajectory, possibly empty file')
end

header.number_analog_measurements_per_frame = double(header_data(3));
header.number_first_frame = double(header_data(4)); %1 based
header.number_last_frame = double(header_data(5)); %1 based
header.max_interpolating_gap = double(header_data(6));
header.data_start_block = double(header_data(9));
header.number_analog_samples_per_frame = double(header_data(10));
header.has_lable_and_range = header_data(148) > 0; %12345 if true, 0 if false
header.label_range_start = double(header_data(149));
header.has_4_char_event_labels = header_data(150) > 0; %12345 if true, 0 if false
header.number_time_events = double(header_data(151));

fseek(fid,12,'bof');
header.scale = fread(fid,1,'float32');

fseek(fid,20,'bof');
header.frame_rate = fread(fid,1,'float32');

header.parameters_start_byte = (header.parameters_start_block-1)*512;
fseek(fid,header.parameters_start_byte,'bof'); %Seek to parameters section

%Check parameters format
paramters_file_format = fread(fid,1,'int16');
if paramters_file_format ~= hex2dec('5001')
    error('The parameters on the file are not of the correct format, possible corrupted file.')
end

header.number_parameters_blocks = fread(fid,1);
header.processor_type = fread(fid,1) - 83;

if header.processor_type ~= 1
    error('Cannot process data floating point format')
end

groups = {};
parameters = {};

while 1
    name_size = fread(fid,1,'int8');
    is_locked = name_size < 0; 
    id = fread(fid,1,'int8');
    name = char(fread(fid,name_size)');
    next_offset = fread(fid,1,'int16');
    next_pos = ftell(fid)+next_offset-2;
    
    if name_size == 0
        break;
    end
    
    if id<0
        %Group
        group.id = abs(id);
        group.name = name;
        
        description_size = fread(fid,1);
        group.description = char(fread(fid,description_size)');
        
        groups{end+1} = group;
    else
        %Parameter
        parameter.group_id = abs(id);
        parameter.name = name;
        parameter.is_locked = is_locked;
        
        %-1 for char, 1 for byte, 2 for int16, 5 for float
        byte_length = fread(fid,1,'int8');
        
        switch byte_length
            case -1
                parameter.data_type = 'char';
            case 1
                parameter.data_type = 'int8';
            case 2
                parameter.data_type = 'int16';
            case 4
                parameter.data_type = 'float32';
            otherwise
                warning('Unrecognized data type found in parameters!')
        end
        
        parameter.n_dimensions = fread(fid,1,'int8');
        parameter.dimensions = fread(fid,parameter.n_dimensions,'int8');
        
        %Sometimes the dimensions of some parameters are 0, and it makes no
        %sense
        %parameter.dimensions(parameter.dimensions==0) = 1; % was trying to
        %corret the issue but this fix does not work, must check the binary
        %file to see if the data is there
        
        parameter.data = fread(fid,prod(parameter.dimensions),parameter.data_type);
        if parameter.n_dimensions > 1 && ~isempty(parameter.dimensions)
            parameter.data = reshape(parameter.data,parameter.dimensions');
        end
        
        if strcmp(parameter.data_type,'char')
            parameter.data = char(parameter.data');
        end
        
        description_size = fread(fid,1);
        parameter.description = char(fread(fid,description_size)');
        
        parameters{end+1} = parameter;
    end
    
    if next_offset == 0
        %Reached the final parameter
        break
    else
        fseek(fid,next_pos,'bof');
    end
end

%Put the parameters into the groups

for g=1:length(groups)
    eval(sprintf('%s = groups{%d};',groups{g}.name,g))
    
    groups{g}.parameters = {};
    
    for p=1:length(parameters)
        
        if parameters{p}.group_id == groups{g}.id
            param = parameters{p};
            param = rmfield(param, 'group_id');
            param = rmfield(param, 'data_type');
            param = rmfield(param, 'n_dimensions');
            param = rmfield(param, 'dimensions');
            
            groups{g}.parameters{end+1} = param;
            eval(sprintf('%s.%s = parameters{%d}.data;',groups{g}.name,parameters{p}.name,p))
        end 
    end
end

header.data_start_byte = (header.data_start_block-1)*512;
fseek(fid,header.data_start_byte,'bof'); %Seek to 3D data section

% if header.number_analog_samples_per_frame > 0
%     error('C3D file contains analog samples! This function is not prepared to read analog sensors')
% end

%Number of labels
n_labels = double(header.number_trajectories);

%Get correct number of frames from POINT.FRAMES or if the varial exists
%TRIAL.ACTUAL_START and ACTUAL_END fields
n_frames = POINT.FRAMES;

if exist('TRIAL','var')
    n_frames = TRIAL.ACTUAL_END_FIELD - TRIAL.ACTUAL_START_FIELD + 1;
    n_frames = n_frames(1);
end

measurements_in_frame = double(n_labels*4+header.number_analog_measurements_per_frame);

% fprintf('%d frames detected.\n',n_frames);
%Read the point data
l_data = fread(fid,measurements_in_frame*n_frames,'float32');
% fprintf('%g MB loaded.\n',measurements_in_frame*n_frames*4/(1024*1024));

%Time variable
time = linspace(0,n_frames/POINT.RATE,n_frames);

%Reshape the data into the 1 frame per line format and add time
data = reshape(l_data,measurements_in_frame,n_frames)';
data = [time' data];
