/
Spectral_Registration.m
executable file
·130 lines (97 loc) · 5.98 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
function [AllFramesFTrealign MRS_struct] = Spectral_Registration(MRS_struct, OnWhat)
%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
if(OnWhat) %Read water data
%First, take the complex data and turn it into a real matrix
flatdata(:,1,:)=real(MRS_struct.data_water);
flatdata(:,2,:)=imag(MRS_struct.data_water);
else % read spectro data
%First, take the complex data and turn it into a real matrix
flatdata(:,1,:)=real(MRS_struct.data);
flatdata(:,2,:)=imag(MRS_struct.data);
end
%Correct to a point 10% into the file (seems better that the actual beginning)
AlignRow=ceil(size(flatdata,3)/10);
MRS_struct.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.sw;
time=((0:1:(MRS_struct.npoints-1)).'/MRS_struct.sw);
%Fitting to determine frequency and phase corrections.
for corrloop=1:size(flatdata,3)
target=MRS_struct.flattarget(:);
transient=squeeze(flatdata(:,:,corrloop));
input.data=transient(:);
parsFit(corrloop,:)=nlinfit(input,target,@FreqPhaseShiftNest,parsGuess);
end
if(OnWhat)
%Applyng frequency and phase corrections.
for corrloop=1:size(flatdata,3)
MRS_struct.data_water(:,corrloop)=MRS_struct.data_water(:,corrloop).*exp(1i*parsFit(corrloop,1)*2*pi*time)*exp(1i*pi/180*parsFit(corrloop,2));
end
FullData = MRS_struct.data_water;
FullData = FullData.* repmat( (exp(-(time)*MRS_struct.LB*pi)), [1 size(MRS_struct.data_water,2)]);
AllFramesFTrealign=fftshift(fft(FullData,MRS_struct.ZeroFillTo,1),1);
else
%Applyng frequency and phase corrections.
for corrloop=1:size(flatdata,3)
MRS_struct.data_align(:,corrloop)=MRS_struct.data(:,corrloop).*exp(1i*parsFit(corrloop,1)*2*pi*time)*exp(1i*pi/180*parsFit(corrloop,2));
end
FullData = MRS_struct.data_align;
FullData = FullData.* repmat( (exp(-(time)*MRS_struct.LB*pi)), [1 size(MRS_struct.data,2)]);
AllFramesFTrealign=fftshift(fft(FullData,MRS_struct.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.freq-ChoCrFitLimHigh);
cclb=find(min(z)==z);
z=abs(MRS_struct.freq-ChoCrFitLimLow);
ccub=find(min(z)==z);
freqrange=MRS_struct.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.LarmorFreq) (MRS_struct.LarmorFreq) (180/pi) 1 1 1];
ChoCrMeanSpecFit = FitChoCr(freqrange, ChoCrMeanSpec, ChoCr_initx,MRS_struct.LarmorFreq);
MRS_struct.ChoCrMeanSpecFit = ChoCrMeanSpecFit./[1 (2*MRS_struct.LarmorFreq) (MRS_struct.LarmorFreq) (180/pi) 1 1 1];
% figure(5)
% plot(freqrange,ChoCrMeanSpec,freqrange,TwoLorentzModel(ChoCrMeanSpecFit./[1 (2*MRS_struct.LarmorFreq) (MRS_struct.LarmorFreq) (180/pi) 1 1 1],freqrange))
AllFramesFTrealign=AllFramesFTrealign*exp(1i*pi/180*(ChoCrMeanSpecFit(4)));%phase
ChoCrFreqShift = ChoCrMeanSpecFit(3);
ChoCrFreqShift = ChoCrFreqShift - 3.02*MRS_struct.LarmorFreq;
ChoCrFreqShift = ChoCrFreqShift ./ (MRS_struct.LarmorFreq*(MRS_struct.freq(2) - MRS_struct.freq(1) ));
ChoCrFreqShift_points = round(ChoCrFreqShift)
AllFramesFTrealign=circshift(AllFramesFTrealign, [-ChoCrFreqShift_points 0]);%freq
%MRS_struct.ChoCrFreqShift_pts(ii) = ChoCrFreqShift_points;
%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.freq-CrFitLimHigh);
clb=find(min(z)==z);
z=abs(MRS_struct.freq-CrFitLimLow);
cub=find(min(z)==z);
freqrange=MRS_struct.freq(clb:cub);
Cr_initx = [ Area_estimate Width_estimate 3.02 0 Baseline_offset 0 ].*[1 (2*MRS_struct.LarmorFreq) (MRS_struct.LarmorFreq) (180/pi) 1 1 ];
CrMeanSpec = mean(AllFramesFTrealign(clb:cub,:),2);
CrMeanSpecFit = FitCr(freqrange, CrMeanSpec, Cr_initx);
%Some Output
MRS_struct.FreqStdevHz(MRS_struct.ii)=std(parsFit(:,1),1);
MRS_struct.CrFWHMHz(MRS_struct.ii)=CrMeanSpecFit(2);
end
% AllFramesFTrealign=reshape(AllFramesFTrealign,[size(AllFramesFTrealign,1) MRS_struct.Navg(ii) MRS_struct.nrows]);
%
% OddFramesFTrealign=AllFramesFTrealign(:,:,2:2:end);
% EvenFramesFTrealign=AllFramesFTrealign(:,:,1:2:end);
%
% OddFramesFTrealign=reshape(OddFramesFTrealign,[size(OddFramesFTrealign,1) size(OddFramesFTrealign,2)*size(OddFramesFTrealign,3) ]);
% EvenFramesFTrealign=reshape(EvenFramesFTrealign,[size(EvenFramesFTrealign,1) size(EvenFramesFTrealign,2)*size(EvenFramesFTrealign,3) ]);
%
% MRS_struct.Navg(ii) = MRS_struct.Navg(ii)*MRS_struct.nrows; % ADH fix to have the correct Navg for .data now
end