Initial commit
Esse commit está contido em:
Diego Alvarez-Estevez
@@ -0,0 +1,30 @@
|
||||
% Created by Diego Alvarez-Estevez (http://dalvarezestevez.com)
|
||||
% Last modified 19/01/2018
|
||||
|
||||
% LogFloat conversion (see www.edfplus.info/specs/edffloat.html for
|
||||
% more information)
|
||||
|
||||
%% This program is free software: you can redistribute it and/or modify
|
||||
%% it under the terms of the GNU General Public License as published by
|
||||
%% the Free Software Foundation, either version 3 of the License, or
|
||||
%% (at your option) any later version.
|
||||
|
||||
%% This program is distributed in the hope that it will be useful,
|
||||
%% but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
%% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
%% GNU General Public License for more details.
|
||||
|
||||
%% You should have received a copy of the GNU General Public License
|
||||
%% along with this program. If not, see <http://www.gnu.org/licenses/>.
|
||||
|
||||
function out = LogFloatVector(vector, Y0, A)
|
||||
|
||||
out = zeros(size(vector));
|
||||
|
||||
r = (log10(vector(vector > Y0)) - log10(Y0)) / A;
|
||||
out(vector > Y0) = int16(round(min(r, intmax('int16'))));
|
||||
|
||||
r = (-log10((-1)*vector(vector < Y0)) + log10(Y0)) / A;
|
||||
out(vector < Y0) = int16(round(max(r, -intmax('int16'))));
|
||||
|
||||
% Those with zero value at input remain zero at output
|
||||
@@ -0,0 +1,127 @@
|
||||
% Created by Diego Alvarez-Estevez (http://dalvarezestevez.com)
|
||||
% Last modified 19/01/2018
|
||||
|
||||
% Extracts EEG data from an OpenBCI exported file and saves it to EDF format (check www.edfplus.info for more information on EDF/EDF+ format)
|
||||
|
||||
%% This program is free software: you can redistribute it and/or modify
|
||||
%% it under the terms of the GNU General Public License as published by
|
||||
%% the Free Software Foundation, either version 3 of the License, or
|
||||
%% (at your option) any later version.
|
||||
|
||||
%% This program is distributed in the hope that it will be useful,
|
||||
%% but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
%% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
%% GNU General Public License for more details.
|
||||
|
||||
%% You should have received a copy of the GNU General Public License
|
||||
%% along with this program. If not, see <http://www.gnu.org/licenses/>.
|
||||
|
||||
%%----------------------------------------------------------------------
|
||||
|
||||
% Clean possible previously-existing results
|
||||
clear all;
|
||||
clc;
|
||||
|
||||
%% USER CONFIGURABLE PARAMETERS
|
||||
|
||||
% Skip first (possible invalid, such as comments, etc.) rows of the OpenBCI
|
||||
% export file. Note this is file-dependent, use a "more" (windows) or "cat" (unix/linux) alike commands
|
||||
% to check the real beginning of the recording
|
||||
skipFirstLines = 0;
|
||||
% Sampling frequency of the signals, depends on OpenBCI configuration, default is 250 Hz
|
||||
fs = 250;
|
||||
% Number of channels (default in OpenBCI is 8)
|
||||
nchannels = 8;
|
||||
% Signal filtering paramaters
|
||||
highpass_f = 0.1; % High-pass cut-off frequency
|
||||
notch_fc = 50; % Frequency of mains interference (50 Hz for Europe)
|
||||
notch_w = 1; % Notch filter bandwith
|
||||
|
||||
%% CONVERSION SCRIPT STARTS
|
||||
|
||||
% Ask the user to select the source file
|
||||
[filename, path] = uigetfile(['.', filesep, '*.*'], 'Select file in OpenBCI format');
|
||||
if (filename == 0)
|
||||
return;
|
||||
end
|
||||
fullfilename = fullfile(path, filename);
|
||||
|
||||
% Asks some necessary data to the user
|
||||
disp('Please provide the following information to complete the header information of the EDF file');
|
||||
disp('Check www.edfplus.info/specs/edf.html for detailed information');
|
||||
% Examples
|
||||
% patId = 'XXX';
|
||||
% recId = 'XXX';
|
||||
% startdate = '23.06.17';
|
||||
% starttime = '00.08.52';
|
||||
patId = input('\nSpecify the local patient identification string\n(max 80 characters, press [ENTER] to continue): ', 's');
|
||||
recId = input('\nSpecify the local recording identification string\n(max 80 characters, press [ENTER] to continue): ' , 's');
|
||||
startdate = input('\nSpecify the Startdate of recording\n(DD.MM.YY, press [ENTER] to continue): ', 's');
|
||||
starttime = input('\nSpecify the Starttime of recording\n(HH.MM.SS, press [ENTER] to continue): ', 's');
|
||||
|
||||
% Read OpenBCI data
|
||||
fprintf(1, '\nReading OpenBCI data...');
|
||||
data = csvread(fullfilename, skipFirstLines, 0);
|
||||
eegdata = data(:, 2:9);
|
||||
|
||||
% Apply signal filtering
|
||||
fprintf(1, '\nFiltering signals...\n');
|
||||
for k = 1:nchannels
|
||||
f_eeg(:, k) = polymanHighPassFilter(eegdata(:, k), fs, 1, highpass_f);
|
||||
|
||||
f_eeg(:, k) = polymanNotchFilter(f_eeg(:, k), fs, notch_fc, notch_w);
|
||||
end
|
||||
|
||||
% Export unfiltered data
|
||||
disp('Writing signals (unfiltered) to EDF file...');
|
||||
filename1 = [filename, '-raw.edf'];
|
||||
signals = [];
|
||||
snames = [];
|
||||
units = [];
|
||||
prefilterings = [];
|
||||
for k = 1:nchannels
|
||||
signals{k} = eegdata(:, k);
|
||||
snames{k} = sprintf('EEG CH%d', k);
|
||||
fsamps(k) = fs;
|
||||
Pmins(k) = min(eegdata(:, k));
|
||||
Pmaxs(k) = max(eegdata(:, k));
|
||||
Dmins(k) = -32768;
|
||||
Dmaxs(k) = 32767;
|
||||
units{k} = 'uV';
|
||||
prefilterings{k} = '';
|
||||
end
|
||||
|
||||
statusOK = signals2EDF(signals, fsamps, patId, recId, snames, startdate, starttime, Pmins, Pmaxs, Dmins, Dmaxs, units, prefilterings, filename1);
|
||||
|
||||
if not(statusOK)
|
||||
error('Problem exporting unfiltered signals to EDF');
|
||||
else
|
||||
disp('...done!');
|
||||
end
|
||||
|
||||
% Export filtered data
|
||||
disp('Writing signals (filtered) to EDF file...');
|
||||
filename1 = [filename, '-filt.edf'];
|
||||
signals = [];
|
||||
snames = [];
|
||||
units = [];
|
||||
prefilterings = [];
|
||||
for k = 1:nchannels
|
||||
signals{k} = f_eeg(:, k);
|
||||
snames{k} = sprintf('EEG CH%d', k);
|
||||
fsamps(k) = fs;
|
||||
Pmins(k) = min(f_eeg(:, k));
|
||||
Pmaxs(k) = max(f_eeg(:, k));
|
||||
Dmins(k) = -32768;
|
||||
Dmaxs(k) = 32767;
|
||||
units{k} = 'uV';
|
||||
prefilterings{k} = sprintf('HP(%g) Notch(%g/%g)', highpass_f, notch_fc, notch_w);
|
||||
end
|
||||
|
||||
statusOK = signals2EDF(signals, fsamps, patId, recId, snames, startdate, starttime, Pmins, Pmaxs, Dmins, Dmaxs, units, prefilterings, filename1);
|
||||
|
||||
if not(statusOK)
|
||||
error('Problem exporting filtered signals to EDF');
|
||||
else
|
||||
disp('...done!')
|
||||
end
|
||||
@@ -0,0 +1,38 @@
|
||||
% Created by Diego Alvarez-Estevez (http://dalvarezestevez.com)
|
||||
% Last modified 19/01/2018
|
||||
|
||||
% Implements high pass filter from Polyman
|
||||
|
||||
%% This program is free software: you can redistribute it and/or modify
|
||||
%% it under the terms of the GNU General Public License as published by
|
||||
%% the Free Software Foundation, either version 3 of the License, or
|
||||
%% (at your option) any later version.
|
||||
|
||||
%% This program is distributed in the hope that it will be useful,
|
||||
%% but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
%% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
%% GNU General Public License for more details.
|
||||
|
||||
%% You should have received a copy of the GNU General Public License
|
||||
%% along with this program. If not, see <http://www.gnu.org/licenses/>.
|
||||
|
||||
function signalFilt = polymanHighPassFilter(signal, sr, gain, hpFreq)
|
||||
|
||||
% Check for valid settings
|
||||
if ((hpFreq <= 0) || (hpFreq > sr/2))
|
||||
disp('Error');
|
||||
return;
|
||||
end
|
||||
|
||||
r = 2 * pi * hpFreq * ((1/sr)/2);
|
||||
t = gain / (r + 1);
|
||||
|
||||
b = [t, - t]; % b coeffs according to filter function description (for x(n))
|
||||
a = [1, -(r-1)/(r+1)]; % a coeffs according to filter function description (for y(n))
|
||||
|
||||
% We multiply a(2:end)*(-1) because filter function is implemented as a direct
|
||||
% form II transposed structure
|
||||
a(2:end) = -a(2:end);
|
||||
|
||||
signalFilt = filter(b, a, signal);
|
||||
|
||||
@@ -0,0 +1,43 @@
|
||||
% Created by Diego Alvarez-Estevez (http://dalvarezestevez.com)
|
||||
% Last modified 19/01/2018
|
||||
|
||||
% Implements Notch filter from Polyman
|
||||
|
||||
%% This program is free software: you can redistribute it and/or modify
|
||||
%% it under the terms of the GNU General Public License as published by
|
||||
%% the Free Software Foundation, either version 3 of the License, or
|
||||
%% (at your option) any later version.
|
||||
|
||||
%% This program is distributed in the hope that it will be useful,
|
||||
%% but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
%% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
%% GNU General Public License for more details.
|
||||
|
||||
%% You should have received a copy of the GNU General Public License
|
||||
%% along with this program. If not, see <http://www.gnu.org/licenses/>.
|
||||
|
||||
function signalFilt = polymanNotchFilter(signal, sr, centerFreq, bandWidth)
|
||||
|
||||
% Check for valid settings
|
||||
if ((centerFreq <= 0) || (centerFreq > sr/2))
|
||||
error('Center frequency is not in the Nyquist interval');
|
||||
elseif ((bandWidth <= 0) || (bandWidth > sr/4))
|
||||
error('Not valid bandwidth');
|
||||
end
|
||||
|
||||
r = 2 * pi * centerFreq/sr;
|
||||
s = 2 * pi * (centerFreq - (bandWidth/2))/sr;
|
||||
t = 1 - sqrt((cos(r)-cos(s))^2 + (sin(r)-sin(s))^2) * sqrt(exp(1)^2-1);
|
||||
x = cos(r);
|
||||
y = sin(r);
|
||||
|
||||
b = [1, -2*x, x^2+y^2]; % b coeffs according to filter function description (for x(n))
|
||||
a = [1, 2*t*x, -((t*x)^2 + (t*y)^2)]; % a coeffs according to filter function description (for y(n))
|
||||
|
||||
% We multiply a(2:end)*(-1) because filter function is implemented as a direct
|
||||
% form II transposed structure
|
||||
a(2:end) = -a(2:end);
|
||||
|
||||
signalFilt = filter(b, a, signal);
|
||||
|
||||
|
||||
+178
@@ -0,0 +1,178 @@
|
||||
% Created by Diego Alvarez-Estevez (http://dalvarezestevez.com)
|
||||
% Last modified 19/01/2018
|
||||
|
||||
% Saves signals contained on a Matlab structure to EDF format
|
||||
|
||||
%% This program is free software: you can redistribute it and/or modify
|
||||
%% it under the terms of the GNU General Public License as published by
|
||||
%% the Free Software Foundation, either version 3 of the License, or
|
||||
%% (at your option) any later version.
|
||||
|
||||
%% This program is distributed in the hope that it will be useful,
|
||||
%% but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
%% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
%% GNU General Public License for more details.
|
||||
|
||||
%% You should have received a copy of the GNU General Public License
|
||||
%% along with this program. If not, see <http://www.gnu.org/licenses/>.
|
||||
|
||||
function statusok = signals2EDF(signals, srs, patId, recId, snames, startdate, starttime, Pmins, Pmaxs, Dmins, Dmaxs, units, prefilterings, filename)
|
||||
|
||||
fid = fopen(filename, 'w', 'ieee-le');
|
||||
|
||||
if (fid == -1)
|
||||
fprintf(1, 'Error creating output file %s\n', filename);
|
||||
return;
|
||||
end
|
||||
|
||||
numSignals = length(srs);
|
||||
|
||||
% Calculate recording length (in seconds)
|
||||
recLength = -1;
|
||||
for k = 1:numSignals
|
||||
recLength = max(recLength, length(signals{k})/srs(k));
|
||||
end
|
||||
recLength = ceil(recLength); % Ceil to an integer number
|
||||
|
||||
% Asuming 1s databloack duration. Check possible incompatibilities
|
||||
blockSizeBytes = sum(2*srs(:));
|
||||
|
||||
if (blockSizeBytes > 61440)
|
||||
error('Yet to be implemented: Signals cannot fit on a 1s datablock. Check for other block size possibilities');
|
||||
else
|
||||
blockSize = 1;
|
||||
numBlocks = recLength;
|
||||
end
|
||||
|
||||
general_header_size = 256; %bytes
|
||||
one_signal_header_size = 256; %bytes
|
||||
|
||||
% Write edf
|
||||
|
||||
% FIXED HEADER
|
||||
header.version = 0;
|
||||
header.local_patient_identification = patId;
|
||||
header.local_recording_identification = recId;
|
||||
header.startdate_recording = startdate;
|
||||
header.starttime_recording = starttime;
|
||||
header.num_signals = numSignals;
|
||||
header.num_bytes_header = general_header_size + one_signal_header_size * numSignals;
|
||||
header.reserved = '';
|
||||
header.duration_data_record = blockSize;
|
||||
header.num_data_records = numBlocks;
|
||||
|
||||
fprintf(fid, trimAndFillWithBlanks(num2str(header.version), 8)); % version
|
||||
fprintf(fid, '%-80s', header.local_patient_identification);
|
||||
fprintf(fid, '%-80s', header.local_recording_identification);
|
||||
fprintf(fid, '%-8s', header.startdate_recording);
|
||||
fprintf(fid, '%-8s', header.starttime_recording);
|
||||
fprintf(fid, trimAndFillWithBlanks(num2str(header.num_bytes_header), 8));
|
||||
fprintf(fid, '%-44s', header.reserved);
|
||||
fprintf(fid, trimAndFillWithBlanks(num2str(header.num_data_records), 8));
|
||||
fprintf(fid, trimAndFillWithBlanks(num2str(header.duration_data_record), 8));
|
||||
fprintf(fid, trimAndFillWithBlanks(num2str(header.num_signals), 4));
|
||||
|
||||
% SIGNAL DEPENDENT HEADER
|
||||
signalOffsets = zeros(1, numSignals); % In bytes
|
||||
for k = 1:numSignals
|
||||
header.signals_info(k).label = snames{k};
|
||||
header.signals_info(k).transducer_type = '';
|
||||
header.signals_info(k).physical_dimension = units{k};
|
||||
header.signals_info(k).physical_min = Pmins(k);
|
||||
header.signals_info(k).physical_max = Pmaxs(k);
|
||||
header.signals_info(k).digital_min = Dmins(k);
|
||||
header.signals_info(k).digital_max = Dmaxs(k);
|
||||
header.signals_info(k).prefiltering = prefilterings{k};
|
||||
header.signals_info(k).num_samples_datarecord = srs(k)*blockSize;
|
||||
header.signals_info(k).reserved = '';
|
||||
% NOTE: The two following are not specific EDF header fields, but are practical for EDF handling
|
||||
header.signals_info(k).sample_rate = header.signals_info(k).num_samples_datarecord / header.duration_data_record;
|
||||
if (k > 1)
|
||||
signalOffsets(k) = signalOffsets(k - 1) + 2 * header.signals_info(k - 1).num_samples_datarecord;
|
||||
end
|
||||
header.signals_info(k).signalOffset = signalOffsets(k);
|
||||
end
|
||||
|
||||
% Write signal-dependent header to file
|
||||
for k = 1:numSignals
|
||||
fprintf(fid, '%-16s', header.signals_info(k).label);
|
||||
end
|
||||
for k = 1:numSignals
|
||||
fprintf(fid, '%-80s', header.signals_info(k).transducer_type);
|
||||
end
|
||||
for k = 1:numSignals
|
||||
fprintf(fid, '%-8s', header.signals_info(k).physical_dimension);
|
||||
end
|
||||
for k = 1:numSignals
|
||||
fprintf(fid, trimAndFillWithBlanks(num2str(header.signals_info(k).physical_min), 8));
|
||||
end
|
||||
for k = 1:numSignals
|
||||
fprintf(fid, trimAndFillWithBlanks(num2str(header.signals_info(k).physical_max), 8));
|
||||
end
|
||||
for k = 1:numSignals
|
||||
fprintf(fid, trimAndFillWithBlanks(num2str(header.signals_info(k).digital_min), 8));
|
||||
end
|
||||
for k = 1:numSignals
|
||||
fprintf(fid, trimAndFillWithBlanks(num2str(header.signals_info(k).digital_max), 8));
|
||||
end
|
||||
for k = 1:numSignals
|
||||
fprintf(fid, '%-80s', header.signals_info(k).prefiltering);
|
||||
end
|
||||
for k = 1:numSignals
|
||||
fprintf(fid, trimAndFillWithBlanks(num2str(header.signals_info(k).num_samples_datarecord), 8));
|
||||
end
|
||||
for k = 1:numSignals
|
||||
fprintf(fid, '%-32s', header.signals_info(k).reserved);
|
||||
end
|
||||
|
||||
% Check data starting point
|
||||
current_position = ftell(fid); % in bytes
|
||||
if ne(header.num_bytes_header, current_position)
|
||||
disp('Something wrong could be happening: unexpected position at the beginning of the first data block');
|
||||
end
|
||||
|
||||
bytes_full_data_record = 2 * sum([header.signals_info.num_samples_datarecord]);
|
||||
|
||||
% DATA WRITING
|
||||
for k = 1:numBlocks
|
||||
|
||||
% Initialize datablock
|
||||
data = zeros(1, bytes_full_data_record/2); % Num samples per data record
|
||||
|
||||
for k1 = 1:numSignals
|
||||
|
||||
offsetSignal = (k - 1) * header.signals_info(k1).num_samples_datarecord + 1;
|
||||
onsetSignal = min(offsetSignal + header.signals_info(k1).num_samples_datarecord - 1, length(signals{k1}));
|
||||
|
||||
offsetDataBlock = header.signals_info(k1).signalOffset/2 + 1;
|
||||
onsetDataBlock = offsetDataBlock + length(offsetSignal:onsetSignal) - 1;
|
||||
|
||||
logConversion = strcmp(header.signals_info(k1).physical_dimension', 'Filtered');
|
||||
if logConversion
|
||||
% Parse prefiltering to set this values
|
||||
tokens = textscan(header.signals_info(k1).prefiltering, 'sign*LN[sign*(at %.1fHz)/(%.5f)]/(%.5f)(Kemp:J Sleep Res 1998-supp2:132)');
|
||||
if isempty(tokens{2})
|
||||
disp('Warning: assigned default values to logConversion');
|
||||
LogFloatY0 = 0.0001; % default value
|
||||
else
|
||||
LogFloatY0 = tokens{2};
|
||||
end
|
||||
if isempty(tokens{3})
|
||||
disp('Warning: assigned default values to logConversion');
|
||||
LogFloatA = 0.001; % default value
|
||||
else
|
||||
LogFloatA = tokens{3};
|
||||
end
|
||||
|
||||
% Actual writing
|
||||
data(offsetDataBlock:onsetDataBlock) = Dmins(k1) + (Dmaxs(k1) - Dmins(k1)) * (LogFloatVector(signals{k1}(offsetSignal:onsetSignal), LogFloatY0, LogFloatA) - Pmins(k1))/(Pmaxs(k1) - Pmins(k1));
|
||||
else
|
||||
data(offsetDataBlock:onsetDataBlock) = Dmins(k1) + (Dmaxs(k1) - Dmins(k1)) * (signals{k1}(offsetSignal:onsetSignal) - Pmins(k1))/(Pmaxs(k1) - Pmins(k1));
|
||||
end
|
||||
end
|
||||
data = typecast(int16(data), 'uint16'); % From double to 16bit unsigned integers
|
||||
|
||||
fwrite(fid, data, 'uint16');
|
||||
end
|
||||
|
||||
statusok = (fclose(fid) == 0);
|
||||
@@ -0,0 +1,34 @@
|
||||
% Created by Diego Alvarez-Estevez (http://dalvarezestevez.com)
|
||||
% Last modified 19/01/2018
|
||||
|
||||
% Trims, justifies, and fills with blanks a character vector
|
||||
|
||||
%% This program is free software: you can redistribute it and/or modify
|
||||
%% it under the terms of the GNU General Public License as published by
|
||||
%% the Free Software Foundation, either version 3 of the License, or
|
||||
%% (at your option) any later version.
|
||||
|
||||
%% This program is distributed in the hope that it will be useful,
|
||||
%% but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||||
%% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||||
%% GNU General Public License for more details.
|
||||
|
||||
%% You should have received a copy of the GNU General Public License
|
||||
%% along with this program. If not, see <http://www.gnu.org/licenses/>.
|
||||
|
||||
function result = trimAndFillWithBlanks(str, maxLength, justify)
|
||||
|
||||
if (nargin == 2)
|
||||
justify = 'left';
|
||||
end
|
||||
|
||||
if length(str) > maxLength
|
||||
result = str(1:maxLength);
|
||||
else
|
||||
if strcmp(justify, 'right')
|
||||
result = [blanks(maxLength - length(str)), str];
|
||||
else
|
||||
% default is left
|
||||
result = [str, blanks(maxLength - length(str))];
|
||||
end
|
||||
end
|
||||
Referência em uma Nova Issue
Bloquear um usuário