function [EEG,hdr] = simplebdfreader(file)

IN = fopen(file);

try EEG = eeg_emptyset;
catch, end

% Header
hdr.fileformat = fread(IN, 8, '*char')';
hdr.subjectid = fread(IN, 80, '*char')';
hdr.recordid = fread(IN, 80, '*char')';
hdr.startdate = fread(IN, 8, '*char')';
hdr.starttime = fread(IN, 8, '*char')';
hdr.hdrnbytes = str2num(fread(IN, 8, '*char')');
hdr.dataformat = fread(IN, 44, '*char')';
hdr.nblocks = str2num(fread(IN, 8, '*char')');
hdr.blockdur = str2num(fread(IN, 8, '*char')');
hdr.nchannels = str2num(fread(IN, 4, '*char')');

% Channel headers
hdr.labels = cellstr(fread(IN, [16 hdr.nchannels], '*char')');
[EEG.chanlocs(1:hdr.nchannels).labels] = deal(hdr.labels{:});
hdr.type = cellstr(fread(IN, [80 hdr.nchannels], '*char')');
hdr.unit = cellstr(fread(IN, [8 hdr.nchannels], '*char')');
hdr.physmin = str2num(fread(IN, [8 hdr.nchannels], '*char')');
hdr.physmax = str2num(fread(IN, [8 hdr.nchannels], '*char')');
hdr.digmin = str2num(fread(IN, [8 hdr.nchannels], '*char')');
hdr.digmax = str2num(fread(IN, [8 hdr.nchannels], '*char')');
hdr.gain = (hdr.physmax - hdr.physmin) ./ (hdr.digmax - hdr.digmin);
hdr.filter = cellstr(fread(IN, [80 hdr.nchannels], '*char')');
hdr.nsamples = str2num(fread(IN, [8 hdr.nchannels], '*char')');
hdr.reserved = cellstr(fread(IN, [32 hdr.nchannels], '*char')');

if length(unique(hdr.nsamples)) == 1
    hdr.nsamples = unique(hdr.nsamples);
else
    error('Different number of samples per channel.')
end
EEG.srate = hdr.nsamples / hdr.blockdur;
EEG.xmin = 0;
EEG.xmax = hdr.nblocks * hdr.blockdur * 1000;
EEG.trials = 1;
EEG.ref = 'common';

% Allocate memory (essential!)
EEG.data = zeros(hdr.nchannels, hdr.nsamples * hdr.nblocks, 'single');

% h = waitbar(0, 'Reading bdf file: 0% done');
tic

for block = 0:hdr.nblocks - 1
    buf = fread(IN, [3, hdr.nsamples * hdr.nchannels], '*uint8');
    EEG.data(1:hdr.nchannels, hdr.nsamples * block + 1:hdr.nsamples * block + hdr.nsamples) = ...
        reshape(single(buf(1, :)) + 256 * single(buf(2, :)) + 65536 * single(buf(3, :)) - 131072 * single(bitand(buf(3, :), 128)), ...
        hdr.nsamples, hdr.nchannels)';
%     if (block + 1) == ceil(floor((block + 1) / (hdr.nblocks / 100)) * (hdr.nblocks / 100))
%         p = (block + 1) / hdr.nblocks;
%         t = toc;
%         waitbar(p, h, ['Reading bdf file: ' num2str(round(p * 100)) '% done, ' num2str(round((1 - p) / p * t)) ' s left']);
%     end
end

toc
% close(h);

for chan = 1:hdr.nchannels
    EEG.data(chan, :) = hdr.gain(chan) * EEG.data(chan, :);
end

toc

fclose(IN);
