-
Notifications
You must be signed in to change notification settings - Fork 16
/
PhilipsRead_data.m
executable file
·222 lines (176 loc) · 11.7 KB
/
PhilipsRead_data.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
function [ MRS_struct ] = PhilipsRead_data(MRS_struct, fname, fname_water )
% RE/CJE Parse SPAR file for header info
% 110825
%ADH 20130701 - modify such that if the water reference scan is loaded in
%the same file, it works. for this case, will only have one input fname
%because all the raw data is in the same place. Next step will be to put
%a flag to see if there is water data in a separate scan, in the same scan
%as a water ref (as called in SDAT) or if its just not been collected at all
% will need to get all the indicies first because will need to know how big
% the waterdata is before loading in the GABAdata.
% work out data header name
sparname = [fname(1:(end-4)) 'list'];
sparheader = textread(sparname, '%s');
%ADH - decide if there is water data as ref data included in the data and
%if so, set a flag to pull it out properly...
sparidx=find(ismember(sparheader, 'number_of_mixes')==1);
nu_mixes = str2num(sparheader{sparidx+2});
% MRS_struct.nu_mixes = nu_mixes;
if nu_mixes == 2
% ADH 20130708 - if have 2 nu_mixes - that mean there is water data
% within the data and need to separate it our properly
MRS_struct.Reference_compound = 'H2O'; % ADH because actually is water, its jsut not in a separate file
sparidx=find(ismember(sparheader, 'F-resolution')==1, 1, 'first'); % gaba_data - ADH 20130701
MRS_struct.npoints = str2num(sparheader{sparidx+2});
sparidx=find(ismember(sparheader, 'F-resolution')==1, 1, 'last'); %water_data- ADH 20130701
MRS_struct.npoints_water = str2num(sparheader{sparidx+2});
sparidx=find(ismember(sparheader, 'number_of_extra_attribute_1_values')==1, 1, 'first');
MRS_struct.nrows = str2num(sparheader{sparidx+2});
sparidx=find(ismember(sparheader, 'number_of_extra_attribute_1_values')==1, 1, 'last');
MRS_struct.nrows_water = str2num(sparheader{sparidx+2}); % this might not be correct because we dont collect all the water data - ADH 20130701
MRS_struct.nrows_water = 1; % ADH 20130701
% set = 1 because only collect 1 water scan but the
% header thinks there will be as many rows as in the
% gaba data
%ADH not completely sure why Navg changes across the structure of data
%reads and others don't but for now keep this characteristic... may
%screw up later though - something to keep an eye out for. - ADH 20130701
sparidx=find(ismember(sparheader, 'number_of_signal_averages')==1, 1, 'first');
MRS_struct.Navg(MRS_struct.ii) = str2num(sparheader{sparidx+2})*MRS_struct.nrows;
sparidx=find(ismember(sparheader, 'number_of_signal_averages')==1, 1, 'last');
MRS_struct.Navg_water(MRS_struct.ii) = str2num(sparheader{sparidx+2}); %Trial SDAT might be average not sum.
%Need to determine the offset of the data - i.e. how many channels are
%there...
sparidx=find(ismember(sparheader, 'NOI')==1); %ADH skip past noise
MRS_struct.coil_chammels=size(sparidx,1)-2;
sparidx=find(ismember(sparheader, 'STD')==1);
%MRS_struct.ptr_offset_water=str2num(sparheader{sparidx(3)+20});
%MRS_struct.ptr_offset_data=str2num(sparheader{sparidx(3+MRS_struct.Navg_water(MRS_struct.ii))+20});
ptr_offset_water=str2num(sparheader{sparidx(3)+20});
ptr_offset_data=str2num(sparheader{sparidx(3+MRS_struct.Navg_water(MRS_struct.ii))+20});
%MRS_struct.sw=2000;
%%%HARDCODED TR - Need to remove in Gannet 2.0
%MRS_struct.TR=2000;
%ADH the actual reading of the data happens with readraw_Gannet - just
% takes in the real data to analyse so need to decide on the offsets
% proior to this. will need to decide it for both the water data and the
% GABA data.
% Use 2 readraw_Gannet - one for water and one for data - need to sort
% out the skips of data...- ADH 20130701
%Need to skip rows associated with the '
MRS_struct.data = readraw_Gannet(fname, 'float', ...
[2 MRS_struct.npoints 1 MRS_struct.Navg(MRS_struct.ii)/MRS_struct.nrows MRS_struct.nrows], 'l',ptr_offset_data);
% Make data complex.
MRS_struct.data = squeeze(MRS_struct.data(1,:,:,:,:)+ 1i*MRS_struct.data(2,:,:,:,:));
MRS_struct.data(5,5)
% then the same for the water data- ADH 20130701
MRS_struct.data_water = readraw_Gannet(fname, 'float', ...
[2 MRS_struct.npoints_water 1 MRS_struct.Navg_water(MRS_struct.ii) MRS_struct.nrows_water], 'l',ptr_offset_water);
% Make data complex.- ADH 20130701
MRS_struct.data_water = squeeze(MRS_struct.data_water(1,:,:,:)+ 1i*MRS_struct.data_water(2,:,:,:));
MRS_struct.data_water(5,5)
% % this part is from the original PhilipsRead_data needs to be
% % removed because not a 2nd waterfile- ADH 20130701
%%'Duplicated' code added in to handle .data water files
% work out water data header name
%
% if nargin >2
% sparnameW = [fname_water(1:(end-4)) 'list'];
% sparheaderW = textread(sparnameW, '%s');
% sparidx=find(ismember(sparheaderW, 'number_of_signal_averages')==1);
% MRS_struct.Nwateravg(MRS_struct.ii) = str2num(sparheaderW{sparidx+2}); %Trial SDAT might be average not sum.
%
%
% %Need to skip rows associated with the '
%
% MRS_struct.data_water = readraw_Gannet(fname_water, 'float', [2 MRS_struct.npoints MRS_struct.Nwateravg(MRS_struct.ii) 1], 'l', MRS_struct.ptr_offset);
%
% % Make data complex.
% MRS_struct.data_water = squeeze(MRS_struct.data_water(1,:,:,:,:)+ 1i*MRS_struct.data_water(2,:,:,:,:));
% end
%
%undo time series phase cycling to match GE
corrph = ones(size(MRS_struct.data));
for jj=1:size(MRS_struct.data,3)
corrph(:,:,jj) = corrph(:,:,jj) * (-1).^(jj+1);
end
MRS_struct.data = MRS_struct.data .* corrph;
%MRS_struct.data = MRS_struct.data .*repmat(conj(MRS_struct.data(1,:))./abs(MRS_struct.data(1,:)),[MRS_struct.npoints 1]);
%Philips data appear to be phased already (ideal case)
MRS_struct.data = conj(MRS_struct.data); %RE 110728 - empirical factor to scale 'like GE'
MRS_struct.data = reshape(MRS_struct.data,[size(MRS_struct.data,1) size(MRS_struct.data,2)*size(MRS_struct.data,3)]);
% Depending on the ordering of OFF and ON, the minus sign here may have
% to be removed.... not sure how to automate/fix... raee 4/9/12
MRS_struct.data = MRS_struct.data .*repmat(conj(MRS_struct.data(1,:))./abs(MRS_struct.data(1,:)),[MRS_struct.npoints 1]);
% if nargin >2 %- ADH 20130701 don't need this if because the Ref data is
% in the same file
MRS_struct.data_water = conj(MRS_struct.data_water); %RE 110728 - empirical factor to scale 'like GE'
MRS_struct.data_water = MRS_struct.data_water .*repmat(conj(MRS_struct.data_water(1,:))./abs(MRS_struct.data_water(1,:)),[MRS_struct.npoints 1]);
% Apply spectral registration to water, use 3-end water scans
[dummy, MRS_struct] = Spectral_Registration(MRS_struct,1);
MRS_struct.tempwater = MRS_struct.data_water;
if(size(MRS_struct.data_water,2)>3)
MRS_struct.data_water = squeeze(mean(MRS_struct.data_water(:,3:end),2)).'; %adh 20130208 sdat is mean of .data - not sum
else
MRS_struct.data_water = squeeze(mean(MRS_struct.data_water(:,end),2)).'; %adh 20130208 sdat is mean of .data - not sum
end
MRS_struct.Nwateravg = 1; %adh 20130208 sdat is mean of .data - not sum - Nwateravg now = to sdat
% not sure if this is needed but will keep
% in for now. - ADH 20130701
% end
else % ADH 20130708 this is the PhilipsRead_data with only 1 mix (so may be separate water files or no water reference
% work out data header name
sparname = [fname(1:(end-4)) 'list'];
sparheader = textread(sparname, '%s');
sparidx=find(ismember(sparheader, 'F-resolution')==1);
MRS_struct.npoints = str2num(sparheader{sparidx+2});
sparidx=find(ismember(sparheader, 'number_of_extra_attribute_1_values')==1);
MRS_struct.nrows = str2num(sparheader{sparidx+2});
sparidx=find(ismember(sparheader, 'number_of_signal_averages')==1);
MRS_struct.Navg(MRS_struct.ii) = str2num(sparheader{sparidx+2});
MRS_struct.Navg(MRS_struct.ii) = MRS_struct.Navg(MRS_struct.ii)*MRS_struct.nrows;
%Need to determine the offset of the data - i.e. how many channels are
%there...
sparidx=find(ismember(sparheader, 'NOI')==1);
MRS_struct.coil_chammels=size(sparidx,1)-2;
sparidx=find(ismember(sparheader, 'STD')==1);
MRS_struct.ptr_offset=str2num(sparheader{sparidx(3)+20});
MRS_struct.sw=2000;
%%%HARDCODED TR - Need to remove in Gannet 2.0
MRS_struct.TR=2000;
%Need to skip rows associated with the '
MRS_struct.data = readraw_Gannet(fname, 'float', [2 MRS_struct.npoints 1 MRS_struct.Navg(MRS_struct.ii)/MRS_struct.nrows MRS_struct.nrows], 'l',MRS_struct.ptr_offset);
% Make data complex.
MRS_struct.data = squeeze(MRS_struct.data(1,:,:,:,:)+ 1i*MRS_struct.data(2,:,:,:,:));
%%'Duplicated' code added in to handle .data water files
% work out water data header name
if nargin >2
sparnameW = [fname_water(1:(end-4)) 'list'];
sparheaderW = textread(sparnameW, '%s');
sparidx=find(ismember(sparheaderW, 'number_of_signal_averages')==1);
MRS_struct.Nwateravg(MRS_struct.ii) = str2num(sparheaderW{sparidx+2}); %Trial SDAT might be average not sum.
%Need to skip rows associated with the '
MRS_struct.data_water = readraw_Gannet(fname_water, 'float', [2 MRS_struct.npoints MRS_struct.Nwateravg(MRS_struct.ii) 1], 'l', MRS_struct.ptr_offset);
% Make data complex.
MRS_struct.data_water = squeeze(MRS_struct.data_water(1,:,:,:,:)+ 1i*MRS_struct.data_water(2,:,:,:,:));
end
%undo time series phase cycling to match GE
corrph = ones(size(MRS_struct.data));
for jj=1:size(MRS_struct.data,3)
corrph(:,:,jj) = corrph(:,:,jj) * (-1).^(jj+1);
end
MRS_struct.data = MRS_struct.data .* corrph;
%MRS_struct.data = MRS_struct.data .*repmat(conj(MRS_struct.data(1,:))./abs(MRS_struct.data(1,:)),[MRS_struct.npoints 1]);
%Philips data appear to be phased already (ideal case)
MRS_struct.data = conj(MRS_struct.data); %RE 110728 - empirical factor to scale 'like GE'
MRS_struct.data = reshape(MRS_struct.data,[size(MRS_struct.data,1) size(MRS_struct.data,2)*size(MRS_struct.data,3)]);
% Depending on the ordering of OFF and ON, the minus sign here may have
% to be removed.... not sure how to automate/fix... raee 4/9/12
MRS_struct.data = MRS_struct.data .*repmat(conj(MRS_struct.data(1,:))./abs(MRS_struct.data(1,:)),[MRS_struct.npoints 1]);
if nargin >2
MRS_struct.data_water = conj(MRS_struct.data_water); %RE 110728 - empirical factor to scale 'like GE'
MRS_struct.data_water = MRS_struct.data_water .*repmat(conj(MRS_struct.data_water(1,:))./abs(MRS_struct.data_water(1,:)),[MRS_struct.npoints 1]);
MRS_struct.data_water = squeeze(mean(MRS_struct.data_water,2)).'; %adh 20130208 sdat is mean of .data - not sum
MRS_struct.Nwateravg = 1; %adh 20130208 sdat is mean of .data - not sum - Nwateravg now = to sdat
end
end