Skip to content
Snippets Groups Projects
Commit c839aa68 authored by Tamas Borbath's avatar Tamas Borbath
Browse files

Updated createSimulatedSpectra and createSummedSpectrum to include stats about... 

Updated createSimulatedSpectra and createSummedSpectrum to include stats about SNR etc. Changed MM spectrum to have 0 em (Lorentz) model
parent 291f45c8
Branches
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
createTestData/createSimulatedSpectra.m
+ 192
− 89
View file @ c839aa68
......@@ -15,7 +15,7 @@ TE = 24 * 1e-3;
numMetabolites = length(metabolites);
saveData = false;
saveFigures = true;
saveFigures = false;
saveRandNumbers = false; %should generally be kept false to have the same random numbers consistently
%% path to the export of the spectra
......@@ -41,9 +41,7 @@ df2_local_std = 3;% 3 Hz coresponds to ~0.0075 ppm.
% pH changes should have a similar order of magnitude.
em_local_scale = [0.2 0.2 0.4 0.4 0.6 0.6 0.8 0.8 1 1 1.2 1.2 1.5 1.5 2];
% when truncating at 200 ms like used in real spectra, the NAA SNR of the defaultSummedFid corresponds to
% Chosen NoiseLevel -> SNR(NAA)
% 1 -> SNR = 207; 3 -> SNR = 158; 7 -> SNR = 110; 15 -> SNR = 44; 25 -> SNR = 13
%check the resulting SNR levels from noiseLevel in snr_NAA_noise in the noise simulation
noiseLevel = [1:25];
numberDfPermutations = length(df2_global);
......@@ -113,6 +111,7 @@ phase_1_f2_mx = ppmVector-phasePivot_ppm;
ppmMask = (ppmVector>0.6) & (ppmVector<4.1);
legendLabels = cell(3,1);
%% set default values for the parameters
default_pc0 = 0;
......@@ -134,23 +133,26 @@ current_em_std_factor = default_em_std_factor;
conc_std_factor = zeros(1,numMetabolites);
% generate the default spectrum
[defaultSummedFid, defaultSummedFidWithNoise, defaultCurrent_em, defaultCurrentConcentrationRm] = ...
[defaultSummedFid, defaultSummedFidWithNoise, defaultCurrent_em, defaultCurrentConcentrationRm, defaultStats] = ...
createSummedSpectrum(...
fids, concentrations, conc_std, T2s, T2_std, TE, metabolites, metaboliteFileNames, timePointsSampling, phase_1_f2_mx, ppmMask, ...
current_pc0, current_pc1, current_df2, current_gm, current_em_std_factor, conc_std_factor, noiseVector, ...
truncPoint, saveData, 'meanConc', [], dwellTime * 1e3, scanFrequency, dataExportPathBase, water_ppm, scalingFactorRefPeak);
plot(ppmVector,real(fftshift(fft(defaultSummedFid))));
% figure(1);
% subplot(5,4,8+7:8+8)
figure(1);
subplot(5,4,8+7:8+8)
subplot(3,4,9:10)%Revision 1
% figure(2);
% subplot(3,4,7:8)
plot(ppmVector,real(fftshift(fft(defaultSummedFid))));
%% generate the series of different concentrations
stdLimConc = 3.5;
numberConcPermutations = 100;
concentrationsRm = zeros(numberConcPermutations, numMetabolites);
paramKeyword = 'conc';
titleStr = '{\boldmath$c$}';
if saveRandNumbers
%generally you should not enter this if. It creates random numbers which can be used for calculating the concentrations
concentrationRandNumbers = rand(numberConcPermutations,numMetabolites);
......@@ -158,6 +160,7 @@ if saveRandNumbers
end
load([dataExportPathBase, 'concentrationRandNumbers.mat'], 'concentrationRandNumbers');
dataExportPath = [dataExportPathBase, paramKeyword, '\'];
snr_NAA_conc = zeros(numberConcPermutations,1);
for indexParam = 1:numberConcPermutations
for indexMetabolite = 1:numMetabolites
% generate a random number between [-stdLim; stdLim]
......@@ -176,28 +179,28 @@ for indexParam = 1:numberConcPermutations
end
conc_std_factor(indexMetabolite) = stdRand;
end
[summedFid, summedFidWithNoise, current_em, currentConcentrationRm] = createSummedSpectrum(...
[summedFid, summedFidWithNoise, current_em, currentConcentrationRm, stats] = createSummedSpectrum(...
fids, concentrations, conc_std, T2s, T2_std, TE, metabolites, metaboliteFileNames, timePointsSampling, phase_1_f2_mx, ppmMask, ...
current_pc0, current_pc1, current_df2, current_gm, current_em_std_factor, conc_std_factor, noiseVector, ...
truncPoint, saveData, paramKeyword, indexParam, dwellTime * 1e3, scanFrequency, dataExportPath, referencePeakPpm, scalingFactorRefPeak);
snr_NAA_conc(indexParam) = stats{2,2};
concentrationsRm(indexParam,:) = currentConcentrationRm;
hold on;
plot(ppmVector,real(fftshift(fft(summedFid))));
plot(ppmVector,real(fftshift(fft(summedFidWithNoise))));
% ImportLCModelBasis(dataExportPath, fileName, true, '1H');
end
titleStr = '{\boldmath$c$}';
plotSimulationSetup(titleStr, paramKeyword, saveFigures, dataExportPathBase);
if saveData
save([dataExportPathBase, 'concentrations_simulated_', paramKeyword, '.mat'], 'concentrationsRm', 'metabolites', 'metaboliteFileNames');
end
for indexMetabolite = 1:numMetabolites
figure;
scatter(concentrationsRm(:,indexMetabolite), concentrationsRm(:,indexMetabolite));
title(metabolites(indexMetabolite))
end
% for indexMetabolite = 1:numMetabolites
% figure;
% scatter(concentrationsRm(:,indexMetabolite), concentrationsRm(:,indexMetabolite));
% title(metabolites(indexMetabolite))
% end
%% #############################################################################################
% generate the series of different parameters, same concentrations
......@@ -206,58 +209,69 @@ conc_std_factor = zeros(1,numMetabolites); % reset value
%% gauss
paramKeyword = 'gauss';
titleStr = '$$\nu_{g}$$';
units = '(Hz)';
dataExportPath = [dataExportPathBase, paramKeyword, '\'];
concentrationsRm = zeros(length(gaussianValues), numMetabolites);
figure(1);
subplot(5,4,5:6)
% figure(1);
% subplot(5,4,5:6)
figure(2);
subplot(2,4,1:2)%Revision 1
% subplot(2,4,5:6)
snr_NAA_Gauss = zeros(length(gaussianValues),1);
for indexParam = 1:length(gaussianValues)
current_gm = gaussianValues(indexParam);
[summedFid, summedFidWithNoise, current_em, currentConcentrationRm] = createSummedSpectrum(...
[summedFid, summedFidWithNoise, current_em, currentConcentrationRm, stats] = createSummedSpectrum(...
fids, concentrations, conc_std, T2s, T2_std, TE, metabolites, metaboliteFileNames, timePointsSampling, phase_1_f2_mx, ppmMask, ...
current_pc0, current_pc1, current_df2, current_gm, current_em_std_factor, conc_std_factor, noiseVector, ...
truncPoint, saveData, paramKeyword, indexParam, dwellTime * 1e3, scanFrequency, dataExportPath, referencePeakPpm, scalingFactorRefPeak);
snr_NAA_Gauss(indexParam) = stats{2,2};
concentrationsRm(indexParam,:) = currentConcentrationRm;
hold on;
plot(ppmVector,real(fftshift(fft(summedFid))));
legendLabels = plotSpectra(summedFidWithNoise, ppmVector, false, legendLabels, indexParam, gaussianValues, '$\rm\nu_{g}$', units);
end
% reset parameter
current_gm = default_gm;
titleStr = '$$\nu_{g}$$';
plotSimulationSetup(titleStr, paramKeyword, saveFigures, dataExportPathBase);
legend(legendLabels, 'Interpreter', 'latex');
plotArrowsPlotOffsets()
%%
if saveData
save([dataExportPathBase, 'concentrations_simulated_', paramKeyword, '.mat'], 'concentrationsRm', 'metabolites', 'metaboliteFileNames');
end
%% pc0 - zero order phase
paramKeyword = 'pc0';
titleStr = '$$\varphi_{0}$$';
units = '(degrees)';
dataExportPath = [dataExportPathBase, paramKeyword, '\'];
concentrationsRm = zeros(length(pc0), numMetabolites);
figure(1);
subplot(5,4,1:2)
% subplot(5,4,1:2)
subplot(3,4,1:2) %Revision 1
% subplot(2,4,1:2)
snr_NAA_pc0 = zeros(length(pc0),1);
for indexParam = 1:length(pc0)
current_pc0 = pc0(indexParam);
[summedFid, summedFidWithNoise, current_em, currentConcentrationRm] = createSummedSpectrum(...
[summedFid, summedFidWithNoise, current_em, currentConcentrationRm, stats] = createSummedSpectrum(...
fids, concentrations, conc_std, T2s, T2_std, TE, metabolites, metaboliteFileNames, timePointsSampling, phase_1_f2_mx, ppmMask, ...
current_pc0, current_pc1, current_df2, current_gm, current_em_std_factor, conc_std_factor, noiseVector, ...
truncPoint, saveData, paramKeyword, indexParam, dwellTime * 1e3, scanFrequency, dataExportPath, referencePeakPpm, scalingFactorRefPeak);
snr_NAA_pc0(indexParam) = stats{2,2};
concentrationsRm(indexParam,:) = currentConcentrationRm;
hold on;
plot(ppmVector,real(fftshift(fft(summedFid))));
legendLabels = plotSpectra(summedFidWithNoise, ppmVector, true, legendLabels, indexParam, pc0, titleStr, units);
end
% reset parameter
current_pc0 = default_pc0;
titleStr = '$$\varphi_{0}$$';
plotSimulationSetup(titleStr, paramKeyword, saveFigures, dataExportPathBase);
if saveData
save([dataExportPathBase, 'concentrations_simulated_', paramKeyword, '.mat'], 'concentrationsRm', 'metabolites', 'metaboliteFileNames');
......@@ -265,28 +279,32 @@ end
%% pc1 - first order phase
paramKeyword = 'pc1';
titleStr = '$$\varphi_{1}$$';
units = '(degrees/ppm)';
dataExportPath = [dataExportPathBase, paramKeyword, '\'];
concentrationsRm = zeros(length(pc1), numMetabolites);
figure(1);
subplot(5,4,3:4)
% subplot(5,4,3:4)
subplot(3,4,3:4)%Revision 1
% subplot(2,4,3:4)
snr_NAA_pc1 = zeros(length(pc1),1);
for indexParam = 1:length(pc1)
current_pc1 = pc1(indexParam);
[summedFid, summedFidWithNoise, current_em, currentConcentrationRm] = createSummedSpectrum(...
[summedFid, summedFidWithNoise, current_em, currentConcentrationRm, stats] = createSummedSpectrum(...
fids, concentrations, conc_std, T2s, T2_std, TE, metabolites, metaboliteFileNames, timePointsSampling, phase_1_f2_mx, ppmMask, ...
current_pc0, current_pc1, current_df2, current_gm, current_em_std_factor, conc_std_factor, noiseVector, ...
truncPoint, saveData, paramKeyword, indexParam, dwellTime * 1e3, scanFrequency, dataExportPath, referencePeakPpm, scalingFactorRefPeak);
snr_NAA_pc1(indexParam) = stats{2,2};
concentrationsRm(indexParam,:) = currentConcentrationRm;
hold on;
plot(ppmVector,real(fftshift(fft(summedFid))));
legendLabels = plotSpectra(summedFidWithNoise, ppmVector, true, legendLabels, indexParam, pc1, titleStr, units);
end
% reset parameter
current_pc1 = default_pc1;
titleStr = '$$\varphi_{1}$$';
plotSimulationSetup(titleStr, paramKeyword, saveFigures, dataExportPathBase);
if saveData
save([dataExportPathBase, 'concentrations_simulated_', paramKeyword, '.mat'], 'concentrationsRm', 'metabolites', 'metaboliteFileNames');
......@@ -294,13 +312,18 @@ end
%% df2 local - frequency shifts
paramKeyword = 'df2_local';
titleStr = '{\boldmath$\omega_{local}$}';
units = '(Hz)';
dataExportPath = [dataExportPathBase, paramKeyword, '\'];
concentrationsRm = zeros(length(df2_local_scale), numMetabolites);
% figure(1);
% subplot(5,4,8+1:8+2)
figure(1);
subplot(5,4,8+1:8+2)
subplot(3,4,7:8)%Revision 1
% figure(2);
% subplot(3,4,1:2)
snr_NAA_df2_local = zeros(length(df2_local_scale),1);
for indexParam = 1:length(df2_local_scale)
for indexMetabolite = 1:numMetabolites
......@@ -320,19 +343,19 @@ for indexParam = 1:length(df2_local_scale)
end
current_df2(indexMetabolite) = df2_local_scale(indexParam) * stdRand;
end
[summedFid, summedFidWithNoise, current_em, currentConcentrationRm] = createSummedSpectrum(...
[summedFid, summedFidWithNoise, current_em, currentConcentrationRm, stats] = createSummedSpectrum(...
fids, concentrations, conc_std, T2s, T2_std, TE, metabolites, metaboliteFileNames, timePointsSampling, phase_1_f2_mx, ppmMask, ...
current_pc0, current_pc1, current_df2, current_gm, current_em_std_factor, conc_std_factor, noiseVector, ...
truncPoint, saveData, paramKeyword, indexParam, dwellTime * 1e3, scanFrequency, dataExportPath, water_ppm, scalingFactorRefPeak);
snr_NAA_df2_local(indexParam) = stats{2,2};
concentrationsRm(indexParam,:) = currentConcentrationRm;
hold on;
plot(ppmVector,real(fftshift(fft(summedFid))));
legendLabels = plotSpectra(summedFidWithNoise, ppmVector, true, legendLabels, indexParam, df2_local_scale, titleStr, units);
end
% reset parameter
current_df2 = default_df2;
titleStr = '{\boldmath$\omega_{local}$}';
plotSimulationSetup(titleStr, paramKeyword, saveFigures, dataExportPathBase);
if saveData
save([dataExportPathBase, 'concentrations_simulated_', paramKeyword, '.mat'], 'concentrationsRm', 'metabolites', 'metaboliteFileNames');
......@@ -340,29 +363,34 @@ end
%% df2 global - frequency shifts
paramKeyword = 'df2_global';
titleStr = '$$\omega_{global}$$';
units = ('Hz');
dataExportPath = [dataExportPathBase, paramKeyword, '\'];
concentrationsRm = zeros(length(df2_global), numMetabolites);
% figure(1);
% subplot(5,4,7:8)
figure(1);
subplot(5,4,7:8)
subplot(3,4,5:6)%Revision 1
% subplot(2,4,7:8)
snr_NAA_df2_global = zeros(length(df2_global),1);
for indexParam = 1:length(df2_global)
for indexMetabolite = 1:numMetabolites
current_df2(indexMetabolite) = df2_global(indexParam);
end
[summedFid, summedFidWithNoise, current_em, currentConcentrationRm] = createSummedSpectrum(...
[summedFid, summedFidWithNoise, current_em, currentConcentrationRm, stats] = createSummedSpectrum(...
fids, concentrations, conc_std, T2s, T2_std, TE, metabolites, metaboliteFileNames, timePointsSampling, phase_1_f2_mx, ppmMask, ...
current_pc0, current_pc1, current_df2, current_gm, current_em_std_factor, conc_std_factor, noiseVector, ...
truncPoint, saveData, paramKeyword, indexParam, dwellTime * 1e3, scanFrequency, dataExportPath, water_ppm, scalingFactorRefPeak);
snr_NAA_df2_global(indexParam) = stats{2,2};
concentrationsRm(indexParam,:) = currentConcentrationRm;
hold on;
plot(ppmVector,real(fftshift(fft(summedFid))));
legendLabels = plotSpectra(summedFidWithNoise, ppmVector, true, legendLabels, indexParam, df2_global, titleStr, units);
end
% reset parameter
current_df2 = default_df2;
titleStr = '$$\omega_{global}$$';
plotSimulationSetup(titleStr, paramKeyword, saveFigures, dataExportPathBase);
if saveData
save([dataExportPathBase, 'concentrations_simulated_', paramKeyword, '.mat'], 'concentrationsRm', 'metabolites', 'metaboliteFileNames');
......@@ -370,13 +398,18 @@ end
%% em - lorentzian components of T2
paramKeyword = 'em';
titleStr = '{\boldmath$\nu_{e}$}';
units = '(Hz)';
dataExportPath = [dataExportPathBase, paramKeyword, '\'];
concentrationsRm = zeros(length(df2_global), numMetabolites);
figure(1);
subplot(5,4,8+3:8+4)
% figure(1);
% subplot(5,4,8+3:8+4)
figure(2);
subplot(2,4,3:4)%Revision 1
% figure(2);
% subplot(3,4,3:4)
snr_NAA_em = zeros(length(df2_global),1);
for indexParam = 1:length(df2_global)
for indexMetabolite = 1:numMetabolites
% generate a random number between [-stdLim; stdLim]
......@@ -395,20 +428,22 @@ for indexParam = 1:length(df2_global)
end
current_em_std_factor(indexMetabolite) = em_local_scale(indexParam) * stdRand;
end
[summedFid, summedFidWithNoise, current_em, currentConcentrationRm] = createSummedSpectrum(...
[summedFid, summedFidWithNoise, current_em, currentConcentrationRm, stats] = createSummedSpectrum(...
fids, concentrations, conc_std, T2s, T2_std, TE, metabolites, metaboliteFileNames, timePointsSampling, phase_1_f2_mx, ppmMask, ...
current_pc0, current_pc1, current_df2, current_gm, current_em_std_factor, conc_std_factor, noiseVector, ...
truncPoint, saveData, paramKeyword, indexParam, dwellTime * 1e3, scanFrequency, dataExportPath, referencePeakPpm, scalingFactorRefPeak);
snr_NAA_em(indexParam) = stats{2,2};
concentrationsRm(indexParam,:) = currentConcentrationRm;
hold on;
plot(ppmVector,real(fftshift(fft(summedFid))));
legendLabels = plotSpectra(summedFidWithNoise, ppmVector, false, legendLabels, indexParam, em_local_scale, '$\rm T_{2,k} = T_{2,k} + X_{-1,1}^{k} \cdot $', '$\rm\cdot std(T_{2,k})$');
end
% reset parameter
current_em_std_factor = default_em_std_factor;
titleStr = '{\boldmath$\nu_{e}$}';
plotSimulationSetup(titleStr, paramKeyword, saveFigures, dataExportPathBase);
legend(legendLabels, 'Interpreter', 'latex');
plotArrowsPlotOffsets()
if saveData
save([dataExportPathBase, 'concentrations_simulated_', paramKeyword, '.mat'], 'concentrationsRm', 'metabolites', 'metaboliteFileNames');
end
......@@ -416,35 +451,44 @@ end
%% noise - add different noise levels
paramKeyword = 'noise';
titleStr = '{\boldmath$noise$}';
units = '';
dataExportPath = [dataExportPathBase, paramKeyword, '\'];
concentrationsRm = zeros(length(noiseLevel), numMetabolites);
figure(1);
subplot(5,4,8+5:8+6)
% figure(1);
% subplot(5,4,8+5:8+6)
figure(2);
subplot(2,4,6:7)%Revision 1
% figure(2);
% subplot(3,4,5:6)
snr_NAA_noise = zeros(length(noiseLevel),1);
for indexParam = 1:length(noiseLevel)
noiseVector = wgn(1,nTimepoints, noiseLevel(indexParam));
[summedFid, summedFidWithNoise, current_em, currentConcentrationRm] = createSummedSpectrum(...
[summedFid, summedFidWithNoise, current_em, currentConcentrationRm, stats] = createSummedSpectrum(...
fids, concentrations, conc_std, T2s, T2_std, TE, metabolites, metaboliteFileNames, timePointsSampling, phase_1_f2_mx, ppmMask, ...
current_pc0, current_pc1, current_df2, current_gm, current_em_std_factor, conc_std_factor, noiseVector, ...
truncPoint, saveData, paramKeyword, indexParam, dwellTime * 1e3, scanFrequency, dataExportPath, referencePeakPpm, scalingFactorRefPeak);
snr_NAA_noise(indexParam) = stats{2,2};
concentrationsRm(indexParam,:) = currentConcentrationRm;
hold on;
plot(ppmVector,real(fftshift(fft(summedFidWithNoise))));
legendLabels = plotSpectra(summedFidWithNoise, ppmVector, false, legendLabels, indexParam, snr_NAA_noise, '$\rm SNR_{NAA}$', units);
end
% reset parameter
noiseVector = wgn(1,nTimepoints, default_noiseLevel);
titleStr = '{\boldmath$noise$}';
plotSimulationSetup(titleStr, paramKeyword, saveFigures, dataExportPathBase);
legend(legendLabels, 'Interpreter', 'latex');
plotArrowsPlotOffsets()
if saveData
save([dataExportPathBase, 'concentrations_simulated_', paramKeyword, '.mat'], 'concentrationsRm', 'metabolites', 'metaboliteFileNames');
end
%% noise - add different noise levels
paramKeyword = 'baseline';
titleStr = '{\boldmath$baseline$}';
units = '';
dataExportPath = [dataExportPathBase, paramKeyword, '\'];
concentrationsRm = zeros(length(noiseLevel), numMetabolites);
......@@ -453,30 +497,63 @@ prefixBaseline = 'Baseline_';
baselines = {'1Flat_drop_H20'; '2Flat'; '3Flat'; '4Flat'; '5Wildish'; '6Wildish'; '7Very_wild';...
'8Very_wild'; '9Wild'; '10Wild'; '11Wildish'; '12Wildish'; '13Normal'; '14Lipid_1.3'; '15Lipid_1.3_phased'; '16Zero'};
% figure(1);
% subplot(5,4,8+10:8+11)
figure(1);
subplot(5,4,8+10:8+11)
subplot(3,4,11:12)
% figure(2);
% subplot(3,4,10:11)
snr_NAA_baselines = zeros(length(baselines),1);
for indexParam = 1:length(baselines)
baselineFileName = strcat(prefixBaseline, baselines{indexParam});
[currentBaseline, ~, ~, ~, ~] = ImportLCModelBasis(baselinesPath, baselineFileName, plotBasis, '1H');
currentBaseline = currentBaseline(1:length(fids{1}));% will fail if the fids are longer than the baseline
[summedFid, summedFidWithNoise, current_em, currentConcentrationRm] = createSummedSpectrum(...
[summedFid, summedFidWithNoise, current_em, currentConcentrationRm, stats] = createSummedSpectrum(...
fids, concentrations, conc_std, T2s, T2_std, TE, metabolites, metaboliteFileNames, timePointsSampling, phase_1_f2_mx, ppmMask, ...
current_pc0, current_pc1, current_df2, current_gm, current_em_std_factor, conc_std_factor, noiseVector, ...
truncPoint, saveData, paramKeyword, indexParam, dwellTime * 1e3, scanFrequency, dataExportPath, referencePeakPpm, scalingFactorRefPeak, currentBaseline);
snr_NAA_baselines(indexParam) = stats{2,2};
concentrationsRm(indexParam,:) = currentConcentrationRm;
hold on;
plot(ppmVector,real(fftshift(fft(summedFidWithNoise))));
legendLabels = plotSpectra(summedFidWithNoise, ppmVector, true, legendLabels, indexParam, baselines, titleStr, units);
end
titleStr = '{\boldmath$baseline$}';
plotSimulationSetup(titleStr, paramKeyword, saveFigures, dataExportPathBase);
if saveData
save([dataExportPathBase, 'concentrations_simulated_', paramKeyword, '.mat'], 'concentrationsRm', 'metabolites', 'metaboliteFileNames');
end
end
function legendLabels = plotSpectra(summedFidWithNoise, ppmVector, plotAll, legendLabels, indexParam, ...
paramValues, titleStr, units)
if (plotAll == true)
%plot each fid
hold on;
plot(ppmVector,real(fftshift(fft(summedFidWithNoise))));
else
if ~contains(titleStr,'=')
titleStr = [titleStr, ' = '];
end
% plot just the max, min and median value
offsetPlot = -1;
if (indexParam == 1)
offsetPlot = 0;
end
if (indexParam == round(length(paramValues)/2))
offsetPlot = 1;
end
if (indexParam == length(paramValues))
offsetPlot = 2;
end
if (offsetPlot ~= -1)
hold on;
plot(ppmVector,real(fftshift(fft(summedFidWithNoise)))-offsetPlot*2e3);
legendLabels{offsetPlot+1} = [titleStr, num2str(paramValues(indexParam)), ' ', units];
end
end
end
function plotSimulationSetup(titleStr, paramKeyword, saveFigures, dataExportPathBase)
FontSize = 12;
LineWidth = 1.5;
......@@ -498,3 +575,29 @@ if saveFigures
savefig(gcf, [fileNameFull, '.fig'],'compact');
end
end
function plotArrowsPlotOffsets()
arrowsColor = [0.4 0.4 0.4];
offSetStep = 2*1e3;
%arrows first - second
yUpper = -0.2*offSetStep;
yLower = -1.2*offSetStep;
xa = [0.6 0.6]-0.05;
doublePointedArrowInside(xa,yUpper,yLower, arrowsColor);
%arrows second - third
yUpper = -1.2*offSetStep;
yLower = -2.2*offSetStep;
xa = [0.6 0.6]-0.1;
doublePointedArrowInside(xa,yUpper,yLower, arrowsColor);
% text box with label "Plot offsets"
xa = [0.6 0]+0.1;
ya = [1e3 0];
[xaf,yaf] = axescoord2figurecoord(xa,ya);
str = {'Plot','offsets'};
annotation('textbox',[xaf(1) yaf(1)-0.02 0.1 0.1],'String',str,'FitBoxToText','On', ...
'EdgeColor','none', 'FontWeight','bold','Color',arrowsColor);
set(gcf,'Color','w')
end
createTestData/createSummedSpectrum.m
+ 90
− 2
View file @ c839aa68
function [summedFid, summedFidWithNoise, current_em, currentConcentrationRm] = createSummedSpectrum( ...
function [summedFid, summedFidWithNoise, current_em, currentConcentrationRm, stats] = createSummedSpectrum( ...
fids, concentrations, conc_std, T2s, T2_std, TE, metabolites, metaboliteFileNames, shift_t2_mx, phase_1_f2_mx, ppmMask, ...
current_pc0, current_pc1, current_df2, current_gm, current_em_std_factor, conc_std_factor, noiseVector, ...
truncPoint, saveData, paramKeyWord, indexParam, dwellTimeMs, scanFrequency, dataExportPath, referencePeakPpm, scalingFactorRefPeak, baseline)
......@@ -13,7 +13,7 @@ currentConcentrationRm = size(1,length(metabolites));
for indexMetabolite = 1:length(metabolites)
current_T2 = T2s(indexMetabolite) + current_em_std_factor(indexMetabolite) * T2_std(indexMetabolite);%in ms
if strcmp(metabolites(indexMetabolite),'MMB')
current_em(indexMetabolite) = 1;
current_em(indexMetabolite) = 0;
%the Gaussian for the acquired MMB should be smaller than the one for
%metabolites, since it was already measured with the inherent
%micro- and macro-susceptibility. It should however not be negative!
......@@ -106,6 +106,9 @@ baselineScaled = baseline * maxValueSpectrum;
summedFidWithNoise = summedFid + noiseVector + baselineScaled;
% figure; plot(real(fftshift(fft(summedFid)))); hold on; plot(real(fftshift(fft(baseline*maxValueSpectrum))))
bandwidth_Hz = 1 / dwellTimeMs * 1e3;
stats = calculateSNR(scanFrequency, bandwidth_Hz, summedFidWithNoise);
if saveData
% export the spectra but also the used concentrations and all parameters
% export not truncated data
......@@ -123,3 +126,88 @@ if saveData
'current_pc0', 'current_pc1', 'current_df2', 'current_gm', 'current_em_std_factor', 'conc_std_factor', 'baselineScaled');
end
end
function stats = calculateSNR(scanFrequency, bandwidth_Hz, fid)
noiseRangePpm = [-4, -1];
number_Pts = size(fid,2);
bw = linspace(-bandwidth_Hz/2,bandwidth_Hz/2,number_Pts);
ppmVector = bw/scanFrequency + 4.66;
indexNoiseStart = find(ppmVector>noiseRangePpm(1),1,'first');
indexNoiseEnd = find(ppmVector>noiseRangePpm(2),1,'first');
noiseRange = [indexNoiseStart indexNoiseEnd];
specReal = real(fftshift(fft(fid)));
%calculate noise
[noise_std] = noiseCalc(specReal, noiseRange);
list = {[1.9 2.1] 'NAA'; [2.9 3.1] 'Cr';[2.25 2.5] 'Glu'; [3.1 3.3] 'tCho'};
stats = cell(size(list,1)+1,1);
stats{1,1} = 'metaboliteName';
stats{1,2} = 'snr';
stats{1,3} = 'linewidth (Hz)';
stats{1,4} = 'amplitude';
for i=1:size(list,1)
metRange = list{i,1};
metRange = convertPPMtoIndex(ppmVector, metRange);
metRange =[find(metRange==1,1,'first') find(metRange==1,1,'last')];
[snrV, amplitude] = snr(specReal,metRange,noise_std);
lineWidthValue = linewidth(specReal,metRange,bandwidth_Hz);
%area = sum(data(metRange(1):metRange(2)));
stats{i+1,1} = list{i,2};
stats{i+1,2} = snrV;
stats{i+1,3} = lineWidthValue;
stats{i+1,4} = amplitude;
%stats{i,4} = area;
end
end
function [y] = convertPPMtoIndex(ppmVect, x)
%CONVERTPPMTOINDEX converts the ppm pair x to indexes
% Returns the vector of the size of the data points,
% with the selected region set to 1
index1 = ppmVect >= x(1);
index2 = ppmVect <= x(2);
y = index1 & index2;
end
function [noise_std] = noiseCalc(fftsReal, noiseRange)
% noise calculation
noise = (fftsReal(noiseRange(1):noiseRange(2)))';
fitobject = fit((0:length(noise)-1)',noise,'poly1'); %fit a poly2 to noise
y = feval(fitobject,(0:length(noise)-1)'); % substract poly1 from noise
noise = noise - y;% now the noise does not have linear component
noise_std = std(noise);
end
function [snrValue, amplitude] = snr(specReal,peakRange,noise_std)
amplitude = max(abs(specReal(peakRange(1):peakRange(2))));
snrValue = amplitude/noise_std;
%snrValue = 20*log(amplitude/noise); %in dB
end
function lineWidthValue = linewidth(specReal, metRange, bandwidth_Hz)
peakRange = metRange;
peak = specReal(peakRange(1):peakRange(2));
interpFactor = 100; %
indeces = [];
i = 1;
while (length(indeces)<2 && i <4)
peakTmp = interp(peak,interpFactor*i);
[~, maxPos] = max(abs(peakTmp));
normValues = peakTmp/(peakTmp(maxPos)/2);
indeces = find(normValues> 0.995 & normValues< 1.005);
i = i+1;
end
interpFactor = interpFactor*(i-1);
freqResolution = bandwidth_Hz/(length(specReal)*interpFactor);
try
lineWidthValue = (indeces(end) - indeces(1))*freqResolution;
catch err
lineWidthValue = [];
end
end
\ No newline at end of file
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment