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