1function [HDR] = save2gdf(arg1,arg2,arg3); 2% SAVE2GDF loads EEG data and saves it in GDF format 3% It has been tested with data of the following formats: 4% Physiobank, BKR, CNT (Neurscan), EDF, 5% 6% HDR = save2gdf(sourcefile [, destfile [, option]]); 7% 8% HDR = save2gdf(HDR,data); 9% 10% 11% see also: SLOAD, SOPEN, SREAD, SCLOSE, SWRITE 12 13% 14% sourcefile sourcefile wildcards are allowed 15% destfile destination file in BKR format 16% if destfile is empty or a directory, sourcefile but with extension .bkr is used. 17% options 18% 'spikes' 19% 'bursts' 20% gain Gain factor for unscaled EEG data (e.g. old Matlab files) 21% 'removeDC' removes mean 22% 'autoscale k:l' uses only channels from k to l for scaling 23% 'detrend k:l' channels from k to l are detrended with an FIR-highpass filter. 24% 'PhysMax=XXX' uses a fixed scaling factor; might be important for concanating BKR files 25% +XXX and -XXX correspond to the maximum and minimum physical value, resp. 26% You can concanate several options by separating with space, komma or semicolon 27% 28% HDR Header, HDR.FileName must contain target filename 29% data data samples 30 31% This program is free software; you can redistribute it and/or 32% modify it under the terms of the GNU General Public License 33% as published by the Free Software Foundation; either version 3 34% of the License, or (at your option) any later version. 35 36 37% $Id$ 38% Copyright (C) 2003-2005,2007,2008 by Alois Schloegl <alois.schloegl@gmail.com> 39% This file is part of the biosig project http://biosig.sf.net/ 40 41 42FLAG_REMOVE_DC = 0; 43FLAG_AUTOSCALE = 0; 44FLAG_DETREND = 0; 45FLAG_PHYSMAX = 0; 46FLAG_removeDrift = 0; 47 48chansel = 0; 49 50if nargin<2, arg2=[]; end; 51if nargin<3, 52 arg3=[]; 53end; 54 55if ischar(arg1), 56 inpath = fileparts(arg1); 57 infile = dir(arg1); % input file 58 if isempty(infile) 59 fprintf(2,'ERROR SAVE2GDF: file %s not found.\n',arg1); 60 return; 61 end; 62 outfile = arg2; 63elseif isstruct(arg1) & isnumeric(arg2), 64 HDR = arg1; 65 data = arg2; 66else %if isstruct(arg1) & isnumeric(arg2), 67 fprintf(2,'Error SAVE2GDF: invalid input arguments\n'); 68 return; 69end; 70 71if isstruct(arg1), 72 %HDR.FileName = destfile; % Assign Filename 73 if isfield(HDR,'NS') 74 if HDR.NS==size(data,2), 75 % It's ok. 76 elseif HDR.NS==size(data,1), 77 warning('data is transposed\n'); 78 data = data'; 79 elseif HDR.NS==size(data,2)+1, 80 HDR.NS = size(data,2); 81 %warning('data is transposed\n'); 82 %data = data'; 83 else 84 fprintf(2,'HDR.NS=%i is not equal to number of data columns %i\n',HDR.NS,size(data,2)); 85 return; 86 end; 87 else 88 HDR.NS = size(data,2); % number of channels 89 end; 90 91 % THRESHOLD, GDFTYP -> Phys/Dig/Min/Max 92 if isa(data,'single') & strncmp(version,'6',1) 93 data = double(data); 94 end; 95 if HDR.NS, 96 if isfield(HDR,'THRESHOLD') 97 HDR.DigMax = HDR.THRESHOLD(1:HDR.NS,2)'; 98 HDR.DigMin = HDR.THRESHOLD(1:HDR.NS,1)'; 99 100 else %if ~isfield(HDR,'THRESHOLD') 101 fprintf(2,'Warning SAVE2GDF: no THRESHOLD value provided - automated overflow detection not supported\n'); 102 103 HDR.DigMax = max(data,[],1); 104 HDR.DigMin = min(data,[],1); 105 end; 106 end; 107 108 HDR.TYPE = 'GDF'; 109 if ~isfield(HDR,'VERSION'); 110 HDR.VERSION = 2.0; 111 end; 112% HDR.FLAG.UCAL = 0; % data is de-calibrated, no rescaling within SWRITE 113 114 if strcmp(HDR.TYPE,'EVENT') 115 HDR.SampleRate = HDR.EVENT.SampleRate; 116 elseif isfield(HDR,'GDFTYP') 117 if any((HDR.GDFTYP>=16) & (HDR.GDFTYP<=18)) 118 if (HDR.VERSION < 1.90) 119 digmin = -(2^30); 120 digmax = 2^30; 121 HDR.PhysMax = [1,HDR.DigMax]*HDR.Calib; 122 HDR.PhysMin = [1,HDR.DigMin]*HDR.Calib; 123 data = (data - repmat(HDR.DigMin(:)',size(data,1),1)); 124 data = data * diag((digmax-digmin)./HDR.Cal) + digmin; 125 HDR.DigMin(:) = digmin; 126 HDR.DigMax(:) = digmax; 127 HDR.Cal = (HDR.PhysMax-HDR.PhysMin)./(HDR.DigMax-HDR.DigMin); 128 HDR.Off = HDR.PhysMin - HDR.Cal .* HDR.DigMin; 129 HDR.Calib = [HDR.Off;diag(HDR.Cal)]; 130 else 131 % do nothing 132 end 133 else 134 HDR.PhysMax = [1,HDR.DigMax]*HDR.Calib; 135 HDR.PhysMin = [1,HDR.DigMin]*HDR.Calib; 136 %bits = ceil(log2(max(HDR.DigMax-HDR.DigMin+1))/8)*8; % allows int8, int16, int24, int32, etc. 137 bits1 = ceil(log2(HDR.DigMax-HDR.DigMin+1)); 138 [datatyp,limits,datatypes] = gdfdatatype(HDR.GDFTYP); 139 bits = log2(limits(:,2)-limits(:,1)+1); 140 fprintf(1,'SAVE2GDF: %i bits needed, %i bits used for file %s\n',max(bits1),max(bits),HDR.FileName); 141 end; 142 143 % re-scale data to account for the scaling factor in the header 144 %HIS = histo3(data); save HIS HIS 145 else 146 tmp = sort(data,1); 147 tmp = diff(tmp); 148 tmp(tmp<8*eps) = NaN; 149 dQ = min(tmp); 150 151 digmax = HDR.DigMax; 152 digmin = HDR.DigMin; 153 154 bits = ceil(log2(max(digmax-digmin+1))); % allows any bit-depth 155 if (min(dQ)<1) & (HDR.VERSION>1.9), GDFTYP = 16; % float32 156 elseif min(dQ)<1, GDFTYP = 5; % int32 157 elseif bits==8, GDFTYP = 1; % int8 158 elseif bits==16, GDFTYP = 3; % int16 159 elseif bits==32, GDFTYP = 5; % int32 160 elseif bits==64, GDFTYP = 7; % int64 161 elseif ~isempty(bits); GDFTYP = 255+bits; % intN 162 else GDFTYP = 3; % int8 163 end; 164 HDR.GDFTYP = GDFTYP; 165 166 if length(HDR.GDFTYP)==HDR.NS, 167 elseif length(HDR.GDFTYP)==1, 168 HDR.GDFTYP = HDR.GDFTYP*ones(1,HDR.NS); % int16 169 else 170 %% PROBLEM 171 end 172 173 [datatyp,limits,datatypes] = gdfdatatype(HDR.GDFTYP); 174 if 1, 175 % here, the data is forced to a different data type 176 % this is useful if data is float 177 % this can cause round-off errors 178 HDR.DigMin = limits(:,1)'; 179 HDR.DigMax = limits(:,2)'; 180 HDR.FLAG.UCAL = 0; % data is calibrated, rescaling within SWRITE 181 182 else 183 % here, the data is of integer type 184 % no round of errors occur. 185 186 c0 = 0; 187 while any(digmin'<limits(:,1)), 188 c = 2^ceil(log2(max(limits(:,1)-digmin'))) 189 digmin = digmin + c; 190 digmax = digmax + c; 191 c0 = c0 + c; 192 end; 193 while any(digmax'>limits(:,2)), 194 c = 2^ceil(log2(max(digmax'-limits(:,2)))) 195 digmin = digmin - c; 196 digmax = digmax - c; 197 c0 = c0 - c; 198 end; 199 while any(digmin'<limits(:,1)), 200 c = 2^ceil(log2(max(limits(:,1)-digmin'))) 201 digmin = digmin + c; 202 digmax = digmax + c; 203 c0 = c0 + c; 204 end; 205 206 data = data + c0; 207 208 HDR.DigMax = digmax; %limits(:,2); %*ones(1,HDR.NS); 209 HDR.DigMin = digmin; %limits(:,1); %*ones(1,HDR.NS); 210 211 %fprintf(1,'Warning SAVE2GDF: overflow detection not implemented, yet.\n'); 212 end; 213 if isfield(HDR,'Calib') & ~isfield(HDR,'PhysMax'); 214 HDR.PhysMax = [1,HDR.DigMax]*HDR.Calib; 215 HDR.PhysMin = [1,HDR.DigMin]*HDR.Calib; 216 end; 217 end; 218 219 if ~isfield(HDR,'Dur'); 220 HDR.Dur = 1/HDR.SampleRate; 221 HDR.SPR = 1; 222 end; 223 224%% [HDR.PhysMax;HDR.PhysMin;HDR.DigMax;HDR.DigMin;max(data);min(data)], 225 226 HDR = sopen(HDR,'w'); 227 if HDR.FILE.FID < 0, 228 fprintf(1,'Error SAVE2GDF: couldnot open file %s.\n',HDR.FileName); 229 return; 230 end; 231 232 if numel(data)>0, 233 HDR = swrite(HDR,data(:,1:HDR.NS)); % WRITE GDF FILE 234 end; 235 HDR = sclose(HDR); 236 237 % final test 238 try 239 240 [y1,H2]=sload(HDR.FileName,0,'UCAL','OFF'); 241 d2 = [ones(size(data,1),1),data]*H2.Calib; 242 if all(all((d2==y1) | (isnan(d2) & isnan(y1)))), 243 fprintf(2,'SAVE2GDF: saving file %s OK.\n',HDR.FileName); 244 else 245 fprintf(2,'SAVE2GDF: file %s saved. Maximum relative roundoff error is %f.\n',HDR.FileName, max(max((d2-y1)./(abs(d2)+abs(y1)))) ); 246 end; 247 catch 248 fprintf(2,'Error SAVE2GDF: saving file %s failed\n',HDR.FileName); 249 end; 250 return; 251end; 252 253if ~isempty(arg3) 254 gain=arg3; 255 [I,J,V]=find(gain); 256 fid = fopen('MM.mmm','w+'); 257 fprintf(fid,'%%%%MatrixMarket matrix coordinate real general\n'); 258 fprintf(fid,'%% generated (C) 2009 by Alois Schloegl\n'); 259 fprintf(fid,'%% selecting channels\n'); 260 fprintf(fid,'%i %i %i\n',size(gain),length(V)); 261 for k=1:length(V), 262 fprintf(fid,'%2i %2i %f\n',I(k),J(k),V(k)); 263 end; 264 fclose(fid); 265end; 266for k=1:length(infile); 267 filename = fullfile(inpath,infile(k).name); 268 [pf,fn,ext] = fileparts(filename); 269 270 % load eeg data 271 %[data,HDR] = sload(filename); 272 HDR = sopen(filename,'r',0); 273 if HDR.FILE.FID<0, 274 fprintf(2,'Error SAVE2GDF: file %s not found\n',filename); 275 return; 276 end; 277 HDR.FLAG.UCAL = 1; 278 HDR.FLAG.OVERFLOWDETECTION = 0; 279 [data,HDR] = sread(HDR,inf); 280 HDR = sclose(HDR); 281 if isfield(HDR,'EDF') 282 if isfield(HDR.EDF,'Annotations') 283 if ~isempty(HDR.EDF.Annotations) 284 fprintf(2,'Warning SAVE2GDF: Annotations in EDF+ are not fully supported.\n'); 285 end; 286 end; 287 end; 288 if ~isfield(HDR,'DigMax'), 289 HDR.DigMax = max(data,[],1); 290 end; 291 if ~isfield(HDR,'DigMin'), 292 HDR.DigMin = min(data,[],1); 293 end; 294 if ~isfield(HDR,'DigMin'), 295 HDR.PhysMax = [1,HDR.DigMax]*HDR.Calib; 296 end; 297 if ~isfield(HDR,'DigMin'), 298 HDR.PhysMin = [1,HDR.DigMin]*HDR.Calib; 299 end; 300 if ~isfield(HDR,'NS'), 301 warning(['number of channels undefined in ',filename]); 302 HDR.NS = size(data,2); 303 end; 304 305 if isempty(outfile), % default destination directory 306 ix = max(find(filename=='.')); 307 HDR.FileName = [HDR.FILE.Name,'.gdf']; % destination directory is current working directory 308 elseif isdir(outfile), % output file 309 HDR.FILE.Path = outfile; 310 HDR.FileName = fullfile(outfile,[HDR.FILE.Name,'.gdf']); 311 else 312 [HDR.FILE.Path,HDR.FILE.Name,Ext] = fileparts(outfile); 313 HDR.FileName = fullfile(HDR.FILE.Path,[HDR.FILE.Name,Ext]); 314 end; 315 HDR=save2gdf(HDR,data); 316end; 317