-
Notifications
You must be signed in to change notification settings - Fork 2
/
read_ctf_res4.m
251 lines (226 loc) · 8.14 KB
/
read_ctf_res4.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
function hdr = read_ctf_res4(fname)
% READ_CTF_RES4 reads the header in RES4 format from a CTF dataset
%
% Use as
% [hdr] = read_ctf_res4(filename)
%
% See also READ_CTF_MEG4
% "VSM MedTech Ltd. authorizes the release into public domain under the
% GPL licence of the Matlab source code files "read_ctf_res4.m" and
% "read_ctf_meg4.m" by the authors of said files from the F.C. Donders
% Centre, Nijmegen, The Netherlands."
% Author(s): Jim McKay November 1999
% Last revision: Jim McKay
% Copyright (c) 1999-2000 CTF Systems Inc. All Rights Reserved.
%
% modifications Copyright (C) 2002, Ole Jensen
% modifications Copyright (C) 2003, Robert Oostenveld
%
% $Log: read_ctf_res4.m,v $
% Revision 1.12 2005/07/28 15:12:27 roboos
% fixed bug in hdr.timeVec (thanks to Sanja Kovacevic)
%
% Revision 1.11 2005/05/26 09:58:08 roboos
% removed the construction of grad, that is now done in a separate fieldtrip function (ctf2grad)
%
% Revision 1.10 2005/05/24 07:34:37 roboos
% added sensType to the output header
%
% Revision 1.9 2005/05/23 11:20:17 roboos
% fixed bug for run description which is longer than 256 characters, improved reading of filter information (thanks to Durk Talsma)
%
% Revision 1.8 2005/04/27 06:18:21 roboos
% added support for MEG42RS format, which seems to work (only tested on a single dataset from MIND Institute in Albuquerque, NM, USA)
%
% Revision 1.7 2005/02/18 13:16:58 roboos
% VSM MedTech Ltd. authorised the release of this code in the public domain
% updated the copyrights, updated the help
%
% Revision 1.6 2004/07/02 11:43:32 roboos
% typographic change in comment
%
% Revision 1.5 2004/06/28 07:34:41 roberto
% added some extra low-level output fields to header, added cm units
%
% Revision 1.4 2003/04/01 07:51:53 roberto
% fixed bug in gradiometer channel labels
%
% Revision 1.3 2003/03/28 17:27:27 roberto
% renamed output gradiometer from gradHC to grad
%
% Revision 1.2 2003/03/24 12:34:33 roberto
% minor changes
%
% Revision 1.1 2003/03/13 14:27:09 roberto
% separated read_ctf_ds into res4-header part and data part
% restructured the res4 part, using (incorrect) CTF documentation
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% read header information
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fid = fopen(fname,'r','ieee-be');
% Check if header file exist
if fid == -1
errMsg = strcat('Could not open header file:',fname);
error(errMsg);
end
% First 8 bytes contain filetype, check is fileformat is correct.
% This function was written for MEG41RS, but also seems to work for MEG42RS.
r_head=setstr(fread(fid,8,'char'))';
if ~strcmp(r_head(1,1:7),'MEG41RS') & ~strcmp(r_head(1,1:7),'MEG42RS')
fclose(fid)
errMsg = sprintf('Resource file format (%s) is not supported for file %s', r_head(1,1:7), fname);
error(errMsg);
end %if
% Read the initial parameters
appName = setstr(fread(fid,256,'char'))' ;
dataOrigin = setstr(fread(fid,256,'char'))' ;
dataDescrip = setstr(fread(fid,256,'char'))' ;
no_trial_avgd = fread(fid,1,'int16') ;
data_time = setstr(fread(fid,255,'char'))';
data_date = setstr(fread(fid,255,'char'))';
fseek(fid,1288,'bof');
% Read the general recording parameters
no_samples = fread(fid,1,'int32');
no_channels = fread(fid,1,'int16');
fseek(fid,2,'cof'); % hole of 2 bytes due to improper alignment
sample_rate = fread(fid,1,'double');
epoch = fread(fid,1,'double');
no_trials = fread(fid,1,'int16');
fseek(fid,2,'cof'); % hole of 2 bytes due to improper alignment
preTrigpts=fread(fid,1,'int32');
fseek(fid,1360,'bof');
% read in the meg4Filesetup structure
run_name = setstr(fread(fid,32,'char')');
run_title = setstr(fread(fid,256,'char')');
instruments = setstr(fread(fid,32,'char')');
coll_desc = setstr(fread(fid,32,'char')');
subj_id = setstr(fread(fid,32,'char')');
operator = setstr(fread(fid,32,'char')') ;
sensFilename = setstr(fread(fid,60,'char')') ;
% not nececssary to seek, the file pointer is already at the desired location
% fseek(fid,1836,'bof');
% Read in the run description length
rd_len=fread(fid,1,'uint32');
% Go to the run description and read it in
fseek(fid,1844,'bof');
run_desc=setstr(fread(fid,rd_len,'char')');
% read in the filter information
num_filt=fread(fid,1,'uint16');
for fi=0:(num_filt-1),
%filt_info=fread(fid,18,'uint8');
filt_freq =fread(fid,1, 'double');
filt_class=fread(fid,1, 'uint32');
filt_type =fread(fid,1, 'uint32');
num_fparm=fread(fid, 1, 'uint16');
%num_fparm=filt_info(18);
if num_fparm ~= 0,
filt_parm=fread(fid,8*num_fparm,'uint8');
end % if
end % for fi
% Read in the channel names
for i=1:no_channels,
temp=fread(fid,32,'char')';
temp(find(temp<32 )) = ' '; % remove non-printable characters
temp(find(temp>126)) = ' '; % remove non-printable characters
endstr = findstr(temp, '-'); temp(endstr:end) = ' '; % cut off at '-'
endstr = findstr(temp, ' '); temp(endstr:end) = ' '; % cut off at ' '
chan_name(i,:) = char(temp); % as char array
chan_label{i} = deblank(char(temp)); % as cell array
end %for
% pre-allocate some memory space
sensGain = zeros([no_channels,1]);
qGain = zeros([no_channels,1]);
ioGain = zeros([no_channels,1]);
sensType = zeros([no_channels,1]);
% Read in the sensor information
fp = ftell(fid);
for chan=1:no_channels,
fread(fid,1,'uint8'); % Read and ignore 1 byte from enum
sensType(chan)=fread(fid,1,'uint8'); % Read sensor type
fread(fid,2,'uint8'); % Read and ignore originalRunNum
fread(fid,4,'uint8'); % Read and ignore coilShape
sensGain(chan)=fread(fid,1,'double'); % Read sensor gain in Phi0/Tesla
qGain(chan)=fread(fid,1,'double'); % Read qxx gain (usually 2^20 for Q20)
ioGain(chan)=fread(fid,1,'double'); % Read i/o gain of special sensors (usually 1.0)
ioOffset(chan)=fread(fid,1,'double');
numCoils(chan)=fread(fid,1,'int16');
grad_order_no(chan)=fread(fid,1,'int16');
fread(fid,4,'char');
% read the coil positions and orientations
for i=1:8
Chan(chan).coil(i).pos = fread(fid,3,'double')';
fread(fid,1,'double');
Chan(chan).coil(i).ori = fread(fid,3,'double')';
fread(fid,3,'double');
end
% read the coil positions and orientations in head coordinates(?)
for i=1:8
Chan(chan).coilHC(i).pos = fread(fid,3,'double')';
fread(fid,1,'double');
Chan(chan).coilHC(i).ori = fread(fid,3,'double')';
fread(fid,3,'double');
end
% jump to the next sensor info record
fseek(fid, fp+chan*1328, 'bof');
end % for chan
% close the header file
fclose(fid);
% determine the different channel types
rowMEG = [];
rowEEG = [];
rowTRIG = [];
rowREF = [];
rowALL = 1:no_channels;
for k=1:no_channels
if strcmp(char(chan_name(k,1)),'M')
rowMEG = [rowMEG k];
end
if strcmp(char(chan_name(k,1:3)),'EEG')
rowEEG = [rowEEG k];
end
if strcmp(char(chan_name(k,1:4)),'STIM')
rowTRIG = [rowTRIG k];
end
if chan_name(k,1)=='B' | chan_name(k,1)=='G' | chan_name(k,1)=='P' | chan_name(k,1)=='Q' | chan_name(k,1)=='R'
rowREF = [rowREF k];
end
end
% FIXME: instead of looking at the name, this should be done using sensType
% according to Tom Holroyd, this should be
%
% UPPT001 is 11
% UTRG001 is 11
% SCLK01 is 17
% meg channels are 5, refmag 0, refgrad 1, adcs 18.
% STIM is 11
% SCLK01 is 17
% EEG057 is 9
% ADC06 is 18
% ADC07 is 18
% ADC16 is 18
% V0 is 15
% assign all the variables that should be outputted as header information
hdr.Fs = sample_rate;
hdr.nChans = no_channels;
hdr.nSamples = no_samples;
hdr.nSamplesPre = preTrigpts;
hdr.timeVec = (1:no_samples)/sample_rate - preTrigpts/sample_rate - 1/sample_rate;
hdr.nTrials = no_trials;
hdr.gainV = ioGain./(qGain.*sensGain);
hdr.ioGain = ioGain;
hdr.qGain = qGain;
hdr.sensGain = sensGain;
hdr.sensType = sensType;
hdr.label = chan_label(:);
hdr.nameALL = chan_name;
hdr.rowMEG = rowMEG;
hdr.rowEEG = rowEEG;
hdr.rowTRIG = rowTRIG;
hdr.rowREF = rowREF;
hdr.Chan = Chan;
% hdr.nameEEG = [];
% hdr.nameMEG = [];
% hdr.nameEOG = [];
% hdr.trigV = [];
% hdr.SwapData = [];