forked from nputs/Gannet2.0_nicksversion
-
Notifications
You must be signed in to change notification settings - Fork 1
/
Spectral_Registration.m
executable file
·227 lines (190 loc) · 11.9 KB
/
Spectral_Registration.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
function [AllFramesFTrealign MRS_struct] = Spectral_Registration(MRS_struct, OnWhat, Dual)
%Spectral Registration is a time-domain frequency-and-phse correction as
%per Near et al. 2014 [under review].
% OnWhat = 0 for spectro data OnWhat=1 for water data
MRS_struct.p.parsFit=[];
%Dual-channel option only applies registration separately to ONs and OFFs
SpecRegLoop=0;
if(nargin==3)
if(Dual==1)
%We want to run this code twice, once for ONs, once for OFFs.
SpecRegLoop=1;
end
end
while(SpecRegLoop>(-1))
if(OnWhat) %Read water data
%First, take the complex data and turn it into a real matrix
flatdata(:,1,:)=real(MRS_struct.fids.data_water);
flatdata(:,2,:)=imag(MRS_struct.fids.data_water);
else % read spectro data
if(nargin==3)
if(Dual==1)
%This code runs twice, first for ONs, second for OFFs.
SpecRegLoop;
size(real(MRS_struct.fids.data(:,(MRS_struct.fids.ON_OFF==SpecRegLoop))));
flatdata(:,1,:)=real(MRS_struct.fids.data(:,(MRS_struct.fids.ON_OFF==SpecRegLoop)));
flatdata(:,2,:)=imag(MRS_struct.fids.data(:,(MRS_struct.fids.ON_OFF==SpecRegLoop)));
end
else
%First, take the complex data and turn it into a real matrix
flatdata(:,1,:)=real(MRS_struct.fids.data);
flatdata(:,2,:)=imag(MRS_struct.fids.data);
end
end
%Correct to a point 10% into the file (seems better that the actual beginning)
AlignRow=ceil(size(flatdata,3)/10);
MRS_struct.fids.flattarget=squeeze(flatdata(:,:,AlignRow));
%Time domain Frequency and Phase Correction
%Preliminary to fitting:
parsGuess=[0 0]; %initial freq and phase guess
parsFit = zeros([size(flatdata,3) 2]);
input.dwelltime=1/MRS_struct.p.sw;
time=((0:1:(MRS_struct.p.npoints-1)).'/MRS_struct.p.sw);
%Fitting to determine frequency and phase corrections.
for corrloop=1:size(flatdata,3)
target=MRS_struct.fids.flattarget(:);
transient=squeeze(flatdata(:,:,corrloop));
input.data=transient(:);
parsFit(corrloop,:)=nlinfit(input,target,@FreqPhaseShiftNest,parsGuess);
parsGuess = parsFit(corrloop,:); %Carry parameters from point to point
end
if(OnWhat)
%Applyng frequency and phase corrections.
for corrloop=1:size(flatdata,3)
MRS_struct.fids.data_water(:,corrloop)=MRS_struct.fids.data_water(:,corrloop).*exp(1i*parsFit(corrloop,1)*2*pi*time)*exp(1i*pi/180*parsFit(corrloop,2));
end
FullData = MRS_struct.fids.data_water;
FullData = FullData.* repmat( (exp(-(time)*MRS_struct.p.LB*pi)), [1 size(MRS_struct.fids.data_water,2)]);
AllFramesFTrealign=fftshift(fft(FullData,MRS_struct.p.ZeroFillTo,1),1);
else
%Applying frequency and phase corrections.
MRS_struct.p.parsFit=[MRS_struct.p.parsFit parsFit];
for corrloop=1:size(flatdata,3)
if(nargin==3)
if(Dual==1)
%Need to get the slot right to put data back into
averages_per_dynamic=find(MRS_struct.fids.ON_OFF~=(MRS_struct.fids.ON_OFF(1)),1)-1;
dyn=floor((corrloop-1)/averages_per_dynamic); %number of cycles in
ind=mod((corrloop-1),averages_per_dynamic)+1; %number in current cycle
if(SpecRegLoop==1)
if(MRS_struct.fids.ON_OFF(1)==1)
corrloop_d = dyn*averages_per_dynamic*2+ind;
else
corrloop_d = dyn*averages_per_dynamic*2+averages_per_dynamic+ind;
end
else
if(MRS_struct.fids.ON_OFF(1)==1)
corrloop_d = dyn*averages_per_dynamic*2+averages_per_dynamic+ind;
else
corrloop_d = dyn*averages_per_dynamic*2+ind;
end
end
MRS_struct.fids.data_align(:,corrloop_d)=MRS_struct.fids.data(:,corrloop_d).*exp(1i*parsFit(corrloop,1)*2*pi*time)*exp(1i*pi/180*parsFit(corrloop,2));
end
CorrPars(corrloop_d,:)=parsFit(corrloop,:);
else
corrloop_d=corrloop;
MRS_struct.fids.data_align(:,corrloop_d)=MRS_struct.fids.data(:,corrloop_d).*exp(1i*parsFit(corrloop,1)*2*pi*time)*exp(1i*pi/180*parsFit(corrloop,2));
CorrPars(corrloop_d,:)=parsFit(corrloop,:);
end
end
if(SpecRegLoop==0)
FullData = MRS_struct.fids.data_align;
FullData = FullData.* repmat( (exp(-(time)*MRS_struct.p.LB*pi)), [1 size(MRS_struct.fids.data,2)]);
AllFramesFTrealign=fftshift(fft(FullData,MRS_struct.p.ZeroFillTo,1),1);
%In FD, move Cr to 3.02 and get phase 'right' as opposed to 'consistent'.
ChoCrFitLimLow=2.6;
ChoCrFitLimHigh=3.6;
%Still need ranges for Creatine align plot
z=abs(MRS_struct.spec.freq-ChoCrFitLimHigh);
cclb=find(min(z)==z);
z=abs(MRS_struct.spec.freq-ChoCrFitLimLow);
ccub=find(min(z)==z);
freqrange=MRS_struct.spec.freq(cclb:ccub);
%Do some detective work to figure out the initial parameters
ChoCrMeanSpec = mean(AllFramesFTrealign(cclb:ccub,:),2);
Baseline_offset=real(ChoCrMeanSpec(1)+ChoCrMeanSpec(end))/2;
Width_estimate=0.05;%ppm
Area_estimate=(max(real(ChoCrMeanSpec))-min(real(ChoCrMeanSpec)))*Width_estimate*4;
ChoCr_initx = [ Area_estimate Width_estimate 3.02 0 Baseline_offset 0 1].*[1 (2*MRS_struct.p.LarmorFreq) (MRS_struct.p.LarmorFreq) (180/pi) 1 1 1];
if(nargin==3)
if(Dual==1)
%This bit is silly - we don't want to do OFF-to-ON based on the CR signal
ChoCrMeanSpecON = mean(AllFramesFTrealign(cclb:ccub,(MRS_struct.fids.ON_OFF==1)),2);
ChoCrMeanSpecOFF = mean(AllFramesFTrealign(cclb:ccub,(MRS_struct.fids.ON_OFF==0)),2);
ChoCrMeanSpecFitON = FitChoCr(freqrange, ChoCrMeanSpecON, ChoCr_initx,MRS_struct.p.LarmorFreq);
ChoCrMeanSpecFitOFF = FitChoCr(freqrange, ChoCrMeanSpecOFF, ChoCr_initx,MRS_struct.p.LarmorFreq);
AllFramesFTrealign(:,(MRS_struct.fids.ON_OFF==1))=AllFramesFTrealign(:,(MRS_struct.fids.ON_OFF==1))*exp(1i*pi/180*(ChoCrMeanSpecFitON(4)));%phase
AllFramesFTrealign(:,(MRS_struct.fids.ON_OFF==0))=AllFramesFTrealign(:,(MRS_struct.fids.ON_OFF==0))*exp(1i*pi/180*(ChoCrMeanSpecFitOFF(4)));%phase
ChoCrFreqShiftON = ChoCrMeanSpecFitON(3);
ChoCrFreqShiftON = ChoCrFreqShiftON - 3.02*MRS_struct.p.LarmorFreq;
ChoCrFreqShiftON = ChoCrFreqShiftON ./ (MRS_struct.p.LarmorFreq*(MRS_struct.spec.freq(2) - MRS_struct.spec.freq(1) ));
ChoCrFreqShift_pointsON = round(ChoCrFreqShiftON);
AllFramesFTrealign(:,(MRS_struct.fids.ON_OFF==1))=circshift(AllFramesFTrealign(:,(MRS_struct.fids.ON_OFF==1)), [-ChoCrFreqShift_pointsON 0]);%freq
ChoCrFreqShiftOFF = ChoCrMeanSpecFitOFF(3);
ChoCrFreqShiftOFF = ChoCrFreqShiftOFF - 3.02*MRS_struct.p.LarmorFreq;
ChoCrFreqShiftOFF = ChoCrFreqShiftOFF ./ (MRS_struct.p.LarmorFreq*(MRS_struct.spec.freq(2) - MRS_struct.spec.freq(1) ));
ChoCrFreqShift_pointsOFF = round(ChoCrFreqShiftOFF);
AllFramesFTrealign(:,(MRS_struct.fids.ON_OFF==0))=circshift(AllFramesFTrealign(:,(MRS_struct.fids.ON_OFF==0)), [-ChoCrFreqShift_pointsOFF 0]);%freq
end
MRS_struct.out.FreqStdevHz(MRS_struct.ii)=std(parsFit(:,1),1);
MRS_struct.out.CrFWHMHz(MRS_struct.ii)=mean([ChoCrMeanSpecFitON(2) ChoCrMeanSpecFitOFF(2)]);
else
ChoCrMeanSpecFit = FitChoCr(freqrange, ChoCrMeanSpec, ChoCr_initx,MRS_struct.p.LarmorFreq);
MRS_struct.out.ChoCrMeanSpecFit = ChoCrMeanSpecFit./[1 (2*MRS_struct.p.LarmorFreq) (MRS_struct.p.LarmorFreq) (180/pi) 1 1 1];
AllFramesFTrealign=AllFramesFTrealign*exp(1i*pi/180*(ChoCrMeanSpecFit(4)));%phase
ChoCrFreqShift = ChoCrMeanSpecFit(3);
ChoCrFreqShift = ChoCrFreqShift - 3.02*MRS_struct.p.LarmorFreq;
ChoCrFreqShift = ChoCrFreqShift ./ (MRS_struct.p.LarmorFreq*(MRS_struct.spec.freq(2) - MRS_struct.spec.freq(1) ));
ChoCrFreqShift_points = round(ChoCrFreqShift);
AllFramesFTrealign=circshift(AllFramesFTrealign, [-ChoCrFreqShift_points 0]);%freq
end
% figure(5)
% plot(freqrange,ChoCrMeanSpec,freqrange,TwoLorentzModel(ChoCrMeanSpecFit./[1 (2*MRS_struct.p.LarmorFreq) (MRS_struct.p.LarmorFreq) (180/pi) 1 1 1],freqrange))
%Fit just the Cr in the aligned mean spectrum to get CrFWHMHz
CrFitLimLow=2.6;
CrFitLimHigh=3.11;
%Still need ranges for Creatine align plot
z=abs(MRS_struct.spec.freq-CrFitLimHigh);
clb=find(min(z)==z);
z=abs(MRS_struct.spec.freq-CrFitLimLow);
cub=find(min(z)==z);
freqrange=MRS_struct.spec.freq(clb:cub);
Cr_initx = [ Area_estimate Width_estimate 3.02 0 Baseline_offset 0 ].*[1 (2*MRS_struct.p.LarmorFreq) (MRS_struct.p.LarmorFreq) (180/pi) 1 1 ];
CrMeanSpec = mean(AllFramesFTrealign(clb:cub,:),2);
CrMeanSpecFit = FitCr(freqrange, CrMeanSpec, Cr_initx);
%Some Output
MRS_struct.out.FreqStdevHz(MRS_struct.ii)=std(parsFit(:,1),1);
% GO 01/29/16: output frequency correction in Hz for every
% average
MRS_struct.out.FreqHz(MRS_struct.ii,:)=parsFit(:,1).';
MRS_struct.out.CrFWHMHz(MRS_struct.ii)=CrMeanSpecFit(2);
%Decide which rows to reject based on 3-sigma
% Reject any point where the fit params - phase
% or freq are > 3stdev away from the mean
MeanFrameParams = mean(CorrPars, 1);
UpperLim = repmat(MeanFrameParams + 3*std(CorrPars,1), [size(AllFramesFTrealign,2) 1]);
LowerLim = repmat(MeanFrameParams - 3*std(CorrPars,1), [size(AllFramesFTrealign,2) 1]);
%but don't reject on linear, const baseline fit vals
rejectframe = gt(CorrPars, UpperLim);
rejectframe = rejectframe + lt(CorrPars, LowerLim);
MRS_struct.out.reject(:,MRS_struct.ii) = max(rejectframe,[],2);
%Balance up rejects
IsRejectOnOrOff=MRS_struct.fids.ON_OFF(MRS_struct.out.reject(:,MRS_struct.ii)==1);
stepsize=find(MRS_struct.fids.ON_OFF~=(MRS_struct.fids.ON_OFF(1)),1)-1;
for kk=1:size(MRS_struct.out.reject(:,MRS_struct.ii),1)
%first find whether reject is ON or OFF
if(MRS_struct.out.reject(kk,MRS_struct.ii)==1)
IsRejectOnOrOff=MRS_struct.fids.ON_OFF(kk);
if IsRejectOnOrOff==MRS_struct.fids.ON_OFF(1)
MRS_struct.out.reject(kk+stepsize,MRS_struct.ii)=1;
else
MRS_struct.out.reject(kk-stepsize,MRS_struct.ii)=1;
end
end
end
end
end
SpecRegLoop=SpecRegLoop-1;
end
end