From 1a38ad14279e20f2e11c3ad8fcbceef61218d04d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Thanh=20Phong=20L=C3=AA?= <43346967+tpkle@users.noreply.github.com> Date: Wed, 7 Jan 2026 20:59:14 +0100 Subject: [PATCH 1/8] Adding required FID-A functions --- code/op_addphase.m | 8 ++-- code/op_concatAverages.m | 65 ++++++++++++++++++++++++++ code/op_leftshift_keepSize.m | 44 ++++++++++++++++++ code/op_takesubspec.m | 88 ++++++++++++++++++++++++++++++++++++ 4 files changed, 200 insertions(+), 5 deletions(-) create mode 100644 code/op_concatAverages.m create mode 100644 code/op_leftshift_keepSize.m create mode 100644 code/op_takesubspec.m diff --git a/code/op_addphase.m b/code/op_addphase.m index 558a1af..b35b292 100644 --- a/code/op_addphase.m +++ b/code/op_addphase.m @@ -12,16 +12,16 @@ % ph0 = zero order phase to add (degrees) % ph1 = 1st order phase to add (in seconds); % ppm0 = (optional) frequency reference point. Default = 4.65; -% suppressPlot = (optional) Boolian to suppress plots. Default = 0; +% suppressPlot = (optional) Boolean to suppress plots. Default = 1; % % OUTPUTS: % out = Phase adjusted output spectrum. -function out=op_addphase(in,ph0,ph1,ppm0,suppressPlot); +function out=op_addphase(in,ph0,ph1,ppm0,suppressPlot) if nargin<5 - suppressPlot=0; + suppressPlot=1; if nargin<4 ppm0=4.65; if nargin<3 @@ -30,8 +30,6 @@ end end -ph0; -ph1; %Add Zero-order phase fids=in.fids.*ones(size(in.fids))*exp(1i*ph0*pi/180); diff --git a/code/op_concatAverages.m b/code/op_concatAverages.m new file mode 100644 index 0000000..375f0ab --- /dev/null +++ b/code/op_concatAverages.m @@ -0,0 +1,65 @@ +% op_concatAverages.m +% Jamie Near, McGill University 2014. +% +% USAGE: +% out=op_concatAverages(in1,in2); +% +% DESCRIPTION: +% Concatenate two scans along the averages dimension. Two scans with 50 +% averages each will now look like a single scan with 100 averages. For +% looping, it is possible to enter a blank structure as the first input. +% In this case, the function will simply return the second input. +% +% INPUTS: +% in1 = first input in matlab structure format. +% in2 = second input in matlab structure format. +% +% OUTPUTS: +% out = Output following concatenation of inputs along the +% averages dimension. + +function out=op_concatAverages(in1,in2); + + +if isempty(fieldnames(in1)) %Check if the first argument to be a null string + out=in2; %If so, return in2; +else + if in1.dims.averages~=0 && in2.dims.averages~=0 && in1.dims.averages ~= in2.dims.averages + %RULE: Averages dimension must be the same for both inputs UNLESS one of + %the inputs has an averages dimension of zero: + error('averages dimensions must be the same for both inputs'); + end + + if in1.dims.averages==0 && in2.dims.averages==0 + dim=2; + elseif in1.dims.averages==0 + dim=in2.dims.averages; + elseif in2.dims.averages==0; + dim=in1.dims.averages; + else + dim=in1.dims.averages; + end + + fids=cat(dim,in1.fids,in2.fids); + specs=cat(dim,in1.specs,in2.specs); + sz=size(fids); + + %FILLING IN DATA STRUCTURE + out=in1; + out.fids=fids; + out.specs=specs; + out.sz=sz; + out.averages=in1.averages+in2.averages; + out.rawAverages=in1.rawAverages+in2.rawAverages; + if in1.dims.averages==0 + out.dims.averages=dim; + end + + %FILLING IN THE FLAGS + out.flags=in1.flags; + out.flags.writtentostruct=1; + out.flags.averaged=0; +end + + + diff --git a/code/op_leftshift_keepSize.m b/code/op_leftshift_keepSize.m new file mode 100644 index 0000000..eda1fba --- /dev/null +++ b/code/op_leftshift_keepSize.m @@ -0,0 +1,44 @@ +% op_leftshift_keepSize.m +% Jamie Near, McGill University 2014. +% Thanh Phong Lê, 2025 +% +% USAGE: +% out=op_leftshift_keepSize(in,ls); +% +% DESCRIPTION: +% Remove leading datapoints from the fid to get rid of 1st order phase +% errors, and fill zeros at the end such that the size does not change. +% +% INPUTS: +% in = input data in matlab strucure format. +% ls = number of points to remove from the beginning of the fid. +% +% OUTPUTS: +% out = Output following left shifting. + +function out=op_leftshift_keepSize(in, ls) + +if in.flags.leftshifted + cont=input('WARNING: Left shifting has already been performed! Continue anyway? (y or n)','s'); + if cont=='y' + %continue; + else + error('STOPPING'); + end +end + +fids = in.fids; +fids(1:end-ls,:) = fids(1+ls:end,:); +fids(end-ls:end,:) = 0; + +specs=fftshift(ifft(fids,[],in.dims.t),in.dims.t); + +%FILLING IN DATA STRUCTURE +out=in; +out.fids=fids; +out.specs=specs; + +%FILLING IN THE FLAGS +out.flags=in.flags; +out.flags.writtentostruct=1; +out.flags.leftshifted=1; diff --git a/code/op_takesubspec.m b/code/op_takesubspec.m new file mode 100644 index 0000000..e2d89ec --- /dev/null +++ b/code/op_takesubspec.m @@ -0,0 +1,88 @@ +% op_takesubspec.m +% Jamie Near, McGill University 2014. +% +% USAGE: +% out=op_takesubspec(in,index); +% +% DESCRIPTION: +% Extract the subspectra with indices corresponding to the 'index' input +% array. +% +% INPUTS: +% in = input data in matlab structure format. +% index = vector indicating the indices of the subspectra you would like +% to extract. +% +% OUTPUTS: +% out = Output dataset consisting of subspectra indices extracted from +% the input. + +function out=op_takesubspec(in,index); + +if in.flags.subtracted + error('ERROR: Subspectra have already been combined! Aborting!'); +end + +if in.dims.subSpecs==0 + %Can't take subspec because there are none: + error('ERROR: There are no subspectra in this dataset! Aborting!'); +elseif in.dims.subSpecs==1 + %SHOULD NEVER HAPPEN (Time dimension should always be dim=1) + error('ERROR: dims.subSpecs==1. This should never happen! Aborting!'); +elseif in.dims.subSpecs==2 + fids=in.fids(:,index); +elseif in.dims.subSpecs==3; + fids=in.fids(:,:,index); +elseif in.dims.subSpecs==4; + fids=in.fids(:,:,:,index); +elseif in.dims.subSpecs==5 + fids=in.fids(:,:,:,:,index); +end + +%re-calculate Specs using fft +specs=fftshift(ifft(fids,[],in.dims.t),in.dims.t); + +%change the dims variables +if in.dims.t>in.dims.subSpecs + dims.t=in.dims.t-1; +else + dims.t=in.dims.t; +end +if in.dims.coils>in.dims.subSpecs + dims.coils=in.dims.coils-1; +else + dims.coils=in.dims.coils; +end +if in.dims.averages>in.dims.subSpecs + dims.averages=in.dims.averages-1; +else + dims.averages=in.dims.averages; +end +if in.dims.extras>in.dims.subSpecs + dims.extras=in.dims.extras-1; +else + dims.extras=in.dims.extras; +end +dims.subSpecs=0; + +%re-calculate the sz variable +sz=size(fids); + + +%FILLING IN DATA STRUCTURE +out=in; +out.fids=fids; +out.specs=specs; +out.sz=sz; +out.dims=dims; +out.subspecs=1; +if out.dims.averages + out.averages=out.sz(out.dims.averages); +else + out.averages=1; +end +out.rawAverages=out.averages; + +%FILLING IN THE FLAGS +out.flags=in.flags; +out.flags.writtentostruct=1; From 3dea6c0f21a8cc43b19f33194b6c475523ce02ed Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Thanh=20Phong=20L=C3=AA?= <43346967+tpkle@users.noreply.github.com> Date: Wed, 7 Jan 2026 21:01:46 +0100 Subject: [PATCH 2/8] Update of Varian data loader Output the number of points to be leftshifted to fix the 1st order phase distortion. --- code/compMRS_loadspecVarian.m | 2 +- code/dependencies/mklassenTools/readprocpar.m | 3 ++- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/code/compMRS_loadspecVarian.m b/code/compMRS_loadspecVarian.m index d59a870..5bcff2d 100644 --- a/code/compMRS_loadspecVarian.m +++ b/code/compMRS_loadspecVarian.m @@ -339,7 +339,7 @@ out.seq=sequence; out.te=te; out.tr=tr; -out.pointsToLeftshift=0; +out.pointsToLeftshift=par.lsfid; out.filepath=filename; %FILLING IN THE FLAGS diff --git a/code/dependencies/mklassenTools/readprocpar.m b/code/dependencies/mklassenTools/readprocpar.m index d2aa9a3..20fd56b 100644 --- a/code/dependencies/mklassenTools/readprocpar.m +++ b/code/dependencies/mklassenTools/readprocpar.m @@ -228,11 +228,12 @@ par.reductionFactor = []; par.referenceUnits = []; par.referenceLines = []; +par.lsfid = []; % added 20260101 Thanh names = fieldnames(par); params = fieldnames(data); -for i = 1:length(params); +for i = 1:length(params) par.(params{i}) = data.(params{i}).value; end From 85e579dddc35db9a0197a01b0f76e4aabe22c10a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Thanh=20Phong=20L=C3=AA?= <43346967+tpkle@users.noreply.github.com> Date: Wed, 7 Jan 2026 21:10:59 +0100 Subject: [PATCH 3/8] Update Bruker loader Update compMRS_parseBrukerFormat.m to properly read refscans in method file. Update compMRS_loadspecBruker.m to get the group delay in all cases (either acqp or acqus files, and from parameter GRPDLY or ACQ_RxFilterInfo). We no longer perform leftshift in compMRS_loadspecBruker; this will be done in compMRS_DPproc.m --- code/compMRS_loadspecBruker.m | 108 ++++++++++++++++++------------- code/compMRS_parseBrukerFormat.m | 24 +++++++ 2 files changed, 88 insertions(+), 44 deletions(-) diff --git a/code/compMRS_loadspecBruker.m b/code/compMRS_loadspecBruker.m index 1ed22e5..969f3f0 100644 --- a/code/compMRS_loadspecBruker.m +++ b/code/compMRS_loadspecBruker.m @@ -4,7 +4,7 @@ %Wendy Oakden, Sunnybrook Research Institute, 2020 %Jamie Near, Sunnybrook Research Institute, 2023 %Georg Oeltzschner, Johns Hopkins University 2025 -%Thanh Phong Le, EPFL, 2025 +%Thanh Phong Lê, EPFL, 2025 %Diana Rotaru, Medical University of Vienna, 2025 % % USAGE: @@ -74,6 +74,12 @@ acqpFile = fullfile(inDir, 'acqp'); headerACQP = compMRS_parseBrukerFormat(acqpFile); +% Populate the header information from the ACQUS file +acqusFile = fullfile(inDir, 'acqus'); +if isfile(acqusFile) + headerACQUS = compMRS_parseBrukerFormat(acqusFile); +end + % Populate the header information from the Method file methodFile = fullfile(inDir, 'method'); headerMethod = compMRS_parseBrukerFormat(methodFile); @@ -168,24 +174,27 @@ end % Group delay (i.e. left shift) -if contains(version, ["PV 5", "PV 6", "PV 7", "PV-7"]) - leftshift = round(headerACQP.GRPDLY); -elseif contains(version,'PV-360') - % I'm not sure this differentiation is necessary but Jessie is getting - % GRPDELY from the ACQUS file - might remove if not needed - if exist('headerACQUS', 'var') - leftshift = round(headerACQUS.GRPDLY); +leftshift=-1; +if isfield(headerACQP,'GRPDLY') + leftshift = headerACQP.GRPDLY; +end + +if leftshift==-1 && exist('headerACQUS', 'var') + leftshift = headerACQUS.GRPDLY; +end + +if leftshift==-1 && isfield(headerACQP, 'ACQ_RxFilterInfo') + if ~iscell(headerACQP.ACQ_RxFilterInfo{1}) + leftshift = str2double(headerACQP.ACQ_RxFilterInfo{1}); else - leftshift = round(headerACQP.GRPDLY); - end - - % Found a couple instances of GRPDLY being -1, in that case, set it to - % the PV360 default which is 77 (from Brayan Alves' script) - if leftshift == -1 - leftshift = 77; + leftshift = str2double(headerACQP.ACQ_RxFilterInfo{1}{1}); end end +if leftshift==-1 + disp('NO GROUP DELAY FOUND') +end + % Spectral width if contains(version, ["PV 5", "PV 6", "PV 7", "PV-7"]) spectralwidth = headerMethod.PVM_DigSw; @@ -419,13 +428,13 @@ error('ERROR: rawData variable not recognized. Options are ''y'' or ''n''.'); end - -%Perform left-shifting to remove points before the echo -fids_trunc=fids_raw(leftshift+1:end,:,:,:); - -%replace the left-shifted points with zeros at the end -fids=padarray(fids_trunc, [leftshift,0],'post'); - +% Removed: we will perform leftshift in compMRS_DPproc - Thanh 20260101 +% %Perform left-shifting to remove points before the echo +% fids_trunc=fids_raw(leftshift+1:end,:,:,:); +% +% %replace the left-shifted points with zeros at the end +% fids=padarray(fids_trunc, [leftshift,0],'post'); +fids=fids_raw; sz=size(fids); %size of the array %calculate the dwelltime: @@ -471,12 +480,19 @@ end if strcmpi(rawData,'y') - dims.averages=2; + if rawRepetitions>1 - dims.subSpecs=4; + if rawAverages>1 + dims.averages=2; + dims.subSpecs=4; + else + dims.averages=0; + dims.subSpecs=3; + end else dims.subSpecs=0; end + elseif strcmpi(rawData,'n') dims.averages=0; if rawRepetitions>1 @@ -516,14 +532,16 @@ % The data stored are already averaged. averages_ref=1; ref.flags.averaged=1; - - fids_ref=reshape(fids_ref,[], Nrcvrs); - fids_ref_trunc=fids_ref(leftshift+1:end,:,:); - reffids=padarray(fids_ref_trunc, [leftshift,0],'post'); - + + reffids=reshape(fids_ref,[], Nrcvrs); + + % Removed: we will perform leftshift in compMRS_DPproc - Thanh 20260101 + % fids_ref_trunc=fids_ref(leftshift+1:end,:,:); + % reffids=padarray(fids_ref_trunc, [leftshift,0],'post'); + % Apply ref freq shift (the difference between txfrq and txfrq_ref) % (but not if it's PV360) - sz_ref=size(fids_ref); %size of the array + sz_ref=size(reffids); %size of the array if ~contains(version, ["PV-360"]) tmat=repmat(t',[1 sz_ref(2:end)]); reffids=reffids.*exp(-1i*tmat*(txfrq_ref-txfrq)*2*pi); @@ -585,13 +603,14 @@ rawAverages_ref=1; end - fids_ref=reshape(fids_ref,[],averages_ref); - fids_ref_trunc=fids_ref(leftshift+1:end,:); - reffids=padarray(fids_ref_trunc, [leftshift,0],'post'); - + reffids=reshape(fids_ref,[],averages_ref); + % Removed: we will perform leftshift in compMRS_DPproc - Thanh 20260101 + % fids_ref_trunc=fids_ref(leftshift+1:end,:); + % reffids=padarray(fids_ref_trunc, [leftshift,0],'post'); + % % Apply ref freq shift (the difference between txfrq and txfrq_ref) % (but not if it's PV360) - sz_ref=size(fids_ref); %size of the array + sz_ref=size(reffids); %size of the array if ~contains(version, ["PV-360"]) tmat=repmat(t',[1 sz_ref(2:end)]); reffids=reffids.*exp(-1i*tmat*(txfrq_ref-txfrq)*2*pi); @@ -670,9 +689,10 @@ averages_nav=rawAverages_nav; if contains(version,'PV 5') - fids_nav=reshape(fids_nav,[],rawAverages_nav); - fids_nav_trunc=fids_nav(leftshift+1:end,:); - navfids=padarray(fids_nav_trunc, [leftshift,0],'post'); + navfids=reshape(fids_nav,[],rawAverages_nav); + % Removed: we will perform leftshift in compMRS_DPproc - Thanh 20260101 + % fids_nav_trunc=fids_nav(leftshift+1:end,:); + % navfids=padarray(fids_nav_trunc, [leftshift,0],'post'); else navfids=reshape(fids_nav,headerMethod.PVM_NavPoints,Nrcvrs, rawAverages, rawRepetitions); %Permute so that time is along 1st dimension, averages is along 2nd @@ -756,7 +776,7 @@ out.seq=sequence; out.te=te; out.tr=tr; -out.pointsToLeftshift=0; +out.pointsToLeftshift=leftshift; out.version=version; out.filepath=fileRaw; @@ -764,7 +784,7 @@ %FILLING IN THE FLAGS FOR THE RAW DATA out.flags.writtentostruct=1; out.flags.gotparams=1; -out.flags.leftshifted=1; +out.flags.leftshifted=0; out.flags.filtered=0; out.flags.zeropadded=0; out.flags.freqcorrected=0; @@ -808,7 +828,7 @@ ref.seq=sequence; ref.te=te; ref.tr=tr; - ref.pointsToLeftshift=0; + ref.pointsToLeftshift=leftshift; ref.version=version; % DGR added file path to FID-A structure; % when no separate water scan is acquired (a reference scan) @@ -823,7 +843,7 @@ %FILLING IN THE FLAGS FOR THE REF DATA ref.flags.writtentostruct=1; ref.flags.gotparams=1; - ref.flags.leftshifted=1; + ref.flags.leftshifted=0; ref.flags.filtered=0; ref.flags.zeropadded=0; ref.flags.freqcorrected=0; @@ -869,13 +889,13 @@ nav.seq=sequence; nav.te=te; nav.tr=tr; - nav.pointsToLeftshift=0; + nav.pointsToLeftshift=leftshift; nav.version=version; %FILLING IN THE FLAGS FOR THE NAV DATA nav.flags.writtentostruct=1; nav.flags.gotparams=1; - nav.flags.leftshifted=1; + nav.flags.leftshifted=0; nav.flags.filtered=0; nav.flags.zeropadded=0; nav.flags.freqcorrected=0; diff --git a/code/compMRS_parseBrukerFormat.m b/code/compMRS_parseBrukerFormat.m index 375509f..696be3b 100644 --- a/code/compMRS_parseBrukerFormat.m +++ b/code/compMRS_parseBrukerFormat.m @@ -124,6 +124,13 @@ % In this case, let's look for recurring parentheses % again: multiLine = erase(multiLine, newline); % remove new line characters + + % For some parameters (for example the reference scan), + % we need to "uncompress" the data + + multiLine = replacePattern(multiLine); + + [tokensParentheses, ~] = regexp(multiLine,'\(([^\)]+)\)','tokens','match'); if ~isempty(tokensParentheses) @@ -201,4 +208,21 @@ fclose(fid); +end + +function out = replacePattern(str) + +expr = '@(\d+)\*\((\d+)\)'; +tokens = regexp(str, expr, 'tokens'); +matches = regexp(str, expr, 'match'); + +out = str; + +for k = 1:numel(matches) + count = str2double(tokens{k}{1}); + value = tokens{k}{2}; + replacement = strjoin(repmat({value}, 1, count), ' '); + out = strrep(out, matches{k}, replacement); +end + end \ No newline at end of file From 80d5a2591e556c2615b9be988c85f7a3f70b6173 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Thanh=20Phong=20L=C3=AA?= <43346967+tpkle@users.noreply.github.com> Date: Wed, 7 Jan 2026 21:17:04 +0100 Subject: [PATCH 4/8] Update of compMRS_DPproc and compMRS_DPload compMRS_DPload is modified to make it compatible with DP36 (no separate refscan, only automatic refscan). Update of compMRS_DPproc. Most DPs can be processed. There are still issues with SPECIAL acquisitions on Bruker systems. The first-order phase correction in case of a non-integer of left-shifted points still has to be checked. New compMRS_runDPproc_onAllDP.m script to process all DPs at once. --- code/compMRS_DPload.m | 38 ++--- code/compMRS_DPproc.m | 308 ++++++++++++++++++++++++++++++----- code/compMRS_DPproc_allDPs.m | 36 ++++ 3 files changed, 314 insertions(+), 68 deletions(-) create mode 100644 code/compMRS_DPproc_allDPs.m diff --git a/code/compMRS_DPload.m b/code/compMRS_DPload.m index e48fbb9..39c20b8 100644 --- a/code/compMRS_DPload.m +++ b/code/compMRS_DPload.m @@ -52,38 +52,22 @@ for n = 1:check.nSes(m) sess=dir([DPid filesep subjs(m).name filesep 'ses*']); - %Find the MRS data path and the REF data path: + %Find the MRS data path and the REF data path (if applicable): svspath = dir([DPid filesep subjs(m).name filesep sess(n).name filesep 'mrs' filesep '*svs']); refpath = dir([DPid filesep subjs(m).name filesep sess(n).name filesep 'mrs' filesep '*mrsref']); - % need to find out whether the ref data is acquired separately - % or included with the metabolite data - % For PV6 and higher where only one directory is expected to - % exist if water data was acquired automatically, or two - % directories if the water data was acquired separately + % Only one directory is expected to exist if water data was + % acquired automatically, or two directories if the water data + % was acquired separately - % The syntax is different for ParaVision 7: it is usually - % 'PV-7' while in version 6 it is PV 6 with a space... - if contains(check.version{1}, ["PV 6", "PV 7", "PV-7", "PV-360", "PV 360"]) - - % refscan may be acquired automatically - [out{m,n}, outw_auto{m,n}]=compMRS_loadspecBruker([svspath(length(svspath)).folder filesep svspath(length(svspath)).name],'y'); - - % If there was a refscan acquired separatly, load it as well - if ~isempty(svspath) && ~isempty(refpath) % refscan acquired separately - [outw{m,n}]=compMRS_loadspecBruker([refpath(length(refpath)).folder filesep refpath(length(refpath)).name],'y'); - end - - %If PV5, then the water reference scans must be collected - %separately: - elseif str2num(extractBetween(check.version{m,n}, 'PV ', '.')) == 5 - if exist([svspath(length(svspath)).folder filesep svspath(length(svspath)).name]) && exist([refpath(length(refpath)).folder filesep refpath(length(refpath)).name]) % refscan acquired separately - [out{m,n}]=compMRS_loadspecBruker([svspath(length(svspath)).folder filesep svspath(length(svspath)).name],'y'); - [outw{m,n}]=compMRS_loadspecBruker([refpath(length(refpath)).folder filesep refpath(length(refpath)).name],'y'); - else - error(['ERROR: For PV5, must have both svs and ref directories!! Aborting!! ']) - end + % refscan may be acquired automatically + [out{m,n}, outw_auto{m,n}]=compMRS_loadspecBruker([svspath(length(svspath)).folder filesep svspath(length(svspath)).name],'y'); + + % If there was a refscan acquired separatly, load it as well + if ~isempty(svspath) && ~isempty(refpath) % refscan acquired separately + [outw{m,n}]=compMRS_loadspecBruker([refpath(length(refpath)).folder filesep refpath(length(refpath)).name],'y'); end + end end elseif strcmp(check.vendor(1),'VARIAN') diff --git a/code/compMRS_DPproc.m b/code/compMRS_DPproc.m index 096538d..7c30825 100644 --- a/code/compMRS_DPproc.m +++ b/code/compMRS_DPproc.m @@ -1,47 +1,69 @@ %compMRS_DPproc.m - +% %Jamie Near, Sunnybrook Research Institute, 2025 %Diana Rotaru, Columbia University, 2025 %Thanh Phong Lê, CIBM Center for Biomedical Imaging and École polytechnique fédérale de Lausanne, 2025 % - - -% Input: -% Output: - % USAGE: +% [out,outw]=compMRS_DPproc(DPid); +% +% DESCRIPTION: +% +% +% INPUTS: +% DPid: Data Packet ID (i.e. DP01, DP02, etc.) +% opt: Optional parameters: +% +% +% OUTPUTS: +% out: m x n cell array where m is the number of subjects in the DP, +% and n is the number of sessions in the DP. Each element of the +% cell array is a water suppressed FID-A data struct. +% outw: m x n cell array where m is the number of subjects in the DP, +% and n is the number of sessions in the DP. Each element of the +% cell array is a water unsuppressed FID-A data struct. +% +% Both outputs are spectrally processed. out{m, n} contains the linewidth and SNR. -% Add the 'code' folder to path then run this code from the data folder. - -% [out,outw]=compMRS_DPproc.m - - -function [out, outw] = compMRS_DPproc(DPid) +function [out, outw] = compMRS_DPproc(DPid, opt) +disp(['Processing ' num2str(DPid)]) -mkdir('plots'); +if ~exist('opt','var') + % Default options for the processing + opt.doECCbeforeAvg = 1; % This is how it usually done with Bruker data + opt.predefCoilAmpl = 1; % Use the predefined coil scaling coefficients, when available + opt.rmBadAvg = 0; + opt.doBlockAveraging = 0; + opt.doDriftCorrection = 0; + opt.iterin = 20; + opt.tmaxin = 0.2; + opt.aaDomain = 'f'; +end -% Options for the processing -opt.doECCbeforeAvg = 1; % This is how it usually done with Bruker data -opt.predefCoilAmpl = 1; % Use the predefined coil scaling coefficients, when available -opt.rmBadAvg = 0; % try + [check]=compMRS_DPcheck(DPid); if check.allSame + + % Create a folder to save the plots + if ~isfolder([pwd filesep 'plots']) + mkdir('plots') + end + [in, inw, inw_auto] = compMRS_DPload(DPid); %Loop through subjects and sessions. for m = 1:check.nSubj for n = 1:check.nSes(m) - disp(['Processing ' DPid ' sub-' num2str(m) ' ses-' num2str(n)]) - + % If separate water scan is available if ~isempty(inw) && ~isempty(inw{m, n}) ident = [DPid '_sub-' num2str(m) '_ses-' num2str(n) '_sepWater']; [out{m, n}, outw{m, n}] = compMRS_DPproc_sub(in{m, n}, inw{m, n}, ident, opt); - + % If automatic water scan is available elseif ~isempty(inw_auto) && ~isempty(inw_auto{m, n}) ident = [DPid '_sub-' num2str(m) '_ses-' num2str(n) '_autoWater']; @@ -57,9 +79,16 @@ end function [out, outw] = compMRS_DPproc_sub(in_mn, inw_mn, ident, opt) + + ls = in_mn.pointsToLeftshift; + frac_ls = ls-floor(ls); + + out_mn = op_leftshift_keepSize(in_mn, floor(ls)); + out_wmn= op_leftshift_keepSize(inw_mn, floor(ls)); + % do ECC before averaging (yes/no) if opt.doECCbeforeAvg - [out_mn, outw_mn]=op_eccKlose(in_mn,inw_mn); + [out_mn, outw_mn]=op_eccKlose(out_mn,out_wmn); end % Do coil combination WITHOUT averaging (if applicable) @@ -97,32 +126,51 @@ end end - % do bad averages removal (yes/no) - not implemented yet + % Process the water + outw_mn=op_alignAverages(outw_mn,0.2,'n'); + outw_mn=op_averaging(outw_mn); + + % do bad averages removal, if applicable (code from Jamie) if opt.rmBadAvg - % op_rmbadaverages + outw_mn=subBadAveragesRemoval(outw_mn, opt); end - % Perform partial averaging + % Perform partial averaging, if applicable + if opt.doBlockAveraging + % Perform partial averaging - % Find the possible blocks sizes for partial averaging - av_block_sizes = divisors(out_mn.averages); + % Find the possible blocks sizes for partial averaging + av_block_sizes = divisors(out_mn.averages); + + % Some data are already partially averaged (Varian datasets) + av_eff_block_sizes = av_block_sizes*out_mn.rawAverages/out_mn.averages; - % Some data are already partially averaged (Varian datasets) - av_eff_block_sizes = av_block_sizes*out_mn.rawAverages/out_mn.averages; - - % Remove block sizes bigger than 32 (does not really make sense to go higher - % than that) - av_block_sizes =av_block_sizes(av_eff_block_sizes<=32); - av_eff_block_sizes=av_eff_block_sizes(av_eff_block_sizes<=32); + % Remove effective block sizes bigger than 32 (does not really make + % sense to go higher than that) + av_block_sizes =av_block_sizes(av_eff_block_sizes<=32); + av_eff_block_sizes=av_eff_block_sizes(av_eff_block_sizes<=32); + else + % Do not perform block averaging + av_block_sizes = 1; + av_eff_block_sizes= out_mn.rawAverages/out_mn.averages; + end - % Iterate along the list of block sizes. The first element is without - % block averaging + % Iterate along the list of block sizes, if applicable + for kk=1:length(av_block_sizes) out_part_avg = op_blockAvg(out_mn,av_block_sizes(kk)); - % do drift correction (if applicable) - [out_part_avg, fs]=op_freqAlignAverages(out_part_avg, 1, 'n'); + % do drift correction (if applicable) (code from Jamie) + + + if opt.doDriftCorrection + out_part_avg = subDriftCorrection(out_part_avg, ident, opt); + end + + + + % do averaging (if applicable) out_part_avg=op_averaging(out_part_avg); @@ -131,8 +179,15 @@ [out_part_avg, ~]=op_eccKlose(out_part_avg, outw_mn); end + compensateFractionalGroupDelay=1; + if compensateFractionalGroupDelay + ph1 = -frac_ls*in_mn.dwelltime; + out_part_avg=op_addphase(out_part_avg,0,ph1,4.65,1); + end + % Compute the quality metrics % Get LW (of NAA) and SNR + out_part_avg = op_autophase(out_part_avg, 1.7, 2.3); [FWHM_NAA] = op_getLW(out_part_avg, 1.8, 2.2, 8, 1); [SNR]=op_getSNR(out_part_avg,1.8,2.2,-2, 0, 1); @@ -144,13 +199,18 @@ out_all{kk}=out_part_avg; end - % Search for the best output + % Search for the best output, if we performed block averaging + % (otherwise return trivial result) SNR_LW_ratios = zeros(length(out_all), 0); for kk=1:length(out_all) SNR_LW_ratios(kk) = out_all{kk}.SNR_LW_ratio; end + [~, index] = max(SNR_LW_ratios); - disp(['The best result is with block size ' num2str(av_eff_block_sizes(index)) '.']) + + if opt.doBlockAveraging + disp(['The best result is with block size ' num2str(av_eff_block_sizes(index)) '.']) + end % Output the output with best SNR/LW out = out_all{index}; @@ -159,7 +219,7 @@ % Plot and save the result to check plotlegend = {}; for ii=1:length(out_all) - plotlegend{ii} = [num2str(out_all{ii}.block_size) ' blocks']; + plotlegend{ii} = [num2str(out_all{ii}.block_size) ' avg/block']; end f=figure ('name', ident); hold on @@ -176,4 +236,170 @@ legend(plotlegend) saveas(f, ['plots/' ident], 'fig') print(f, ['plots/' ident], '-dpng', '-r600'); -end \ No newline at end of file + +end + + +function output = subBadAveragesRemoval(input, ident, opt) + nBadAvgTotal=0; + nBadAverages=1; + + sat='n'; + output = input; + while sat=='n' || sat=='N' + nsd=4; %Setting the number of standard deviations; + iter=1; + nBadAverages=1; + nBadAvgTotal=0; + while nBadAverages>0 + [output,metric{iter},badAverages]=op_rmbadaverages(output,nsd,'t'); + nBadAverages=length(badAverages); + nBadAvgTotal=nBadAvgTotal+nBadAverages; + iter=iter+1; + disp([num2str(nBadAverages) ' bad averages removed on this iteration.']); + disp([num2str(nBadAvgTotal) ' bad averages removed in total.']); + close all; + end + %figure('position',[0 50 560 420]); + %Make figure to show pre-post removal of averages + h=figure('visible','off'); + subplot(1,2,1); + plot(input.ppm,real(input.specs(:,:)));xlim([1 5]); + set(gca,'FontSize',8); + set(gca,'XDir','reverse'); + xlabel('Frequency (ppm)','FontSize',10); + ylabel('Amplitude(a.u.)','FontSize',10); + title('Before','FontSize',12); + subplot(1,2,2); + plot(output.ppm,real(output.specs(:,:)));xlim([1 5]); + set(gca,'FontSize',8); + set(gca,'XDir','reverse'); + xlabel('Frequency (ppm)','FontSize',10); + ylabel('Amplitude(a.u.)','FontSize',10); + title('After','FontSize',12); + set(h,'PaperUnits','centimeters'); + set(h,'PaperPosition',[0 0 20 15]); + saveas(h,['plots' filesep ident '_rmBadAvg_prePostFig'],'jpg'); + saveas(h,['plots' filesep ident '_rmBadAvg_prePostFig'],'fig'); + close(h); + + %figure('position',[0 550 560 420]); + h=figure('visible','off'); + plot([1:length(metric{1})],metric{1},'.r',[1:length(metric{iter-1})],metric{iter-1},'x','MarkerSize',16); + set(gca,'FontSize',8); + xlabel('Scan Number','FontSize',10); + ylabel('Deviation Metric','FontSize',10); + legend('Before rmBadAv','After rmBadAv'); + legend boxoff; + title('Deviation Metric','FontSize',12); + set(h,'PaperUnits','centimeters'); + set(h,'PaperPosition',[0 0 20 10]); + saveas(h,['plots' filesep ident '_rmBadAvg_scatterFig'],'png'); + saveas(h,['plots' filesep ident '_rmBadAvg_scatterFig'],'fig'); + close(h); + + %sat1=input('are you satisfied with the removal of bad averages? ','s'); + sat='y'; + end +end + + +function output = subDriftCorrection(input, ident, opt); + sat='n'; + out_rm2=input; + while sat=='n' || sat=='N' + fsPoly=100; + phsPoly=1000; + fscum=zeros(out_rm2.sz(out_rm2.dims.averages),1); + phscum=zeros(out_rm2.sz(out_rm2.dims.averages),1); + iter=1; + while (abs(fsPoly(1))>0.001 || abs(phsPoly(1))>0.01) && iter Date: Fri, 9 Jan 2026 12:34:19 +0100 Subject: [PATCH 5/8] Correction of small bugs All datasets can be processed now (there are still some phasing issues), except DP22 and DP23 (different sequence, and how data are arranged). --- code/compMRS_DPproc.m | 23 +++++---- code/compMRS_loadspecBruker.m | 47 ++++++++++-------- code/op_combinesubspecs.m | 92 +++++++++++++++++++++++++++++++++++ 3 files changed, 129 insertions(+), 33 deletions(-) create mode 100644 code/op_combinesubspecs.m diff --git a/code/compMRS_DPproc.m b/code/compMRS_DPproc.m index 7c30825..a51a27e 100644 --- a/code/compMRS_DPproc.m +++ b/code/compMRS_DPproc.m @@ -31,8 +31,8 @@ if ~exist('opt','var') % Default options for the processing - opt.doECCbeforeAvg = 1; % This is how it usually done with Bruker data - opt.predefCoilAmpl = 1; % Use the predefined coil scaling coefficients, when available + opt.doECCbeforeAvg = 0; % 1 is how it usually done with Bruker data + opt.predefCoilAmpl = 0; % Use the predefined coil scaling coefficients, when available opt.rmBadAvg = 0; opt.doBlockAveraging = 0; opt.doDriftCorrection = 0; @@ -63,9 +63,9 @@ if ~isempty(inw) && ~isempty(inw{m, n}) ident = [DPid '_sub-' num2str(m) '_ses-' num2str(n) '_sepWater']; [out{m, n}, outw{m, n}] = compMRS_DPproc_sub(in{m, n}, inw{m, n}, ident, opt); - + end % If automatic water scan is available - elseif ~isempty(inw_auto) && ~isempty(inw_auto{m, n}) + if ~isempty(inw_auto) && ~isempty(inw_auto{m, n}) ident = [DPid '_sub-' num2str(m) '_ses-' num2str(n) '_autoWater']; [out{m, n}, outw{m, n}] = compMRS_DPproc_sub(in{m, n}, inw_auto{m, n}, ident, opt); end @@ -84,11 +84,13 @@ frac_ls = ls-floor(ls); out_mn = op_leftshift_keepSize(in_mn, floor(ls)); - out_wmn= op_leftshift_keepSize(inw_mn, floor(ls)); + outw_mn= op_leftshift_keepSize(inw_mn, floor(ls)); + + % Here check whether subspectra have to be combined first. % do ECC before averaging (yes/no) if opt.doECCbeforeAvg - [out_mn, outw_mn]=op_eccKlose(out_mn,out_wmn); + [out_mn, outw_mn]=op_eccKlose(out_mn,outw_mn); end % Do coil combination WITHOUT averaging (if applicable) @@ -122,7 +124,7 @@ if contains(in_mn.version, ["PV 360", "PV-360"]) % divide by number of channels to achieve averaging out_mn = op_ampScale(out_mn, 1.0/length(coilcombos_to_apply.ph)); - outw_mn = op_ampScale(ref_proc, 1.0/length(coilcombos_to_apply.ph)); + outw_mn = op_ampScale(outw_mn, 1.0/length(coilcombos_to_apply.ph)); end end @@ -168,9 +170,6 @@ out_part_avg = subDriftCorrection(out_part_avg, ident, opt); end - - - % do averaging (if applicable) out_part_avg=op_averaging(out_part_avg); @@ -179,10 +178,10 @@ [out_part_avg, ~]=op_eccKlose(out_part_avg, outw_mn); end - compensateFractionalGroupDelay=1; + compensateFractionalGroupDelay=1; if compensateFractionalGroupDelay ph1 = -frac_ls*in_mn.dwelltime; - out_part_avg=op_addphase(out_part_avg,0,ph1,4.65,1); + out_part_avg=op_addphase(out_part_avg,0,ph1,4.7,1); end % Compute the quality metrics diff --git a/code/compMRS_loadspecBruker.m b/code/compMRS_loadspecBruker.m index 969f3f0..7133c09 100644 --- a/code/compMRS_loadspecBruker.m +++ b/code/compMRS_loadspecBruker.m @@ -211,7 +211,7 @@ txfrq = headerMethod.PVM_FrqWork(1)*1e6; txfrq_ref = headerMethod.PVM_FrqRef(1)*1e6; end -Bo=txfrq/42577000; +Bo=txfrq_ref/42577000; % spectralwidthppm is stored in the header BUT it is slightly different % than if we calculate it from txfrq and spectralwidth - this suggests we @@ -222,15 +222,14 @@ spectralwidthppm=spectralwidth/(txfrq_ref/1e6); % old method spectralwidthppm=headerMethod.PVM_SpecSW; % better method? -% Center frequency is *explicitly* stored in the headers from PV6 onwards - -% mostly it's been 4.7 so we'll assume that for PV5 too. +% Center frequency is *explicitly* stored in the headers from PV6 onwards % For compatibility with the other FID-A functions, we should probably % store centerfreq in the FID-A header? -if contains(version, ["PV 5", "PV 6", "PV 7", "PV-7"]) - centerfreq = 4.7; +if contains(version, ["PV 5"]) + centerfreq = 4.7 + headerMethod.PVM_SpecOffsetppm; centerfreq_ref = centerfreq; -elseif contains(version,'PV-360') - centerfreq = headerMethod.PVM_FrqRefPpm(1); +elseif contains(version,["PV 6", "PV 7", "PV-7" "PV-360"]) + centerfreq = 4.7; %headerMethod.PVM_FrqWorkPpm(1); end % Receiver gain which we seem to need when reconciling raw and processed data @@ -356,19 +355,22 @@ %dimension, and coils is along 3rd dimension: fids_raw=permute(fids_raw,[1,3,2,4]); - elseif contains(version,'PV 5') && multiRcvrs + elseif contains(version,'PV 5') % && multiRcvrs % PV 5 does not appear to store individual coil channels in the % job0 file, but still preserves individual transients % Reshape into a Npts x Naverages x Nrepetitions array fids_raw=reshape(fids_raw,rawDataPoints,rawAverages,rawRepetitions); - elseif ~contains(version,'PV 5') && ~multiRcvrs - % Found a dataset with only one receiver and one repetition - fids_raw=reshape(fids_raw,rawDataPoints,rawAverages,rawRepetitions); + % elseif ~contains(version,'PV 5') && ~multiRcvrs + % % Found a dataset with only one receiver and one repetition + % fids_raw=reshape(fids_raw,rawDataPoints,rawAverages,rawRepetitions); end + + + % Remove singletons fids_raw = squeeze(fids_raw); @@ -449,7 +451,7 @@ ppm=f/(txfrq/1e6)+centerfreq; % Apply ref freq shift (the difference between txfrq and txfrq_ref) -if contains(version, ["PV-360"]) || (contains(version, ["PV 6.0.1", "PV-7"]) && strcmpi(rawData,'y')) +if contains(version, ["PV-360"]) || (contains(version, ["PV 6", "PV-7"]) && strcmpi(rawData,'y')) tmat=repmat(t',[1 sz(2:end)]); fids=fids.*exp(-1i*tmat*(txfrq_ref-txfrq)*2*pi); end @@ -480,13 +482,15 @@ end if strcmpi(rawData,'y') - + if rawAverages>1 + dims.averages=2; + else + dims.averages=0; + end if rawRepetitions>1 if rawAverages>1 - dims.averages=2; dims.subSpecs=4; else - dims.averages=0; dims.subSpecs=3; end else @@ -508,6 +512,7 @@ rawSubspecs=rawRepetitions; + %% NOW TRY LOADING IN THE RAW REFERENCE SCAN DATA (IF IT EXISTS) % Added by Thanh 18.04.2025, read the RefScan from method file @@ -732,12 +737,12 @@ %Do the fourier transform specs_nav=fftshift(ifft(navfids,[],1),1); -% % Apply ref freq shift (the difference between txfrq and txfrq_ref) -% % (but not if it's PV360) -% if ~contains(version, ["PV-360"]) -% tmat=repmat(t',[1 sz(2:end)]); -% navfids=navfids.*exp(-1i*tmat*(txfrq_ref-txfrq)*2*pi); -% end + % Apply ref freq shift (the difference between txfrq and txfrq_ref) + % (but not if it's PV360) + % if ~contains(version, ["PV-360"]) + % tmat=repmat(t',[1 sz(2:end)]); + % navfids=navfids.*exp(-1i*tmat*(txfrq_ref-txfrq)*2*pi); + % end nav.flags.averaged=0; %specify the dims diff --git a/code/op_combinesubspecs.m b/code/op_combinesubspecs.m new file mode 100644 index 0000000..2a45ae7 --- /dev/null +++ b/code/op_combinesubspecs.m @@ -0,0 +1,92 @@ +% op_combinesubspecs.m +% Jamie Near, McGill University 2014. +% +% USAGE: +% out=op_combinesubspecs(in,mode); +% +% DESCRIPTION: +% Combine the subspectra in an acquisition either by addition or +% subtraction. +% +% INPUTS: +% in = input data in matlab structure format. +% mode = -"diff" adds the subspectra together (this is counter intuitive, +% but the reason is that many "difference editing" sequences use phase +% cycling of the readout ADC to achieve "subtraction by addition". +% -"summ" performs a subtraction of the subspectra. +% +% OUTPUTS: +% out = Output following combination of subspectra. + +function out=op_combinesubspecs(in,mode) + +if in.flags.subtracted + error('ERROR: Subspectra have already been combined! Aborting!'); +end +if in.flags.isFourSteps + error('ERROR: data with four steps must first be converted using op_fourStepCombine.m! Aborting!'); +end + +% if ~in.flags.freqcorrected +% disp('WARNING: Frequency correction has not yet been performed!'); +% end +% if ~in.flags.phasecorrected +% disp('WARNING: Phase correction has not yet been performed!'); +% end + + +switch mode + case 'diff' + %add the spectrum along the subSpecs dimension; + fids=sum(in.fids,in.dims.subSpecs); + case 'summ' + %subtract the spectrum along the subSpecs dimension; + fids=diff(in.fids,1,in.dims.subSpecs); +end + +fids=fids/in.sz(in.dims.subSpecs); %divide by nymber of subspecs so that this is an averaging operation; +fids=squeeze(fids); + +%re-calculate Specs using fft +specs=fftshift(ifft(fids,[],in.dims.t),in.dims.t); + +%change the dims variables +if in.dims.t>in.dims.subSpecs + dims.t=in.dims.t-1; +else + dims.t=in.dims.t; +end +if in.dims.coils>in.dims.subSpecs + dims.coils=in.dims.coils-1; +else + dims.coils=in.dims.coils; +end +if in.dims.averages>in.dims.subSpecs + dims.averages=in.dims.averages-1; +else + dims.averages=in.dims.averages; +end +dims.subSpecs=0; +if in.dims.extras>in.dims.subSpecs + dims.extras=in.dims.extras-1; +else + dims.extras=in.dims.extras; +end + +%re-calculate the sz variable +sz=size(fids); + + +%FILLING IN DATA STRUCTURE +out=in; +out.fids=fids; +out.specs=specs; +out.sz=sz; +out.dims=dims; +out.subspecs=1; +out.averages=in.averages/2; + +%FILLING IN THE FLAGS +out.flags=in.flags; +out.flags.writtentostruct=1; +out.flags.subtracted=1; From 57a449d0bb05a3de1aa97f0c0c889b3ce0048c39 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Thanh=20Phong=20L=C3=AA?= <43346967+tpkle@users.noreply.github.com> Date: Fri, 9 Jan 2026 14:14:34 +0100 Subject: [PATCH 6/8] Correction of another bug --- code/compMRS_loadspecBruker.m | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/code/compMRS_loadspecBruker.m b/code/compMRS_loadspecBruker.m index 7133c09..050ce62 100644 --- a/code/compMRS_loadspecBruker.m +++ b/code/compMRS_loadspecBruker.m @@ -353,7 +353,11 @@ fids_raw=reshape(fids_raw,rawDataPoints,Nrcvrs,rawAverages,rawRepetitions); %Permute so that time is along 1st dimension, averages is along 2nd %dimension, and coils is along 3rd dimension: - fids_raw=permute(fids_raw,[1,3,2,4]); + fids_raw=permute(fids_raw,[1,3,2,4]); + + elseif ~contains(version,'PV 5') && ~multiRcvrs + % Found a dataset with only one receiver and one repetition + fids_raw=reshape(fids_raw,rawDataPoints,rawAverages,rawRepetitions); elseif contains(version,'PV 5') % && multiRcvrs % PV 5 does not appear to store individual coil channels in the @@ -362,9 +366,7 @@ % Reshape into a Npts x Naverages x Nrepetitions array fids_raw=reshape(fids_raw,rawDataPoints,rawAverages,rawRepetitions); - % elseif ~contains(version,'PV 5') && ~multiRcvrs - % % Found a dataset with only one receiver and one repetition - % fids_raw=reshape(fids_raw,rawDataPoints,rawAverages,rawRepetitions); + end From 4a218ac039eedd0ae03b2b9d374a882e3d73e271 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Thanh=20Phong=20L=C3=AA?= <43346967+tpkle@users.noreply.github.com> Date: Fri, 16 Jan 2026 16:35:14 +0100 Subject: [PATCH 7/8] Adding support for SPECIAL and correct frequency shifts --- code/compMRS_DPproc.m | 71 ++++++++++++-------- code/compMRS_DPproc_allDPs.m | 4 +- code/compMRS_loadspecBruker.m | 122 +++++++++++++++++++++------------- 3 files changed, 120 insertions(+), 77 deletions(-) diff --git a/code/compMRS_DPproc.m b/code/compMRS_DPproc.m index a51a27e..747bac4 100644 --- a/code/compMRS_DPproc.m +++ b/code/compMRS_DPproc.m @@ -26,19 +26,21 @@ % Both outputs are spectrally processed. out{m, n} contains the linewidth and SNR. -function [out, outw] = compMRS_DPproc(DPid, opt) +function [out, outw, out_auto, outw_auto] = compMRS_DPproc(DPid, opt) disp(['Processing ' num2str(DPid)]) if ~exist('opt','var') % Default options for the processing - opt.doECCbeforeAvg = 0; % 1 is how it usually done with Bruker data - opt.predefCoilAmpl = 0; % Use the predefined coil scaling coefficients, when available - opt.rmBadAvg = 0; - opt.doBlockAveraging = 0; - opt.doDriftCorrection = 0; - opt.iterin = 20; - opt.tmaxin = 0.2; - opt.aaDomain = 'f'; + opt.doECCbeforeAvg = 0; % 1 is how it usually done with Bruker data + opt.predefCoilAmpl = 0; % Use the predefined coil scaling coefficients, when available + opt.rmBadAvg = 1; + opt.doBlockAveraging = 0; + opt.doDriftCorrection = 1; + opt.iterin = 20; + opt.tmaxin = 0.2; + opt.aaDomain = 'f'; + opt.autophase = 1; + opt.compFracGroupDelay=1; end @@ -53,21 +55,27 @@ end [in, inw, inw_auto] = compMRS_DPload(DPid); + autoWaterExists=~isempty(inw_auto); %Loop through subjects and sessions. + + out = cell(check.nSubj); + outw = cell(check.nSubj); + out_auto = cell(check.nSubj); + outw_auto = cell(check.nSubj); + for m = 1:check.nSubj for n = 1:check.nSes(m) disp(['Processing ' DPid ' sub-' num2str(m) ' ses-' num2str(n)]) - + ident = [DPid '_sub-' num2str(m) '_ses-' num2str(n)]; % If separate water scan is available if ~isempty(inw) && ~isempty(inw{m, n}) - ident = [DPid '_sub-' num2str(m) '_ses-' num2str(n) '_sepWater']; - [out{m, n}, outw{m, n}] = compMRS_DPproc_sub(in{m, n}, inw{m, n}, ident, opt); + + [out{m}{n}, outw{m}{n}] = compMRS_DPproc_sub(in{m, n}, inw{m, n}, [ident '_sepWater'], opt); end % If automatic water scan is available - if ~isempty(inw_auto) && ~isempty(inw_auto{m, n}) - ident = [DPid '_sub-' num2str(m) '_ses-' num2str(n) '_autoWater']; - [out{m, n}, outw{m, n}] = compMRS_DPproc_sub(in{m, n}, inw_auto{m, n}, ident, opt); + if autoWaterExists && ~isempty(inw_auto{m, n}) + [out_auto{m}{n}, outw_auto{m}{n}] = compMRS_DPproc_sub(in{m, n}, inw_auto{m, n}, [ident '_autoWater'], opt); end end end @@ -86,11 +94,12 @@ out_mn = op_leftshift_keepSize(in_mn, floor(ls)); outw_mn= op_leftshift_keepSize(inw_mn, floor(ls)); - % Here check whether subspectra have to be combined first. + % Average the water + outw_mn=op_averaging(outw_mn); % do ECC before averaging (yes/no) if opt.doECCbeforeAvg - [out_mn, outw_mn]=op_eccKlose(out_mn,outw_mn); + [out_mn, outw_mn]=op_eccKlose(out_mn, outw_mn); end % Do coil combination WITHOUT averaging (if applicable) @@ -99,7 +108,7 @@ % The code below is mostly from op_combineRcvrs %first find the weights using the water unsuppressed data: - coilcombos_to_apply=op_getcoilcombos(outw_mn,2,'h'); + coilcombos_to_apply=op_getcoilcombos(outw_mn, 2,'h'); % If we choose to use the predefined weights then set the weights read from the method file (when available) if opt.predefCoilAmpl && isfield(out_mn, 'coilcombos') @@ -128,13 +137,17 @@ end end - % Process the water - outw_mn=op_alignAverages(outw_mn,0.2,'n'); - outw_mn=op_averaging(outw_mn); + + + % Combine subspectra for SPECIAL + + if strcmp(out_mn.seq, 'SPECIAL') && out_mn.dims.subSpecs>0 + out_mn=op_combinesubspecs(out_mn,"diff"); + end % do bad averages removal, if applicable (code from Jamie) if opt.rmBadAvg - outw_mn=subBadAveragesRemoval(outw_mn, opt); + out_mn=subBadAveragesRemoval(out_mn, ident, opt); end % Perform partial averaging, if applicable @@ -178,16 +191,17 @@ [out_part_avg, ~]=op_eccKlose(out_part_avg, outw_mn); end - compensateFractionalGroupDelay=1; - if compensateFractionalGroupDelay + + if opt.compFracGroupDelay && abs(frac_ls)>0.00001 ph1 = -frac_ls*in_mn.dwelltime; - out_part_avg=op_addphase(out_part_avg,0,ph1,4.7,1); + out_part_avg=op_addphase(out_part_avg, 0, ph1, out_part_avg.centerfreq, 1); end - % Compute the quality metrics % Get LW (of NAA) and SNR - out_part_avg = op_autophase(out_part_avg, 1.7, 2.3); + if opt.autophase + out_part_avg = op_autophase(out_part_avg, 2.8, 3.2); + end [FWHM_NAA] = op_getLW(out_part_avg, 1.8, 2.2, 8, 1); [SNR]=op_getSNR(out_part_avg,1.8,2.2,-2, 0, 1); @@ -220,6 +234,7 @@ for ii=1:length(out_all) plotlegend{ii} = [num2str(out_all{ii}.block_size) ' avg/block']; end + f=figure ('name', ident); hold on box on @@ -234,7 +249,7 @@ end legend(plotlegend) saveas(f, ['plots/' ident], 'fig') - print(f, ['plots/' ident], '-dpng', '-r600'); + print(f, ['plots/' ident], '-dpng', '-r300'); end diff --git a/code/compMRS_DPproc_allDPs.m b/code/compMRS_DPproc_allDPs.m index 46739df..9494520 100644 --- a/code/compMRS_DPproc_allDPs.m +++ b/code/compMRS_DPproc_allDPs.m @@ -23,7 +23,7 @@ % n is the number of sessions in the DP{k} % Each element {k}{m, n} is a water unsuppressed FID-A data struct. -function [out,outw]=compMRS_DPproc_allDPs() +function [out, outw, out_auto, outw_auto]=compMRS_DPproc_allDPs() % Look for all DPs in the current folder res = dir('DP*'); @@ -31,6 +31,6 @@ % run compMRS_DPproc on all DPs for ii=1:length(res) - [out{ii},outw{ii}]=compMRS_DPproc(res(ii).name); + [out{ii}, outw{ii}, out_auto{ii}, outw_auto{ii}]=compMRS_DPproc(res(ii).name); end end \ No newline at end of file diff --git a/code/compMRS_loadspecBruker.m b/code/compMRS_loadspecBruker.m index 050ce62..9feb220 100644 --- a/code/compMRS_loadspecBruker.m +++ b/code/compMRS_loadspecBruker.m @@ -102,6 +102,7 @@ headerRECO = compMRS_parseBrukerFormat(recoFile); end + % Get a few important bits % Software version version = headerACQP.ACQ_sw_version; @@ -228,8 +229,10 @@ if contains(version, ["PV 5"]) centerfreq = 4.7 + headerMethod.PVM_SpecOffsetppm; centerfreq_ref = centerfreq; -elseif contains(version,["PV 6", "PV 7", "PV-7" "PV-360"]) - centerfreq = 4.7; %headerMethod.PVM_FrqWorkPpm(1); +elseif contains(version,["PV 6", "PV 7", "PV-7"]) + centerfreq = 4.7 + headerACQP.ACQ_O1B_list/headerACQP.BF1; +elseif contains(version,["PV-360"]) + centerfreq = headerMethod.PVM_FrqWorkPpm(1); end % Receiver gain which we seem to need when reconciling raw and processed data @@ -342,7 +345,7 @@ end averages=rawAverages; %since these data are uncombined; - out.flags.averaged=0; %make the flags structure + %If there are multiple receivers *I think* that these always get stored %separately by default in the fid.raw file. Therefore, at this stage, @@ -352,7 +355,7 @@ % Reshape into a Npts x Ncoils x Naverages x Nrepetitions array fids_raw=reshape(fids_raw,rawDataPoints,Nrcvrs,rawAverages,rawRepetitions); %Permute so that time is along 1st dimension, averages is along 2nd - %dimension, and coils is along 3rd dimension: + %dimension, and coils is along 3rd dimension, repetitions along 4th dimension: fids_raw=permute(fids_raw,[1,3,2,4]); elseif ~contains(version,'PV 5') && ~multiRcvrs @@ -365,14 +368,10 @@ % Reshape into a Npts x Naverages x Nrepetitions array fids_raw=reshape(fids_raw,rawDataPoints,rawAverages,rawRepetitions); - - - end - % Remove singletons fids_raw = squeeze(fids_raw); @@ -453,10 +452,10 @@ ppm=f/(txfrq/1e6)+centerfreq; % Apply ref freq shift (the difference between txfrq and txfrq_ref) -if contains(version, ["PV-360"]) || (contains(version, ["PV 6", "PV-7"]) && strcmpi(rawData,'y')) - tmat=repmat(t',[1 sz(2:end)]); - fids=fids.*exp(-1i*tmat*(txfrq_ref-txfrq)*2*pi); -end +% if contains(version, ["PV-360"]) || (contains(version, ["PV 6", "PV-7"]) && strcmpi(rawData,'y')) +% tmat=repmat(t',[1 sz(2:end)]); +% fids=fids.*exp(-1i*tmat*(txfrq_ref-txfrq)*2*pi); +% end %Do the fourier transform specs=fftshift(ifft(fids,[],1),1); @@ -472,7 +471,7 @@ %Coils dimension should normally be after the averages dimension, %unless there are no averages, in which case the coils dimension will %be after the time dimension. - if rawAverages==1 && rawRepetitions==1 + if rawAverages==1 % && rawRepetitions==1 dims.coils=2; elseif rawAverages>1 dims.coils=3; @@ -524,10 +523,10 @@ % This will only work for PV 6 and above if (isfield(headerMethod, 'PVM_RefScanYN') && headerMethod.PVM_RefScanYN=='Yes' && isfield(headerMethod, 'PVM_RefScan')) isRef=true; - fids_ref=headerMethod.PVM_RefScan; - real_ref = fids_ref(1:2:length(fids_ref)); - imag_ref = fids_ref(2:2:length(fids_ref)); - fids_ref=complex(real_ref, imag_ref); + reffids=headerMethod.PVM_RefScan; + real_ref = reffids(1:2:length(reffids)); + imag_ref = reffids(2:2:length(reffids)); + reffids=complex(real_ref, imag_ref); % Find the number of averages that were acquired if isfield(headerMethod, 'PVM_RefScanNA') @@ -540,19 +539,19 @@ averages_ref=1; ref.flags.averaged=1; - reffids=reshape(fids_ref,[], Nrcvrs); + reffids=reshape(reffids,[], Nrcvrs); % Removed: we will perform leftshift in compMRS_DPproc - Thanh 20260101 - % fids_ref_trunc=fids_ref(leftshift+1:end,:,:); - % reffids=padarray(fids_ref_trunc, [leftshift,0],'post'); + % reffids_trunc=reffids(leftshift+1:end,:,:); + % reffids=padarray(reffids_trunc, [leftshift,0],'post'); % Apply ref freq shift (the difference between txfrq and txfrq_ref) % (but not if it's PV360) sz_ref=size(reffids); %size of the array - if ~contains(version, ["PV-360"]) - tmat=repmat(t',[1 sz_ref(2:end)]); - reffids=reffids.*exp(-1i*tmat*(txfrq_ref-txfrq)*2*pi); - end + % if ~contains(version, ["PV-360"]) + % tmat=repmat(t',[1 sz_ref(2:end)]); + % reffids=reffids.*exp(-1i*tmat*(txfrq_ref-txfrq)*2*pi); + % end refspecs=fftshift(ifft(reffids,[],1),1); @@ -570,7 +569,7 @@ refdims.extras=0; %Specify the number of subspecs. For ref data, I think these are always - %1 + %1 (as the data are already averaged). subspecs_ref=1; rawSubspecs_ref=1; end @@ -591,13 +590,13 @@ if isRef if contains(version, ["PV 5"]) fileRef = fullfile(inDir, 'fid.refscan'); - fids_ref = compMRS_readBrukerRaw(fileRef, 'int'); % Note: never tested for PV5 as we don't have datasets available + reffids = compMRS_readBrukerRaw(fileRef, 'int'); % Note: never tested for PV5 as we don't have datasets available elseif contains(version, ["PV 6", "PV 7", "PV-7", "PV-360.1", "PV-360.2"]) fileRef = fullfile(inDir, 'fid.refscan'); - fids_ref = compMRS_readBrukerRaw(fileRef, 'int32'); + reffids = compMRS_readBrukerRaw(fileRef, 'int32'); elseif contains(version,'PV-360.3') fileRef = fullfile(inDir, 'pdata', '1', 'fid_refscan.64'); - fids_ref = compMRS_readBrukerRaw(fileRef, 'float64'); + reffids = compMRS_readBrukerRaw(fileRef, 'float64'); end % The reference scan stored will be already combined and averaged @@ -610,18 +609,18 @@ rawAverages_ref=1; end - reffids=reshape(fids_ref,[],averages_ref); + reffids=reshape(reffids,[],averages_ref); % Removed: we will perform leftshift in compMRS_DPproc - Thanh 20260101 - % fids_ref_trunc=fids_ref(leftshift+1:end,:); - % reffids=padarray(fids_ref_trunc, [leftshift,0],'post'); + % reffids_trunc=reffids(leftshift+1:end,:); + % reffids=padarray(reffids_trunc, [leftshift,0],'post'); % % Apply ref freq shift (the difference between txfrq and txfrq_ref) % (but not if it's PV360) sz_ref=size(reffids); %size of the array - if ~contains(version, ["PV-360"]) - tmat=repmat(t',[1 sz_ref(2:end)]); - reffids=reffids.*exp(-1i*tmat*(txfrq_ref-txfrq)*2*pi); - end + % if ~contains(version, ["PV-360"]) + % tmat=repmat(t',[1 sz_ref(2:end)]); + % reffids=reffids.*exp(-1i*tmat*(txfrq_ref-txfrq)*2*pi); + % end refspecs=fftshift(ifft(reffids,[],1),1); @@ -670,10 +669,12 @@ isNav=true; data = fopen(fullfile(inDir, 'rawdata.job1')); fids_nav=fread(data,'int32'); + fclose(data); elseif exist(fullfile(inDir, 'rawdata.Navigator')); isNav=true; data = fopen(fullfile(inDir, 'rawdata.Navigator')); fids_nav=fread(data,'int32'); + fclose(data); end % From the PV 6 manual: rawdata.job1: Contains serially stored FIDs of @@ -692,7 +693,7 @@ end %Find the number of averages in the ref dataset: - rawAverages_nav=rawAverages; + rawAverages_nav=rawAverages*rawRepetitions; averages_nav=rawAverages_nav; if contains(version,'PV 5') @@ -731,21 +732,14 @@ ppm_nav=f/(txfrq/1e6)+centerfreq; % Apply ref freq shift (the difference between txfrq and txfrq_ref) - if contains(version, ["PV-360"]) || (contains(version, ["PV 6", "PV-7"]) && strcmpi(rawData,'y')) - tmat_nav=repmat(t_nav',[1 length_nav(2:end)]); - navfids=navfids.*exp(-1i*tmat_nav*(txfrq_ref-txfrq)*2*pi); - end + % if ~contains(version, ["PV-360"]) || (contains(version, ["PV 6", "PV-7"]) && strcmpi(rawData,'y')) + % tmat_nav=repmat(t_nav',[1 length_nav(2:end)]); + % navfids=navfids.*exp(-1i*tmat_nav*(txfrq_ref-txfrq)*2*pi); + % end %Do the fourier transform specs_nav=fftshift(ifft(navfids,[],1),1); - % Apply ref freq shift (the difference between txfrq and txfrq_ref) - % (but not if it's PV360) - % if ~contains(version, ["PV-360"]) - % tmat=repmat(t',[1 sz(2:end)]); - % navfids=navfids.*exp(-1i*tmat*(txfrq_ref-txfrq)*2*pi); - % end - nav.flags.averaged=0; %specify the dims navdims.t=1; @@ -764,6 +758,38 @@ disp('WARNING NAVIGATOR SCANS NOT FOUND. RETURNING EMPTY NAV STRUCTURE.'); end + +% Handle data where repetitions are used instead of averages, we'll need to +% swap both dimensions +if rawAverages==1 && rawRepetitions > 1 + dims.averages = dims.subSpecs; + dims.subSpecs = 0; + rawAverages = rawRepetitions; + rawRepetitions = 1; + + % handle subspecs parameter + % if length(size(fids_raw))==3 + % fids_raw = permute(fids_raw, [1 3 2]); + % elseif length(size(fids_raw))==4 + % fids_raw = permute(fids_raw, [1 4 3 2]); + % end + +end + +% Handle SPECIAL acquisitions; split the ISIS on/off into two subspecs and averages. +if strcmp(sequence, 'SPECIAL') && (dims.subSpecs==0) % Here we checked that it was not already split into subspecs. + % We need to arrange the data into subspecs + + fids = cat(4, fids(:, :, 1:2:end,: ), fids(:, :, 2:2:end, :)); + specs=fftshift(ifft(fids,[],1),1); + rawRepetitions=2; + averages = size(fids, 3); + dims.subSpecs = 4; + sz = size(fids); +end + + +out.flags.averaged=(rawAverages==1); %make the flags structure %FILLING IN DATA STRUCTURE FOR THE RAW DATA out.fids=fids; out.specs=specs; @@ -776,6 +802,7 @@ out.date=date; out.dims=dims; out.Bo=Bo; +out.centerfreq=centerfreq; out.averages=averages; out.rawAverages=rawAverages; out.subspecs=subspecs; @@ -968,6 +995,7 @@ out.coilcombos=coilcombos; out.isECCed=isECCed; out.isRFLed=isRFLed; + end % THE FUNCTION BELOW WAS COMMENTED OUT AS A SEPARATE FUNCTION FILE NAMED From b00842318eef50a4bbb0adcce8c152e827c17aee Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Thanh=20Phong=20L=C3=AA?= <43346967+tpkle@users.noreply.github.com> Date: Fri, 16 Jan 2026 16:35:14 +0100 Subject: [PATCH 8/8] Adding support for SPECIAL and correct frequency shifts There are still few phasing issues. Most DPs can be phased using op_autophase at the end, but for some datasets (DP 06, 07, 08, 36), the result is better without autophase. There are several issues with DP17, and the SNR seems very poor. --- code/compMRS_DPproc.m | 71 ++++++++++++-------- code/compMRS_DPproc_allDPs.m | 4 +- code/compMRS_loadspecBruker.m | 122 +++++++++++++++++++++------------- 3 files changed, 120 insertions(+), 77 deletions(-) diff --git a/code/compMRS_DPproc.m b/code/compMRS_DPproc.m index a51a27e..747bac4 100644 --- a/code/compMRS_DPproc.m +++ b/code/compMRS_DPproc.m @@ -26,19 +26,21 @@ % Both outputs are spectrally processed. out{m, n} contains the linewidth and SNR. -function [out, outw] = compMRS_DPproc(DPid, opt) +function [out, outw, out_auto, outw_auto] = compMRS_DPproc(DPid, opt) disp(['Processing ' num2str(DPid)]) if ~exist('opt','var') % Default options for the processing - opt.doECCbeforeAvg = 0; % 1 is how it usually done with Bruker data - opt.predefCoilAmpl = 0; % Use the predefined coil scaling coefficients, when available - opt.rmBadAvg = 0; - opt.doBlockAveraging = 0; - opt.doDriftCorrection = 0; - opt.iterin = 20; - opt.tmaxin = 0.2; - opt.aaDomain = 'f'; + opt.doECCbeforeAvg = 0; % 1 is how it usually done with Bruker data + opt.predefCoilAmpl = 0; % Use the predefined coil scaling coefficients, when available + opt.rmBadAvg = 1; + opt.doBlockAveraging = 0; + opt.doDriftCorrection = 1; + opt.iterin = 20; + opt.tmaxin = 0.2; + opt.aaDomain = 'f'; + opt.autophase = 1; + opt.compFracGroupDelay=1; end @@ -53,21 +55,27 @@ end [in, inw, inw_auto] = compMRS_DPload(DPid); + autoWaterExists=~isempty(inw_auto); %Loop through subjects and sessions. + + out = cell(check.nSubj); + outw = cell(check.nSubj); + out_auto = cell(check.nSubj); + outw_auto = cell(check.nSubj); + for m = 1:check.nSubj for n = 1:check.nSes(m) disp(['Processing ' DPid ' sub-' num2str(m) ' ses-' num2str(n)]) - + ident = [DPid '_sub-' num2str(m) '_ses-' num2str(n)]; % If separate water scan is available if ~isempty(inw) && ~isempty(inw{m, n}) - ident = [DPid '_sub-' num2str(m) '_ses-' num2str(n) '_sepWater']; - [out{m, n}, outw{m, n}] = compMRS_DPproc_sub(in{m, n}, inw{m, n}, ident, opt); + + [out{m}{n}, outw{m}{n}] = compMRS_DPproc_sub(in{m, n}, inw{m, n}, [ident '_sepWater'], opt); end % If automatic water scan is available - if ~isempty(inw_auto) && ~isempty(inw_auto{m, n}) - ident = [DPid '_sub-' num2str(m) '_ses-' num2str(n) '_autoWater']; - [out{m, n}, outw{m, n}] = compMRS_DPproc_sub(in{m, n}, inw_auto{m, n}, ident, opt); + if autoWaterExists && ~isempty(inw_auto{m, n}) + [out_auto{m}{n}, outw_auto{m}{n}] = compMRS_DPproc_sub(in{m, n}, inw_auto{m, n}, [ident '_autoWater'], opt); end end end @@ -86,11 +94,12 @@ out_mn = op_leftshift_keepSize(in_mn, floor(ls)); outw_mn= op_leftshift_keepSize(inw_mn, floor(ls)); - % Here check whether subspectra have to be combined first. + % Average the water + outw_mn=op_averaging(outw_mn); % do ECC before averaging (yes/no) if opt.doECCbeforeAvg - [out_mn, outw_mn]=op_eccKlose(out_mn,outw_mn); + [out_mn, outw_mn]=op_eccKlose(out_mn, outw_mn); end % Do coil combination WITHOUT averaging (if applicable) @@ -99,7 +108,7 @@ % The code below is mostly from op_combineRcvrs %first find the weights using the water unsuppressed data: - coilcombos_to_apply=op_getcoilcombos(outw_mn,2,'h'); + coilcombos_to_apply=op_getcoilcombos(outw_mn, 2,'h'); % If we choose to use the predefined weights then set the weights read from the method file (when available) if opt.predefCoilAmpl && isfield(out_mn, 'coilcombos') @@ -128,13 +137,17 @@ end end - % Process the water - outw_mn=op_alignAverages(outw_mn,0.2,'n'); - outw_mn=op_averaging(outw_mn); + + + % Combine subspectra for SPECIAL + + if strcmp(out_mn.seq, 'SPECIAL') && out_mn.dims.subSpecs>0 + out_mn=op_combinesubspecs(out_mn,"diff"); + end % do bad averages removal, if applicable (code from Jamie) if opt.rmBadAvg - outw_mn=subBadAveragesRemoval(outw_mn, opt); + out_mn=subBadAveragesRemoval(out_mn, ident, opt); end % Perform partial averaging, if applicable @@ -178,16 +191,17 @@ [out_part_avg, ~]=op_eccKlose(out_part_avg, outw_mn); end - compensateFractionalGroupDelay=1; - if compensateFractionalGroupDelay + + if opt.compFracGroupDelay && abs(frac_ls)>0.00001 ph1 = -frac_ls*in_mn.dwelltime; - out_part_avg=op_addphase(out_part_avg,0,ph1,4.7,1); + out_part_avg=op_addphase(out_part_avg, 0, ph1, out_part_avg.centerfreq, 1); end - % Compute the quality metrics % Get LW (of NAA) and SNR - out_part_avg = op_autophase(out_part_avg, 1.7, 2.3); + if opt.autophase + out_part_avg = op_autophase(out_part_avg, 2.8, 3.2); + end [FWHM_NAA] = op_getLW(out_part_avg, 1.8, 2.2, 8, 1); [SNR]=op_getSNR(out_part_avg,1.8,2.2,-2, 0, 1); @@ -220,6 +234,7 @@ for ii=1:length(out_all) plotlegend{ii} = [num2str(out_all{ii}.block_size) ' avg/block']; end + f=figure ('name', ident); hold on box on @@ -234,7 +249,7 @@ end legend(plotlegend) saveas(f, ['plots/' ident], 'fig') - print(f, ['plots/' ident], '-dpng', '-r600'); + print(f, ['plots/' ident], '-dpng', '-r300'); end diff --git a/code/compMRS_DPproc_allDPs.m b/code/compMRS_DPproc_allDPs.m index 46739df..9494520 100644 --- a/code/compMRS_DPproc_allDPs.m +++ b/code/compMRS_DPproc_allDPs.m @@ -23,7 +23,7 @@ % n is the number of sessions in the DP{k} % Each element {k}{m, n} is a water unsuppressed FID-A data struct. -function [out,outw]=compMRS_DPproc_allDPs() +function [out, outw, out_auto, outw_auto]=compMRS_DPproc_allDPs() % Look for all DPs in the current folder res = dir('DP*'); @@ -31,6 +31,6 @@ % run compMRS_DPproc on all DPs for ii=1:length(res) - [out{ii},outw{ii}]=compMRS_DPproc(res(ii).name); + [out{ii}, outw{ii}, out_auto{ii}, outw_auto{ii}]=compMRS_DPproc(res(ii).name); end end \ No newline at end of file diff --git a/code/compMRS_loadspecBruker.m b/code/compMRS_loadspecBruker.m index 050ce62..9feb220 100644 --- a/code/compMRS_loadspecBruker.m +++ b/code/compMRS_loadspecBruker.m @@ -102,6 +102,7 @@ headerRECO = compMRS_parseBrukerFormat(recoFile); end + % Get a few important bits % Software version version = headerACQP.ACQ_sw_version; @@ -228,8 +229,10 @@ if contains(version, ["PV 5"]) centerfreq = 4.7 + headerMethod.PVM_SpecOffsetppm; centerfreq_ref = centerfreq; -elseif contains(version,["PV 6", "PV 7", "PV-7" "PV-360"]) - centerfreq = 4.7; %headerMethod.PVM_FrqWorkPpm(1); +elseif contains(version,["PV 6", "PV 7", "PV-7"]) + centerfreq = 4.7 + headerACQP.ACQ_O1B_list/headerACQP.BF1; +elseif contains(version,["PV-360"]) + centerfreq = headerMethod.PVM_FrqWorkPpm(1); end % Receiver gain which we seem to need when reconciling raw and processed data @@ -342,7 +345,7 @@ end averages=rawAverages; %since these data are uncombined; - out.flags.averaged=0; %make the flags structure + %If there are multiple receivers *I think* that these always get stored %separately by default in the fid.raw file. Therefore, at this stage, @@ -352,7 +355,7 @@ % Reshape into a Npts x Ncoils x Naverages x Nrepetitions array fids_raw=reshape(fids_raw,rawDataPoints,Nrcvrs,rawAverages,rawRepetitions); %Permute so that time is along 1st dimension, averages is along 2nd - %dimension, and coils is along 3rd dimension: + %dimension, and coils is along 3rd dimension, repetitions along 4th dimension: fids_raw=permute(fids_raw,[1,3,2,4]); elseif ~contains(version,'PV 5') && ~multiRcvrs @@ -365,14 +368,10 @@ % Reshape into a Npts x Naverages x Nrepetitions array fids_raw=reshape(fids_raw,rawDataPoints,rawAverages,rawRepetitions); - - - end - % Remove singletons fids_raw = squeeze(fids_raw); @@ -453,10 +452,10 @@ ppm=f/(txfrq/1e6)+centerfreq; % Apply ref freq shift (the difference between txfrq and txfrq_ref) -if contains(version, ["PV-360"]) || (contains(version, ["PV 6", "PV-7"]) && strcmpi(rawData,'y')) - tmat=repmat(t',[1 sz(2:end)]); - fids=fids.*exp(-1i*tmat*(txfrq_ref-txfrq)*2*pi); -end +% if contains(version, ["PV-360"]) || (contains(version, ["PV 6", "PV-7"]) && strcmpi(rawData,'y')) +% tmat=repmat(t',[1 sz(2:end)]); +% fids=fids.*exp(-1i*tmat*(txfrq_ref-txfrq)*2*pi); +% end %Do the fourier transform specs=fftshift(ifft(fids,[],1),1); @@ -472,7 +471,7 @@ %Coils dimension should normally be after the averages dimension, %unless there are no averages, in which case the coils dimension will %be after the time dimension. - if rawAverages==1 && rawRepetitions==1 + if rawAverages==1 % && rawRepetitions==1 dims.coils=2; elseif rawAverages>1 dims.coils=3; @@ -524,10 +523,10 @@ % This will only work for PV 6 and above if (isfield(headerMethod, 'PVM_RefScanYN') && headerMethod.PVM_RefScanYN=='Yes' && isfield(headerMethod, 'PVM_RefScan')) isRef=true; - fids_ref=headerMethod.PVM_RefScan; - real_ref = fids_ref(1:2:length(fids_ref)); - imag_ref = fids_ref(2:2:length(fids_ref)); - fids_ref=complex(real_ref, imag_ref); + reffids=headerMethod.PVM_RefScan; + real_ref = reffids(1:2:length(reffids)); + imag_ref = reffids(2:2:length(reffids)); + reffids=complex(real_ref, imag_ref); % Find the number of averages that were acquired if isfield(headerMethod, 'PVM_RefScanNA') @@ -540,19 +539,19 @@ averages_ref=1; ref.flags.averaged=1; - reffids=reshape(fids_ref,[], Nrcvrs); + reffids=reshape(reffids,[], Nrcvrs); % Removed: we will perform leftshift in compMRS_DPproc - Thanh 20260101 - % fids_ref_trunc=fids_ref(leftshift+1:end,:,:); - % reffids=padarray(fids_ref_trunc, [leftshift,0],'post'); + % reffids_trunc=reffids(leftshift+1:end,:,:); + % reffids=padarray(reffids_trunc, [leftshift,0],'post'); % Apply ref freq shift (the difference between txfrq and txfrq_ref) % (but not if it's PV360) sz_ref=size(reffids); %size of the array - if ~contains(version, ["PV-360"]) - tmat=repmat(t',[1 sz_ref(2:end)]); - reffids=reffids.*exp(-1i*tmat*(txfrq_ref-txfrq)*2*pi); - end + % if ~contains(version, ["PV-360"]) + % tmat=repmat(t',[1 sz_ref(2:end)]); + % reffids=reffids.*exp(-1i*tmat*(txfrq_ref-txfrq)*2*pi); + % end refspecs=fftshift(ifft(reffids,[],1),1); @@ -570,7 +569,7 @@ refdims.extras=0; %Specify the number of subspecs. For ref data, I think these are always - %1 + %1 (as the data are already averaged). subspecs_ref=1; rawSubspecs_ref=1; end @@ -591,13 +590,13 @@ if isRef if contains(version, ["PV 5"]) fileRef = fullfile(inDir, 'fid.refscan'); - fids_ref = compMRS_readBrukerRaw(fileRef, 'int'); % Note: never tested for PV5 as we don't have datasets available + reffids = compMRS_readBrukerRaw(fileRef, 'int'); % Note: never tested for PV5 as we don't have datasets available elseif contains(version, ["PV 6", "PV 7", "PV-7", "PV-360.1", "PV-360.2"]) fileRef = fullfile(inDir, 'fid.refscan'); - fids_ref = compMRS_readBrukerRaw(fileRef, 'int32'); + reffids = compMRS_readBrukerRaw(fileRef, 'int32'); elseif contains(version,'PV-360.3') fileRef = fullfile(inDir, 'pdata', '1', 'fid_refscan.64'); - fids_ref = compMRS_readBrukerRaw(fileRef, 'float64'); + reffids = compMRS_readBrukerRaw(fileRef, 'float64'); end % The reference scan stored will be already combined and averaged @@ -610,18 +609,18 @@ rawAverages_ref=1; end - reffids=reshape(fids_ref,[],averages_ref); + reffids=reshape(reffids,[],averages_ref); % Removed: we will perform leftshift in compMRS_DPproc - Thanh 20260101 - % fids_ref_trunc=fids_ref(leftshift+1:end,:); - % reffids=padarray(fids_ref_trunc, [leftshift,0],'post'); + % reffids_trunc=reffids(leftshift+1:end,:); + % reffids=padarray(reffids_trunc, [leftshift,0],'post'); % % Apply ref freq shift (the difference between txfrq and txfrq_ref) % (but not if it's PV360) sz_ref=size(reffids); %size of the array - if ~contains(version, ["PV-360"]) - tmat=repmat(t',[1 sz_ref(2:end)]); - reffids=reffids.*exp(-1i*tmat*(txfrq_ref-txfrq)*2*pi); - end + % if ~contains(version, ["PV-360"]) + % tmat=repmat(t',[1 sz_ref(2:end)]); + % reffids=reffids.*exp(-1i*tmat*(txfrq_ref-txfrq)*2*pi); + % end refspecs=fftshift(ifft(reffids,[],1),1); @@ -670,10 +669,12 @@ isNav=true; data = fopen(fullfile(inDir, 'rawdata.job1')); fids_nav=fread(data,'int32'); + fclose(data); elseif exist(fullfile(inDir, 'rawdata.Navigator')); isNav=true; data = fopen(fullfile(inDir, 'rawdata.Navigator')); fids_nav=fread(data,'int32'); + fclose(data); end % From the PV 6 manual: rawdata.job1: Contains serially stored FIDs of @@ -692,7 +693,7 @@ end %Find the number of averages in the ref dataset: - rawAverages_nav=rawAverages; + rawAverages_nav=rawAverages*rawRepetitions; averages_nav=rawAverages_nav; if contains(version,'PV 5') @@ -731,21 +732,14 @@ ppm_nav=f/(txfrq/1e6)+centerfreq; % Apply ref freq shift (the difference between txfrq and txfrq_ref) - if contains(version, ["PV-360"]) || (contains(version, ["PV 6", "PV-7"]) && strcmpi(rawData,'y')) - tmat_nav=repmat(t_nav',[1 length_nav(2:end)]); - navfids=navfids.*exp(-1i*tmat_nav*(txfrq_ref-txfrq)*2*pi); - end + % if ~contains(version, ["PV-360"]) || (contains(version, ["PV 6", "PV-7"]) && strcmpi(rawData,'y')) + % tmat_nav=repmat(t_nav',[1 length_nav(2:end)]); + % navfids=navfids.*exp(-1i*tmat_nav*(txfrq_ref-txfrq)*2*pi); + % end %Do the fourier transform specs_nav=fftshift(ifft(navfids,[],1),1); - % Apply ref freq shift (the difference between txfrq and txfrq_ref) - % (but not if it's PV360) - % if ~contains(version, ["PV-360"]) - % tmat=repmat(t',[1 sz(2:end)]); - % navfids=navfids.*exp(-1i*tmat*(txfrq_ref-txfrq)*2*pi); - % end - nav.flags.averaged=0; %specify the dims navdims.t=1; @@ -764,6 +758,38 @@ disp('WARNING NAVIGATOR SCANS NOT FOUND. RETURNING EMPTY NAV STRUCTURE.'); end + +% Handle data where repetitions are used instead of averages, we'll need to +% swap both dimensions +if rawAverages==1 && rawRepetitions > 1 + dims.averages = dims.subSpecs; + dims.subSpecs = 0; + rawAverages = rawRepetitions; + rawRepetitions = 1; + + % handle subspecs parameter + % if length(size(fids_raw))==3 + % fids_raw = permute(fids_raw, [1 3 2]); + % elseif length(size(fids_raw))==4 + % fids_raw = permute(fids_raw, [1 4 3 2]); + % end + +end + +% Handle SPECIAL acquisitions; split the ISIS on/off into two subspecs and averages. +if strcmp(sequence, 'SPECIAL') && (dims.subSpecs==0) % Here we checked that it was not already split into subspecs. + % We need to arrange the data into subspecs + + fids = cat(4, fids(:, :, 1:2:end,: ), fids(:, :, 2:2:end, :)); + specs=fftshift(ifft(fids,[],1),1); + rawRepetitions=2; + averages = size(fids, 3); + dims.subSpecs = 4; + sz = size(fids); +end + + +out.flags.averaged=(rawAverages==1); %make the flags structure %FILLING IN DATA STRUCTURE FOR THE RAW DATA out.fids=fids; out.specs=specs; @@ -776,6 +802,7 @@ out.date=date; out.dims=dims; out.Bo=Bo; +out.centerfreq=centerfreq; out.averages=averages; out.rawAverages=rawAverages; out.subspecs=subspecs; @@ -968,6 +995,7 @@ out.coilcombos=coilcombos; out.isECCed=isECCed; out.isRFLed=isRFLed; + end % THE FUNCTION BELOW WAS COMMENTED OUT AS A SEPARATE FUNCTION FILE NAMED