diff --git a/MM_T1processing/plot_mean_std_of_spectra_MMT1.m b/MM_T1processing/plot_mean_std_of_spectra_MMT1.m
index 14fbe603c6d433b6e22e9b1584066c42a23c7f99..26d7983bffc3013bf3a9dbb2688e94104142e394 100644
--- a/MM_T1processing/plot_mean_std_of_spectra_MMT1.m
+++ b/MM_T1processing/plot_mean_std_of_spectra_MMT1.m
@@ -3,9 +3,11 @@ function plot_mean_std_of_spectra_MM
pathName = 'MM data path';
[path] = pathToDataFolder(pathName);
-fileNameBase = [path, 'Summed_spectra_TE'];
+fileNameBase = [path, 'Summed_spectra_'];
TEs = [24; 32; 40; 52; 60];
+TI1 = [];
+TI2 = [];
colors = {[0 0.7 0]; [1 0 0]; [0 0 1]; [0 0 0]; [0 0.5 0.5];};
numberOfTEs = length(TEs);
diff --git a/MM_T1processing/plot_mean_std_of_spectra_T1.m b/MM_T1processing/plot_mean_std_of_spectra_T1.m
new file mode 100644
index 0000000000000000000000000000000000000000..f0ad46b3608b789e19fd73309db157b347166e1c
--- /dev/null
+++ b/MM_T1processing/plot_mean_std_of_spectra_T1.m
@@ -0,0 +1,172 @@
+function plot_mean_std_of_spectra_MM
+
+downField_MM_Met = 'MM';
+switch downField_MM_Met
+ case 'DF'
+ pathName = 'DF data path';
+ TEs = [24; 24; 32; 40; 52; 60; ];
+ colors = {[0 0.7 0]; [0 0.7 0]; [0 0 0]; [0 0 1]; [1 0 0]; [0 0.5 0.5];};
+ offsetStep = 2.2e-5;
+ case 'MM'
+ pathName = 'MM T1 data path';
+ TEs = [1; 2; 3; 4; 5; 6; 7; 8; 9; 10; 11; ];
+ colors = {[76 187 23]/255; [76 187 23]/255; [76 187 23]/255; [76 187 23]/255; [76 187 23]/255; [1 121 111]/255; [75 0 130]/255; [75 0 130]/255; [75 0 130]/255; [75 0 130]/255; [75 0 130]/255;};
+ offsetStep = 7.0e-5;
+ case 'MM WM'
+ pathName = 'MM WM T1 data path';
+ TEs = [1; 2; 3; 4; 5; 6; 7; 8; 9; 10; 11; ];
+ colors = {[76 187 23]/255; [76 187 23]/255; [76 187 23]/255; [76 187 23]/255; [76 187 23]/255; [1 121 111]/255; [75 0 130]/255; [75 0 130]/255; [75 0 130]/255; [75 0 130]/255; [75 0 130]/255;};
+ offsetStep = 7.0e-5;
+ case 'Met'
+ pathName = 'DF data path';
+ TEs = [24; 32; 40; 52; 60; ];
+ colors = {[0 0.7 0]; [0 0 0]; [0 0 1]; [1 0 0]; [0 0.5 0.5];};
+ offsetStep = 2.3e-4;
+end
+
+[path] = pathToDataFolder(pathName);
+
+fileNameBase = [path, 'Summed_Averaged_TI1_'];
+
+
+numberOfTEs = length(TEs);
+
+offSetBase = offsetStep * (numberOfTEs + 7); %+7 for WM
+
+for indexTE = 1:numberOfTEs
+ fileName = [fileNameBase num2str(TEs(indexTE))];
+ offset = offSetBase - indexTE * offsetStep;
+ if (indexTE == 1) && (strcmp(downField_MM_Met,'DF'))
+ plotOneTE(fileName, TEs(indexTE), colors{indexTE}, offset+offsetStep/2, 'DF_up');
+ else
+ plotOneTE(fileName, TEs(indexTE), colors{indexTE}, offset, downField_MM_Met);
+ end
+end
+
+end
+
+function plotOneTE(fileName, TE, colorSpectrum, offset, downField_MM_Met)
+load([fileName '.mat'], 'summedSpectraTI')
+
+zeroPPM = 4.67;
+nTimepoints = size(summedSpectraTI.Data{1},summedSpectraTI.kx_dim);
+nAverages = size(summedSpectraTI.Data{1},summedSpectraTI.meas_dim);
+waterFreq = summedSpectraTI.Parameter.Headers.ScanFrequency * 1E-6;
+bandwidth = summedSpectraTI.Parameter.Headers.Bandwidth_Hz;
+frequencyRange = linspace(-bandwidth/2,bandwidth/2,nTimepoints);
+
+if strcmp(downField_MM_Met,'DF_up')
+ fidSpectra = summedSpectraTI.Data{1}(:,1,1,1,1,1,1,1,1,1,1,:) * 5;
+else
+ fidSpectra = summedSpectraTI.Data{1}(:,1,1,1,1,1,1,1,1,1,1,:) * 1;
+end
+
+fidWater = summedSpectraTI.Data{1}(:,1,1,1,1,1,1,1,2,1,1,:);
+spectra = squeeze(fftshift(fft(fidSpectra)));
+water = squeeze(fftshift(fft(fidWater)));
+meanSpectrum = mean(spectra,2);
+stdSpectrum = std(spectra,[],2);
+relativeStd = stdSpectrum ./ max(meanSpectrum) * 100;
+
+ppmVector = (frequencyRange)/waterFreq+zeroPPM;
+
+% figure(101)
+% plot(ppmVector,relativeStd)
+% xlim([0.5 4.2]);
+% set(gca,'xDir','reverse')
+
+
+%%
+scalingRange = (ppmVector < 1.8) & (ppmVector > 0);
+for indexAverage = 1:nAverages
+ [maxWater,indexMaxWater] = max(abs(real(water(:, indexAverage))));
+
+ [~, indexMaxWaterInRange] = max(abs(real((water(scalingRange, indexAverage)-mean(water(scalingRange, indexAverage))))));
+
+ [maxWaterInRange, ~] = max(real(water(scalingRange, indexAverage)));
+ minWaterInRange = min(real(water(scalingRange, indexAverage)));
+
+ maxSpectraInRange = max(real(spectra(scalingRange, indexAverage)));
+ minSpectraInRange = min(real(spectra(scalingRange, indexAverage)));
+
+ indexStart = find(scalingRange);
+ indexMaxWaterInRange = indexMaxWaterInRange + indexStart(1);
+ referenceLines = zeros(nTimepoints,1);
+ referenceLines(:) = -maxWaterInRange*2;
+ referenceLines(indexMaxWaterInRange:indexMaxWater + (indexMaxWater - indexMaxWaterInRange)) = maxWaterInRange*2;
+
+ plotOffset = maxSpectraInRange;
+% figWaterSpectra = figure(2);
+% hold on
+% plot(ppmScale, real(water(:,indexAverage)) - minWaterInRange + plotOffset,'b')
+% plotScale = (maxWaterInRange-minWaterInRange);
+% ylim([minSpectraInRange-(plotOffset*0.3) plotScale*1.2+plotOffset]);
+% hold on
+% pSpectra = plot(ppmScale, real(spectra(:,indexAverage)),'r');
+% hold on
+% % figure
+% plot(ppmScale, referenceLines, 'k--')
+% xlim([0 8.5]);
+% set(gca,'xDir','reverse')
+% set(gca,'LineWidth',1);
+% % set(gca,'ytick',[])
+% xlabel('[ppm]');
+% makedatatip(pSpectra,[indexMaxWaterInRange])
+%
+% print(['WaterSpectraComparisonPlots\' num2str(indexAverage) '_' fileName],'-dpng')
+% close(figWaterSpectra)
+end
+%%
+% figure
+% plot(ppmScale, real(spectra))
+
+% xlim([0.5 4.2]);
+% set(gca,'xDir','reverse')
+% set(gca,'LineWidth',1);
+% % set(gca,'ytick',[])
+% xlabel('[ppm]');
+% title(['Individual Spectra TE ' num2str(TE) ' ms'])
+
+%%
+figure(1)
+stdAlpha = 0.4;
+hold on
+if strcmp(downField_MM_Met,'DF_up')
+ indexDFCut = find(ppmVector > 8.6, 1, 'first');
+ ppmVector = ppmVector(indexDFCut:end);
+ meanSpectrum = meanSpectrum(indexDFCut:end);
+ stdSpectrum = stdSpectrum(indexDFCut:end);
+end
+
+h = area(ppmVector,[(meanSpectrum - stdSpectrum) + offset (stdSpectrum * 2)]);
+
+h(1).FaceColor = 'none';
+h(2).FaceColor = colorSpectrum;
+h(2).FaceAlpha = stdAlpha;
+h(1).EdgeColor = 'none';
+h(2).EdgeColor = 'none';
+hold on
+pSpectrum = plot(ppmVector,meanSpectrum + offset, 'LineWidth',2, 'Color', colorSpectrum);
+
+switch downField_MM_Met
+ case 'DF'
+ xlim([5.5 9.4]);
+ case 'DF_up'
+ xlim([5.5 9.4]);
+ case 'MM'
+ xlim([0.5 4.2]);
+% ylim([4.5e-6 4.7e-5]);
+ case 'MM WM'
+ xlim([0.5 4.2]);
+ %ylim([2e-6 1.3e-4]); %TODO: check why scaling is different, then MM
+ case 'Met'
+ xlim([0.5 4.2]);
+ ylim([7e-5 1.65e-3]);
+end
+set(gca,'xDir','reverse')
+set(gca,'LineWidth',1);
+set(gca,'ytick',[])
+xlabel('[ppm]');
+title('Spectra mean and standard deviation')
+% title(['Spectra TE ' num2str(TE) ' ms'])
+end
\ No newline at end of file
diff --git a/MM_T1processing/reconstruct_all_MMT1.m b/MM_T1processing/reconstruct_all_MMT1.m
index 3c679cc14b7e25594c4552e89ea11dbb27206525..969d64f796a833cc1b1a89739c5dec4f14c85162 100644
--- a/MM_T1processing/reconstruct_all_MMT1.m
+++ b/MM_T1processing/reconstruct_all_MMT1.m
@@ -1,19 +1,18 @@
function reconstruct_all_MM()
-pathName = 'MM T1 data path';
+pathName = 'MM WM T1 data path';
subjects = { ...
-% '1717' ...
-% '2016' ...
-% '9810' ...
-% '3490' ...
-% '2017' ...
- '7338'...
-% '1658'...
-% '2349'...
-% '7782'...
-% '3373'...
-% '2020'...
- };
+ '1717' ...
+ '2016' ...
+ '9810' ...
+ '3490' ...
+ '2017' ...
+ '7338'...
+ '1658'...
+ '2349'...
+ '2020'...
+ }; % '7782'...
+ % '3373'...
pathBase = pathToDataFolder(pathName, subjects);
diff --git a/MM_processing/plot_mean_std_of_spectra_MM.m b/MM_processing/plot_mean_std_of_spectra_MM.m
index f06457ca18d6b6413ad64b2b3b2002ac0b58f73f..bfbec1a23982461390aeb6dbe2f5c6ddf8061932 100644
--- a/MM_processing/plot_mean_std_of_spectra_MM.m
+++ b/MM_processing/plot_mean_std_of_spectra_MM.m
@@ -1,6 +1,6 @@
function plot_mean_std_of_spectra_MM
-downField_MM_Met = 'MM WM';
+downField_MM_Met = 'MM';
switch downField_MM_Met
case 'DF'
pathName = 'DF data path';
diff --git a/MM_processing/reconstruct_MM.m b/MM_processing/reconstruct_MM.m
index 2f7fd4d7dbc5c75d197079fbe175a20d9b880dd0..581135c80518b902a3735273e3c912eb1aea6d63 100644
--- a/MM_processing/reconstruct_MM.m
+++ b/MM_processing/reconstruct_MM.m
@@ -5,9 +5,9 @@ end
%%
a = MR_spectroS(fid_id,1);
%%
-trunc_ms = 250;
-truncPoint = floor( trunc_ms/(a.Parameter.Headers.DwellTimeSig_ns*1e-6));
-a = a.Truncate(truncPoint);
+% trunc_ms = 250;
+% truncPoint = floor( trunc_ms/(a.Parameter.Headers.DwellTimeSig_ns*1e-6));
+% a = a.Truncate(truncPoint);
a = a.FrequencyAlign;
a = a.DeleteMovedAverages;
@@ -22,10 +22,10 @@ a.Parameter.FreqAlignFreqDomainSettings.selectedMix = 1;
%indices of Interest: given in ppm
% water peak
-a.Parameter.FreqAlignFreqDomainSettings.peaksInPpm = [3.925];
+% a.Parameter.FreqAlignFreqDomainSettings.peaksInPpm = [3.925];
%set zeroFillingParameter to get smooth approximations
-a.Parameter.FreqAlignFreqDomainSettings.zeroFillFactor = 50;
+% a.Parameter.FreqAlignFreqDomainSettings.zeroFillFactor = 50;
%search area (+-) in ppm
a.Parameter.FreqAlignFreqDomainSettings.searchArea = 0.1;
@@ -37,20 +37,20 @@ a.Parameter.FreqAlignFreqDomainSettings.splineSmoothingCoeff = 0.01;
% do the actual Frequency Alignment
a = a.FrequencyAlignFreqDomain;
%%
-a.Parameter.HsvdSettings.bound = [-160 250];
-a.Parameter.HsvdSettings.n = size(a.Data{1},1);
-a.Parameter.HsvdSettings.p = 25;
-[~, a] = a.Hsvd;
+% a.Parameter.HsvdSettings.bound = [-160 160];
+% a.Parameter.HsvdSettings.n = 4096; %size(a.Data{1},1);
+% a.Parameter.HsvdSettings.p = 25;
+% [~, a] = a.Hsvd;
+% %%
+% if(doReverse)
+% a = a.ReverseMC;
+% a.Parameter.HsvdSettings.bound = [-160 160];
+% [~, a] = a.Hsvd;
+% end
%%
-if(doReverse)
- a = a.ReverseMC;
- a.Parameter.HsvdSettings.bound = [-160 160];
- [~, a] = a.Hsvd;
-end
-%%
-trunc_ms = 150;
-truncPoint = floor( trunc_ms/(a.Parameter.Headers.DwellTimeSig_ns*1e-6));
-a = a.Truncate(truncPoint);
+% trunc_ms = 150;
+% truncPoint = floor( trunc_ms/(a.Parameter.Headers.DwellTimeSig_ns*1e-6));
+% a = a.Truncate(truncPoint);
%% Phase correction
if ~exist('phase0','var')
phase0 = 0;
diff --git a/Met_T1processing/plot_mean_std_of_spectra_metT1.m b/Met_T1processing/plot_mean_std_of_spectra_metT1.m
new file mode 100644
index 0000000000000000000000000000000000000000..18b35c7933b685c5e26cd378cd6969a01d2c3280
--- /dev/null
+++ b/Met_T1processing/plot_mean_std_of_spectra_metT1.m
@@ -0,0 +1,172 @@
+function plot_mean_std_of_spectra_metT1
+
+downField_MM_Met = 'Met';
+switch downField_MM_Met
+ case 'DF'
+ pathName = 'DF data path';
+ TEs = [24; 24; 32; 40; 52; 60; ];
+ colors = {[0 0.7 0]; [0 0.7 0]; [0 0 0]; [0 0 1]; [1 0 0]; [0 0.5 0.5];};
+ offsetStep = 2.2e-5;
+ case 'MM'
+ pathName = 'MM T1 data path';
+ TEs = [1; 2; 3; 4; 5; 6; 7; 8; 9; 10; 11; ];
+ colors = {[76 187 23]/255; [76 187 23]/255; [1 121 111]/255; [75 0 130]/255; [75 0 130]/255; [75 0 130]/255; [75 0 130]/255;};
+ offsetStep = 7.0e-5;
+ case 'MM WM'
+ pathName = 'MM WM T1 data path';
+ TEs = [1; 2; 3; 4; 5; 6; 7; 8; 9; 10; 11; ];
+ colors = {[76 187 23]/255; [76 187 23]/255; [76 187 23]/255; [76 187 23]/255; [76 187 23]/255; [1 121 111]/255; [75 0 130]/255; [75 0 130]/255; [75 0 130]/255; [75 0 130]/255; [75 0 130]/255;};
+ offsetStep = 7.0e-5;
+ case 'Met'
+ pathName = 'Met T1 WM data path';
+ TEs = [2500; 1400; 1000; 700; 400; 100; 20; ];
+ colors = {[76 187 23]/255; [76 187 23]/255; [1 121 111]/255; [75 0 130]/255; [75 0 130]/255; [75 0 130]/255; [75 0 130]/255;};%{[0 0.7 0]; [0 0 0]; [0 0 1]; [1 0 0]; [0 0.5 0.5];};
+ offsetStep = 40e-4;
+end
+
+[path] = pathToDataFolder(pathName);
+
+fileNameBase = [path, 'Summed_Averaged_TE'];
+
+
+numberOfTEs = length(TEs);
+
+offSetBase = offsetStep * (numberOfTEs); %+7 for WM
+
+for indexTE = 1:numberOfTEs
+ fileName = [fileNameBase num2str(TEs(indexTE))];
+ offset = offSetBase - indexTE * offsetStep;
+ if (indexTE == 1) && (strcmp(downField_MM_Met,'DF'))
+ plotOneTE(fileName, TEs(indexTE), colors{indexTE}, offset+offsetStep/2, 'DF_up');
+ else
+ plotOneTE(fileName, TEs(indexTE), colors{indexTE}, offset, downField_MM_Met);
+ end
+end
+
+end
+
+function plotOneTE(fileName, TE, colorSpectrum, offset, downField_MM_Met)
+load([fileName '.mat'], 'summedSpectraTE')
+
+zeroPPM = 4.67;
+nTimepoints = size(summedSpectraTE.Data{1},summedSpectraTE.kx_dim);
+nAverages = size(summedSpectraTE.Data{1},summedSpectraTE.meas_dim);
+waterFreq = summedSpectraTE.Parameter.Headers.ScanFrequency * 1E-6;
+bandwidth = summedSpectraTE.Parameter.Headers.Bandwidth_Hz;
+frequencyRange = linspace(-bandwidth/2,bandwidth/2,nTimepoints);
+
+if strcmp(downField_MM_Met,'DF_up')
+ fidSpectra = summedSpectraTE.Data{1}(:,1,1,1,1,1,1,1,1,1,1,:) * 5;
+else
+ fidSpectra = summedSpectraTE.Data{1}(:,1,1,1,1,1,1,1,1,1,1,:) * 1;
+end
+
+fidWater = summedSpectraTE.Data{1}(:,1,1,1,1,1,1,1,2,1,1,:);
+spectra = squeeze(fftshift(fft(fidSpectra)));
+water = squeeze(fftshift(fft(fidWater)));
+meanSpectrum = mean(spectra,2);
+stdSpectrum = std(spectra,[],2);
+relativeStd = stdSpectrum ./ max(meanSpectrum) * 100;
+
+ppmVector = (frequencyRange)/waterFreq+zeroPPM;
+
+% figure(101)
+% plot(ppmVector,relativeStd)
+% xlim([0.5 4.2]);
+% set(gca,'xDir','reverse')
+
+
+%%
+scalingRange = (ppmVector < 1.8) & (ppmVector > 0);
+for indexAverage = 1:nAverages
+ [maxWater,indexMaxWater] = max(abs(real(water(:, indexAverage))));
+
+ [~, indexMaxWaterInRange] = max(abs(real((water(scalingRange, indexAverage)-mean(water(scalingRange, indexAverage))))));
+
+ [maxWaterInRange, ~] = max(real(water(scalingRange, indexAverage)));
+ minWaterInRange = min(real(water(scalingRange, indexAverage)));
+
+ maxSpectraInRange = max(real(spectra(scalingRange, indexAverage)));
+ minSpectraInRange = min(real(spectra(scalingRange, indexAverage)));
+
+ indexStart = find(scalingRange);
+ indexMaxWaterInRange = indexMaxWaterInRange + indexStart(1);
+ referenceLines = zeros(nTimepoints,1);
+ referenceLines(:) = -maxWaterInRange*2;
+ referenceLines(indexMaxWaterInRange:indexMaxWater + (indexMaxWater - indexMaxWaterInRange)) = maxWaterInRange*2;
+
+ plotOffset = maxSpectraInRange;
+% figWaterSpectra = figure(2);
+% hold on
+% plot(ppmScale, real(water(:,indexAverage)) - minWaterInRange + plotOffset,'b')
+% plotScale = (maxWaterInRange-minWaterInRange);
+% ylim([minSpectraInRange-(plotOffset*0.3) plotScale*1.2+plotOffset]);
+% hold on
+% pSpectra = plot(ppmScale, real(spectra(:,indexAverage)),'r');
+% hold on
+% % figure
+% plot(ppmScale, referenceLines, 'k--')
+% xlim([0 8.5]);
+% set(gca,'xDir','reverse')
+% set(gca,'LineWidth',1);
+% % set(gca,'ytick',[])
+% xlabel('[ppm]');
+% makedatatip(pSpectra,[indexMaxWaterInRange])
+%
+% print(['WaterSpectraComparisonPlots\' num2str(indexAverage) '_' fileName],'-dpng')
+% close(figWaterSpectra)
+end
+%%
+% figure
+% plot(ppmScale, real(spectra))
+
+% xlim([0.5 4.2]);
+% set(gca,'xDir','reverse')
+% set(gca,'LineWidth',1);
+% % set(gca,'ytick',[])
+% xlabel('[ppm]');
+% title(['Individual Spectra TE ' num2str(TE) ' ms'])
+
+%%
+figure(1)
+stdAlpha = 0.4;
+hold on
+if strcmp(downField_MM_Met,'DF_up')
+ indexDFCut = find(ppmVector > 8.6, 1, 'first');
+ ppmVector = ppmVector(indexDFCut:end);
+ meanSpectrum = meanSpectrum(indexDFCut:end);
+% stdSpectrum = stdSpectrum(indexDFCut:end);
+end
+
+% h = area(ppmVector,[(meanSpectrum - stdSpectrum) + offset (stdSpectrum * 2)]);
+
+h(1).FaceColor = 'none';
+h(2).FaceColor = colorSpectrum;
+h(2).FaceAlpha = stdAlpha;
+h(1).EdgeColor = 'none';
+h(2).EdgeColor = 'none';
+hold on
+pSpectrum = plot(ppmVector,meanSpectrum + offset, 'LineWidth',2, 'Color', colorSpectrum);
+
+switch downField_MM_Met
+ case 'DF'
+ xlim([5.5 9.4]);
+ case 'DF_up'
+ xlim([5.5 9.4]);
+ case 'MM'
+ xlim([0.5 4.2]);
+% ylim([4.5e-6 4.7e-5]);
+ case 'MM WM'
+ xlim([0.5 4.2]);
+ %ylim([2e-6 1.3e-4]); %TODO: check why scaling is different, then MM
+ case 'Met'
+ xlim([0.5 4.2]);
+% ylim([5000 5000]);
+end
+set(gca,'xDir','reverse')
+set(gca,'LineWidth',1);
+set(gca,'ytick',[])
+xlabel('[ppm]');
+title('Spectra mean and standard deviation')
+% title(['Spectra TE ' num2str(TE) ' ms'])
+end
\ No newline at end of file
diff --git a/Met_T1processing/reconstruct_all_MetT1.m b/Met_T1processing/reconstruct_all_MetT1.m
index bb75c2afe9a798584e86231be05dd5496465c9d5..9d271c4b9b6bdaf4236f857bfb516a1a7775fb8e 100644
--- a/Met_T1processing/reconstruct_all_MetT1.m
+++ b/Met_T1processing/reconstruct_all_MetT1.m
@@ -1,6 +1,6 @@
function reconstruct_all_MetT1()
-pathName = 'MET_T1 data path';
+pathName = 'Met T1 WM data path';
subjects = { ...
'1405' ...
'1706' ...
@@ -18,7 +18,7 @@ subjects = { ...
pathBase = pathToDataFolder(pathName, subjects);
TIs = [20; 100; 400; 700; 1000; 1400; 2500];
-
+%TIs = [1000];
numberOfSubjects = length(subjects);
paths = cell(1,numberOfSubjects);
for indexSubj = 1:numberOfSubjects
@@ -48,7 +48,7 @@ addSinglet0ppm = 'Yes';
numberOfTIs = length(TIs);
filesPath = cell(numberOfTIs,numberOfSubjects);
-waterFilesPath = cell(1,numberOfSubjects);
+% waterFilesPath = cell(1,numberOfSubjects);
files = [files_TI20; files_TI100; files_TI400; files_TI700; files_TI1000; files_TI1400; files_TI2500];
for indexSubj = 1:numberOfSubjects
for indexTI = 1:numberOfTIs
diff --git a/algorithms/@MR_spectroS/EddyCurrentCorrection.m b/algorithms/@MR_spectroS/EddyCurrentCorrection.m
index 8a7952288d9154c2f8533eec805c67c0711b442e..9b0cc36a2eb9cdba3319508dae390e5bd45a5fb9 100644
--- a/algorithms/@MR_spectroS/EddyCurrentCorrection.m
+++ b/algorithms/@MR_spectroS/EddyCurrentCorrection.m
@@ -85,9 +85,9 @@ end
nBlocks = this.Parameter.AverageSettings.nBlocks;
% average the data based on the number of blocks
nAverages = size(waterReferenceData, this.meas_dim);
- if mod(nAverages, nBlocks)~=0
- error('The number of measurements of the water reference should be multiple of the number of blocks')
- end
+% if mod(nAverages, nBlocks)~=0
+% error('The number of measurements of the water reference should be multiple of the number of blocks')
+% end
averagesPerBlock = nAverages/nBlocks;
for indexOfBlocks=1:nBlocks
@@ -98,7 +98,7 @@ end
%% TODO: add combine coils with the weights - one has to add the weights as a parameter
%% Apply the actual Eddy Current Correction using the phase of the water reference (Klose correction)
- data = data.*exp(-1j*angle(averagedWaterReferenceData));
+ data = data.*exp(-1*j*angle(averagedWaterReferenceData)); %
end
this.Data{1} = data;
this.Parameter.ReconFlags.isEddyCurrentCorrected= true;
diff --git a/algorithms/@MR_spectroS/Sorting.m b/algorithms/@MR_spectroS/Sorting.m
index fe41cb7c8d546008ba16a8baffd481e537a2d42d..b1c5f4ad7fb243867819f7252d6e4c96640882a2 100644
--- a/algorithms/@MR_spectroS/Sorting.m
+++ b/algorithms/@MR_spectroS/Sorting.m
@@ -9,6 +9,7 @@ else
nAverages = size(this.Data{1},this.meas_dim);
newDataSize = size(this.Data{1});
newDataSize(this.dyn_dim) = 1;
+ newDataSize(this.extr2_dim) = 1;
newDataSize(this.meas_dim) = NDyns*nAverages;
newData = zeros(newDataSize);
oldData = this.Data{1};
diff --git a/postprocessing/calculate_FWHM.m b/postprocessing/calculate_FWHM.m
index 7b51f6eb4849f2cb48cb3a6587ae09e1d26fbd71..d0d3633a039b1e5bc0cad0e0371c97492db4d2b6 100644
--- a/postprocessing/calculate_FWHM.m
+++ b/postprocessing/calculate_FWHM.m
@@ -12,7 +12,7 @@ end
if ~exist('downField_MM_Met', 'var')
- downField_MM_Met = 'DF';
+ downField_MM_Met = 'MM';
end
switch downField_MM_Met
@@ -20,10 +20,21 @@ switch downField_MM_Met
keyMM_ppmNaming = {'DF58' 'DF60' 'DF61' 'DF68' 'DF70' 'DF73' 'DF75' 'DF82' 'DF83' 'DF835' 'DF85' 'NAAB' 'hCs' 'NAA_DF' 'NAD' 'Cr'};%'hist'
orderMet = {'DF58';'DF60';'DF61';'DF68';'DF70';'DF73';'DF75';'DF82';'DF83';'DF835';'DF85';'NAA Broad';'hCs';'NAA'; 'NAD'; 'Cr'};%;'hist'
case 'MM'
+% keyMM_ppmNaming = {'MM09' 'MM12' 'MM14' 'MM17' 'MM20' 'MM22' 'MM26' 'MM27' ...
+% 'MM30' 'MM32' 'MM36' 'MM37' 'MM38' 'Cr39' 'Cr30' 'NAA'};
+% orderMet = { 'M_{0.92}';'M_{1.21}';'M_{1.39}';'M_{1.67}';'M_{2.04}';'M_{2.26}';'M_{2.56}';'M_{2.70}';...
+% 'M_{2.99}';'M_{3.21}';'M_{3.62}';'M_{3.75}';'M_{3.86}'; 'tCr(CH_2)'; 'tCr(CH_3)'; 'NAA'};
+
+% keyMM_ppmNaming = {'MM09' 'MM12' 'MM14' 'MM17' 'MM20' 'MM22' 'MM26' 'MM27' ...
+% 'MM30' 'MM32' 'MM36' 'MM37' 'MM38' 'NAA' 'Creat' 'Cho' 'Glu' 'Tau' 'Glycin' 'PCr'};
+% orderMet = { 'M_{0.92}';'M_{1.21}';'M_{1.39}';'M_{1.67}';'M_{2.04}';'M_{2.26}';'M_{2.56}';'M_{2.70}';...
+% 'M_{2.99}';'M_{3.21}';'M_{3.62}';'M_{3.75}';'M_{3.86}'; 'NAA'; 'tCr(CH_3)'; 'tCho'; 'Glu'; 'Tau'; 'Glyc'; 'tCr(CH_2)'};
+
keyMM_ppmNaming = {'MM09' 'MM12' 'MM14' 'MM17' 'MM20' 'MM22' 'MM26' 'MM27' ...
- 'MM30' 'MM32' 'MM36' 'MM37' 'MM38' 'MM40' 'Cre'};
+ 'MM30' 'MM32' 'MM36' 'MM37' 'MM38' 'MM39' 'Cr39' 'NAA' 'NAA_as' 'Cho' 'mI'};
orderMet = { 'M_{0.92}';'M_{1.21}';'M_{1.39}';'M_{1.67}';'M_{2.04}';'M_{2.26}';'M_{2.56}';'M_{2.70}';...
- 'M_{2.99}';'M_{3.21}';'M_{3.62}';'M_{3.75}';'M_{3.86}';'M_{4.03}'; 'tCr(CH_2)'};
+ 'M_{2.99}';'M_{3.21}';'M_{3.62}';'M_{3.75}';'M_{3.86}'; 'M_{4.03}'; 'tCr(CH_2)'; 'NAA'; 'NAA_{asp}';'tCho'; 'mI'}
+
case 'MM WM'
keyMM_ppmNaming = {'MM09' 'MM12' 'MM14' 'MM17' 'MM20' 'MM22' 'MM26' 'MM27' ...
'MM30' 'MM32' 'MM36' 'MM37' 'MM38' 'MM40' 'Cre'};
@@ -43,15 +54,15 @@ numberOfMet = length(orderMet);
%threshold for crlbs
threshold_CRLB = 900;
FontSize = 14;
-LineWidth = 1.5;
+LineWidth = 2.5;
ppmEnd = 4.2;
ppmStart = 0.2;
interpolationFactor = 20;
%plotting offsets
-offsetMetabolite = 0.1;
+offsetMetabolite = 0.085;
offsetResidual = offsetMetabolite * (numberOfMet+1);
-offsetBaseline = offsetResidual + 0.05;
+offsetBaseline = offsetResidual + 0.042;
%% do the absolute quantification on all files
if ~iscell(FileNames)
@@ -93,7 +104,7 @@ for index =1:length(FileNames)
if strcmp(downField_MM_Met,'DF')
xlim([5.5 9.4]);
else
- xlim([0.5 4.097]);
+ xlim([0.5 4.02]);
end
for plots = 1:length(p)
set(p(plots),'LineWidth',LineWidth);
diff --git a/postprocessing/extractFits.m b/postprocessing/extractFits.m
index d44a9314dc0e1c3e08d82760500d50624dad70eb..72b3ddf32c6436d068f057f75dbf5f59d15ec814 100644
--- a/postprocessing/extractFits.m
+++ b/postprocessing/extractFits.m
@@ -109,8 +109,10 @@ indexOfValue = find(strcmpi(c,'Ph:'))+1; % 'Ph: xx deg xx.x deg/ppm'
phase0 = str2double( c(indexOfValue,1));
phase1 = str2double( c(indexOfValue+2,1));
+
tempX = strfind(c,'ppmgap(2,1)=');
indexOfValue = find(~cellfun(@isempty,tempX));
+
if isempty(indexOfValue)
ppmGap = [];
else
diff --git a/postprocessing/plot_All_Fit_Metabolites.m b/postprocessing/plot_All_Fit_Metabolites.m
index c7c9e7399628138e14ec5db498ac19583867c49a..eea39bf0636841e0d8807047a64f68b8ac47433f 100644
--- a/postprocessing/plot_All_Fit_Metabolites.m
+++ b/postprocessing/plot_All_Fit_Metabolites.m
@@ -8,10 +8,13 @@ if ~iscell(FileName)
FileName= {FileName};
end
+
+
if ~exist('downField_MM_Met','var')
downField_MM_Met = 'Met moiety';
downField_MM_Met = 'DF';
- % downField_MM_Met = 'MM';
+% downField_MM_Met = 'MM';
+ downField_MM_Met = 'fit';
end
overviewPlots = false;
@@ -70,15 +73,25 @@ switch downField_MM_Met
end
case '31P'
- metaboliteNames = {'PCr' 'ATP-g' 'ATP-a' 'ATP-b' 'GPC' 'GPE' 'Pi_in' 'Pi_ext' 'PC' 'PE' 'NADH' 'NAD+' 'UDPG'};
- metaboliteLabels = {'PCr' 'ATP-gamma' 'ATP-alpha' 'ATP-beta' 'GPC' 'GPE' 'Pi intra' 'Pi extra' 'PC' 'PE' 'NADH' 'NAD+' 'UDPG'};
- xLimits = [-20 10];
- offsetMetabolite = 0.1;
- case 'MM'
- metaboliteNames = {'MM09' 'MM12' 'MM14' 'MM17' 'MM20' 'MM22' 'MM26' 'MM27' 'MM30' 'MM32' 'MM36' 'MM37' 'MM38' 'MM40' 'Cre'};
- metaboliteLabels = {'MM09' 'MM12' 'MM14' 'MM17' 'MM20' 'MM22' 'MM26' 'MM27' 'MM30' 'MM32' 'MM36' 'MM37' 'MM38' 'MM40' 'Cre'};
- xLimits = [0.5 4.097];
- offsetMetabolite = 0.1;
+
+ metaboliteNames = {'Cr' 'ATP-g' 'ATP-a' 'GPC' 'GPE' 'Pi_in' 'Pi_ext' 'PC' 'PE' 'NADH' 'NAD' 'UDPG'};
+ metaboliteLabels = {'PCr' 'ATP-gamma' 'ATP-alpha' 'GPC' 'GPE' 'Pi intra' 'Pi extra' 'PC' 'PE' 'NADH' 'NAD+' 'UDPG'};
+ xLimits = [-10 10];
+ case 'fit'
+
+ metaboliteNames = {}; %, 'NAA', 'Cr39', 'MMB_si', 'NAA', 'Cr39', 'Cr30', 'mI' 'combin'
+ metaboliteLabels = {}; %, 'NAA', 'Cr', 'MM' , 'NAA', 'Cr-CH_2', 'Cr-CH_3', 'mI' 'MM'
+ xLimits = [0.5 4.05]
+ offsetMetabolite = 0.15; %0.25
+
+ case 'MM'
+
+ metaboliteNames = {'MM09', 'MM12', 'MM14', 'MM17', 'MM20', 'MM22', 'MM26', 'MM27', 'MM30', 'MM32', 'MM36', 'MM37', 'MM38', 'Cr39', 'NAA', 'mI' ,'GPC', 'Glu', 'NAA_as'};% %, 'Cho', 'Glu', 'Gln', 'Tau', 'Glycin' 'PCr' %, 'NAA_si', 'Creat', 'PCreat', 'GPC','Glu', 'Gln', 'GABA','NAA_as','Asp','Tau', 'mI', 'Glycin'
+ metaboliteLabels = {'M_{0.92}', 'M_{1.21}', 'M_{1.39}', 'M_{1.67}', 'M_{2.04}', 'M_{2.26}', 'M_{2.56}', 'M_{2.70}', 'M_{2.99}', 'M_{3.21}', 'M_{3.62}', 'M_{3.75}', 'M_{3.86}' , 'tCr(CH_2)', 'NAA(CH_3)', 'mI', 'GPC', 'Glu', 'NAA_(asp)'}; %%'NAA', 'tCr(CH_{3})', 'tCho', 'Glu', 'Gln', 'Tau', 'Glyc', %, , 'NAA(CH_3)','tCr(CH_3)', 'tCr(CH_2)','GPC','Glu','Gln','GABA','NAA_{asp}', 'Asp', 'Tau', 'mI', 'Glycin'
+ xLimits = [0.5 4.097]
+ offsetMetabolite = 0.60; %0.25
+
+
otherwise
error('Not known downField_MM_Met');
end
@@ -91,6 +104,7 @@ end
numberOfMet = length(metaboliteNames);
%plotting offsets
+
offsetResidual = offsetMetabolite * (numberOfMet+1);
offsetBaseline = offsetResidual + offsetMetabolite;
@@ -98,6 +112,7 @@ FontSize = 18;
% LineWidth = 1.5;
FontSize = 12;
LineWidth = 2;
+
%
ppm = {};
pData = {};
@@ -216,6 +231,7 @@ for index =1:length(FileName)
metaboliteSpectrum = metaboliteSpectrum + str2double(c(indexOfMet:endOfMet,1)) - bData{index};
end
end
+
if MMB_available
pMetabolite = plot(ppm{index}, metaboliteSpectrum ./ scale - indexMetabolite * offsetMetabolite);
text(xLimitsText, indexMetabolite * -offsetMetabolite,metaboliteLabels{indexMetabolite}, 'FontSize', FontSize, 'FontWeight', 'bold');
@@ -223,6 +239,7 @@ for index =1:length(FileName)
pMetabolite = plot(ppm{index}, metaboliteSpectrum ./ scale - (indexMetabolite-1) * offsetMetabolite);
text(xLimitsText, (indexMetabolite-1) * -offsetMetabolite,metaboliteLabels{indexMetabolite}, 'FontSize', FontSize, 'FontWeight', 'bold');
end
+
set(pMetabolite,'LineWidth',LineWidth);
end
@@ -239,13 +256,16 @@ for index =1:length(FileName)
% ylim([-0.21 0.1]);
% ylim([-0.16 0.115]);
end
+
set(gca,'xDir','reverse')
set(gca,'ytick',[]);
% title(['Summed Metabolite Spectrum with Fitted Metabolites TE = 24 ms'])
% title(FileName{index}, 'Interpreter', 'none')
+
set(gca,'fontsize',FontSize);
set(gca,'FontWeight','bold');
+
clear nOfPoints indexOfMet i
if saveFigures
diff --git a/postprocessing/subtractMetabolitesFromMMB.m b/postprocessing/subtractMetabolitesFromMMB.m
index 9f3a0e6b31266049a98afd760e49db286c80aa74..50a3bffb658bea89c79b44286b5b219b11ceca63 100644
--- a/postprocessing/subtractMetabolitesFromMMB.m
+++ b/postprocessing/subtractMetabolitesFromMMB.m
@@ -1,16 +1,16 @@
function subtractMetabolitesFromMMB()
-pathName = 'MM data path';
-sampleCriteria = {'Output'};
+pathName = 'MM T1 data path';
+sampleCriteria = {'Output_GM_FinalCorrected'};
[localFilePathBase] = pathToDataFolder(pathName, sampleCriteria);
-LCModelOutputPath = [localFilePathBase, 'Output/'];
+LCModelOutputPath = [localFilePathBase, 'Output_GM_FinalCorrected/'];
-outputFileNameBase = 'Summed_Averaged_TE';
-fileNameBaseWithoutMetabolite = 'MMB_without_metabolite_TE';
-orderTEs = {'24' '32' '40' '52' '60'};
+outputFileNameBase = 'Summed_Averaged_TI1_';
+fileNameBaseWithoutMetabolite = 'MMB_without_metabolite_TI1_';
+orderTEs = {'1900'}; %'2360' '2000' '1900' '1800'
extensionCoord = '.coord';
extensionMat = '.mat';
-metabolitesToSubtract = {'Cre'};
+metabolitesToSubtract = {'Cr39', 'Cr30', 'NAA', 'Glx'};
doNullingOfDownField = true;
LCModelTableFitsCoord = strcat(outputFileNameBase, orderTEs, extensionCoord);
reconstructedSpectra = strcat(localFilePathBase, outputFileNameBase, orderTEs, extensionMat);
@@ -28,25 +28,25 @@ LineWidth = 1.5;
for indexTE = 1:numberOfTEs
currentReconstructedSpectra = reconstructedSpectra{indexTE};
- load(currentReconstructedSpectra, 'summedSpectraTE');
+ load(currentReconstructedSpectra, 'summedSpectraTI');
currentLCModelCoordFit = LCModelTableFitsCoord{indexTE};
[nOfPoints, ppmVector, phasedData, fitData, baselineData, residualData, ...
metaboliteNames, metaboliteData, tableConcentrations, ppmStart, ppmEnd, ...
scanFrequency, frequencyShift, ppmGap] = ...
extractFits(LCModelOutputPath,currentLCModelCoordFit);
-
+ ppmGap = 4.1
% align the two spectra
- fid = summedSpectraTE.Data{1};
+ fid = summedSpectraTI.Data{1};
spectrum = fftshift(fft(fid(:,1,1,1,1,1,1,1,1,1,1,1),size(fid,1)*zeroFillFactor));
- ppmSpectrum = summedSpectraTE.getPpmVector(zeroFillFactor);
+ ppmSpectrum = summedSpectraTI.getPpmVector(zeroFillFactor);
[~, maxPpmValueSpectrum] = getMaxPeak(ppmSpectrum', real(spectrum), 3.925, 0.1);
[~, maxPpmValueLCModel] = getMaxPeak(ppmVector, phasedData, 3.925, 0.1);
frequencyShiftDiff = maxPpmValueSpectrum - maxPpmValueLCModel;
- summedSpectraTE = summedSpectraTE.ApplyFrequencyShift(frequencyShiftDiff);
+ summedSpectraTI = summedSpectraTI.ApplyFrequencyShift(frequencyShiftDiff);
%get the aligned reconstructed spectra
- fid = summedSpectraTE.Data{1};
+ fid = summedSpectraTI.Data{1};
spectrum = fftshift(fft(fid(:,1,1,1,1,1,1,1,1,1,1,1),size(fid,1)*zeroFillFactor));
% max of the reconstructed spectrum
if isempty(ppmGap)
@@ -122,7 +122,7 @@ for indexTE = 1:numberOfTEs
% xlim([0 4.2])
% title(sprintf('Imag MMB TE %s', orderTEs{indexTE}))
- MMB_without_metabolite = summedSpectraTE;
+ MMB_without_metabolite = summedSpectraTI;
newFid = ifft(ifftshift(newSpectrum));
MMB_without_metabolite.Data{1} = newFid;
diff --git a/reconstruct_all_Met.m b/reconstruct_all_Met.m
new file mode 100644
index 0000000000000000000000000000000000000000..11fb0de70cd4fee488a9baa0376a0ee850b1feb1
--- /dev/null
+++ b/reconstruct_all_Met.m
@@ -0,0 +1,161 @@
+function reconstruct_all_Met()
+
+pathName = 'Met T1 GM data path';
+subjects = { ...
+ '1405' ...
+ '1706' ...
+ '1717' ...
+ '2016' ...
+ '2017' ...
+ '3373' ...
+ '3490' ...
+ '6971' ...
+ '7338' ...
+ '7782' ...
+ '9810' ...
+ };
+
+pathBase = pathToDataFolder(pathName, subjects);
+
+TI = [20; 100; 400; 700; 1000; 1400; 2500];
+
+numberOfSubjects = length(subjects);
+paths = cell(1,numberOfSubjects);
+for indexSubj = 1:numberOfSubjects
+ paths{indexSubj} = [pathBase subjects{indexSubj} '\\'];
+end
+
+files_TI_20 = cell(1, numberOfSubjects);
+files_TI_100 = cell(1, numberOfSubjects);
+files_TI_400 = cell(1, numberOfSubjects);
+files_TI_700 = cell(1, numberOfSubjects);
+files_TI_1000 = cell(1, numberOfSubjects);
+files_TI_1400 = cell(1, numberOfSubjects);
+files_TI_2500 = cell(1, numberOfSubjects);
+
+% water_files = cell(1, numberOfSubjects);
+for indexSubject = 1:numberOfSubjects
+ files_TI_20{indexSubject} = ls([paths{indexSubject}, '*_semiLASER_2360_625_*.dat']);
+ files_TI_100{indexSubject} = ls([paths{indexSubject}, '*_semiLASER_1800_525_*.dat']);
+ files_TI_400{indexSubject} = ls([paths{indexSubject}, '*_semiLASER_1900_550_*.dat']);
+ files_TI_700{indexSubject} = ls([paths{indexSubject}, '*_semiLASER_2000_575_*.dat']);
+ files_TI_1000{indexSubject} = ls([paths{indexSubject}, '*_semiLASER_2150_600_*.dat']);
+ files_TI_1400{indexSubject} = ls([paths{indexSubject}, '*_semiLASER_1200_20_*.dat']);
+ files_TI_2500{indexSubject} = ls([paths{indexSubject}, '*_semiLASER_1250_20_*.dat']);
+% files_TI_1300_80{indexSubject} = ls([paths{indexSubject}, '*_semiLASER_1300_80_*.dat']);
+% files_TI_1300_60{indexSubject} = ls([paths{indexSubject}, '*_semiLASER_1300_60_*.dat']);
+% files_TI_1300_20{indexSubject} = ls([paths{indexSubject}, '*_semiLASER_1300_20_*.dat']);
+% files_TI_1050_238{indexSubject} = ls([paths{indexSubject}, '*_semiLASER_1050_238_*.dat']);
+% water_files{indexSubject} = ls([paths{indexSubject}, '*_scout_wref_*.dat']);
+end
+
+addSinglet0ppm = 'Yes';
+numberOfTIs = length(TI);
+
+filesPath = cell(numberOfTIs,numberOfSubjects);
+waterFilesPath = cell(1,numberOfSubjects);
+files = [files_TI_20; files_TI_100; files_TI_400;files_TI_700;...
+ files_TI_1000; files_TI_1400; files_TI_2500;];
+for indexSubj = 1:numberOfSubjects
+ for indexTI = 1:numberOfTIs
+ filesPath{indexTI, indexSubj} = [paths{indexSubj} files{indexTI, indexSubj}];
+ end
+% waterFilesPath{indexSubj} = [paths{indexSubj} water_files{indexSubj}];
+end
+
+summedSpectra = cell(numberOfTIs,1);
+indexSubjTIs = zeros(numberOfTIs,1);
+for indexSubj = 1:numberOfSubjects
+ %% water data
+% [a] = reconstructWater(waterFilesPath{indexSubj});
+% filename = [paths{indexSubj} subjects{indexSubj} '_water.mat'];
+% save(filename, 'a');
+% a.ExportLcmRaw(paths{indexSubj}, strcat(subjects{indexSubj}, '_water'), addSinglet0ppm);
+ %calculate scaling factor the summedSpectra individually according to water peak
+% waterData = a.Data{1};
+% maxWaterPeak = max(real(fftshift(fft(waterData))));
+% %scale the data for equalizing between subjects
+% a.Data{1} = a.Data{1}./ maxWaterPeak;
+ %create summed spectra
+% if(indexSubj == 1)
+% summedWater = a;
+% newDataWater= zeros(length(a.Data{1}),1,1,1,1,1,1,1,1,1,1,numberOfSubjects);
+% else
+% oldDataWater = summedWater.Data{1};
+% if (ndims(oldDataWater) ~= 12)
+% newDataWater(:,1,1,1,1,1,1,1,1,1,1,1) = oldDataWater;
+% else
+% newDataWater(:,1,1,1,1,1,1,1,1,1,1,:) = oldDataWater;
+% end
+% newDataWater(:,1,1,1,1,1,1,1,1,1,1,indexSubj) = a.Data{1};
+% summedWater.Data{1} = newDataWater;
+% end
+
+
+
+ %% reconstruct Macromolecule data
+ for indexTI = 1:numberOfTIs
+ if((TI(indexTI) == 1300)| (TI(indexTI) ==1250) | (TI(indexTI)==1200))
+ doReverse = false;
+ else
+ doReverse = false;
+ end
+ a = reconstruct(filesPath{indexTI, indexSubj}, 1, 1);
+ filename = [paths{indexSubj} subjects{indexSubj} '_TI_' num2str(TI(indexTI)) '.mat'];
+ save(filename, 'a');
+ a.ExportLcmRaw(paths{indexSubj}, strcat(subjects{indexSubj}, '_TI1_', num2str(TI(indexTI))) , addSinglet0ppm);
+ %scale the data for equalizing between subjects
+% a.Data{1} = a.Data{1}./ maxWaterPeak;
+
+ % eliminate data with too many deleted averages from the summedSpectra
+
+ deletedAverages = a.Parameter.DeleteMovedAveragesSettings.deletedAverages;
+ if ~isempty(deletedAverages)
+ numDeletedAverages = length(deletedAverages);
+ else
+ numDeletedAverages = 0;
+ end
+
+ if numDeletedAverages <= 2
+ indexSubjTIs(indexTI) = indexSubjTIs(indexTI) + 1;
+ %create summed spectra
+ if(indexSubjTIs(indexTI) == 1)
+ summedSpectra{indexTI} = a;
+ newData= zeros(length(a.Data{1}),1,1,1,1,1,1,1,2,1,1,numberOfSubjects);
+ else
+ oldData = summedSpectra{indexTI}.Data{1};
+ if (ndims(oldData) ~= 12)
+ newData(:,1,1,1,1,1,1,1,:,1,1,1) = oldData;
+ else
+ newData(:,1,1,1,1,1,1,1,:,1,1,:) = oldData;
+ end
+ newData(:,1,1,1,1,1,1,1,:,1,1,indexSubjTIs(indexTI)) = a.Data{1};
+ summedSpectra{indexTI}.Data{1} = newData;
+ end
+ end
+ end
+end
+
+%
+
+% filename = [pathBase, 'Summed_water.mat'];
+% save(filename, 'summedWater');
+% summedWater = summedWater.AverageData;
+% save([pathBase, 'Summed_Averaged_water.mat'], 'summedWater');
+% summedWater.ExportLcmRaw(pathBase, 'Summed_Averaged_water', addSinglet0ppm);
+for indexTI = 1:numberOfTIs
+ filename = [pathBase, 'Summed_spectra_TI' num2str(TI(indexTI)) '.mat'];
+ summedSpectraTI = summedSpectra{indexTI};
+ %store only the data, which was kept after eliminating the deleted ones
+ actualData = summedSpectraTI.Data{1}./indexSubjTIs(indexTI);
+ summedSpectraTI.Data{1} = actualData(:,1,1,1,1,1,1,1,:,1,1,1:indexSubjTIs(indexTI));
+ save(filename, 'summedSpectraTI');
+ %average across subjects
+ summedSpectraTI = summedSpectraTI.AverageData;
+ filename = [pathBase, 'Summed_Averaged_TI_' num2str(TI(indexTI)) '.mat'];
+ save(filename, 'summedSpectraTI');
+ summedSpectraTI.ExportLcmRaw(pathBase, strcat('Summed_Averaged_TI_', num2str(TI(indexTI))), addSinglet0ppm);
+end
+
+%summarize_Deleted_Averages(subjects, pathBase, TI1)
+close all