From 4793f7065cb9ca224388bf61580ff8b4f1284c7b Mon Sep 17 00:00:00 2001
From: Helge <47348963+HJZollner@users.noreply.github.com>
Date: Wed, 23 Mar 2022 17:55:24 -0400
Subject: [PATCH 1/4] Update SpecReg

Update dimensions for all spectral registration functions
---
 .../processingTools/op_SpecRegFreqRestrict.m  |  15 +-
 .../FID-A/processingTools/op_probabSpecReg.m  | 357 ++++++++++++++++++
 .../FID-A/processingTools/op_robustSpecReg.m  |   5 +-
 3 files changed, 372 insertions(+), 5 deletions(-)
 create mode 100644 libraries/FID-A/processingTools/op_probabSpecReg.m

diff --git a/libraries/FID-A/processingTools/op_SpecRegFreqRestrict.m b/libraries/FID-A/processingTools/op_SpecRegFreqRestrict.m
index 003c2f9d..731c0a47 100644
--- a/libraries/FID-A/processingTools/op_SpecRegFreqRestrict.m
+++ b/libraries/FID-A/processingTools/op_SpecRegFreqRestrict.m
@@ -41,6 +41,12 @@
     end
 end
 
+if strcmp(seqType,'HERMES') || strcmp(seqType,'HERMES')
+    [in, ~] = osp_onOffClassifyHERMES(in);
+    ppmMin = [0.5,0.5,1.85,1.85];
+    ppmMax = [3.9,3.8,3.9,3.8];
+end
+
 % Check whether data is coil-combined. If not, throw error.
 if ~in.flags.addedrcvrs
     error('ERROR:  I think it only makes sense to do this after you have combined the channels using op_addrcvrs.  ABORTING!!');
@@ -86,11 +92,12 @@
 % Initialize common variables and loop over sub-spectra
 t                 = in.t;
 input.dwelltime   = in.dwelltime;
+
+for mm=1:numSubSpecs
+
 in_restrict = op_freqrange(in,ppmMin,ppmMax);
 t_restrict                 = in_restrict.t;
-freq        = in_restrict.ppm;
-for mm=1:numSubSpecs
-       
+freq        = in_restrict.ppm;       
     
         DataToAlign = in_restrict.fids(:,:,mm);
         
@@ -261,7 +268,7 @@
 
 % Perform weighted averaging
 % No need for 'mean', since the weights vector w is normalized
-fids = sum(fids_out, in.dims.averages);
+fids = squeeze(sum(fids_out, in.dims.averages));
     
 %%% 3. WRITE THE NEW STRUCTURE %%%
 %re-calculate Specs using fft
diff --git a/libraries/FID-A/processingTools/op_probabSpecReg.m b/libraries/FID-A/processingTools/op_probabSpecReg.m
new file mode 100644
index 00000000..1329bba1
--- /dev/null
+++ b/libraries/FID-A/processingTools/op_probabSpecReg.m
@@ -0,0 +1,357 @@
+function [out, fs, phs, w, driftPre, driftPost] = op_probabSpecReg(in, seqType, echo, F0, noAlign)
+% Spectral registration is a time-domain frequency-and-phase correction
+% routine as per Near et al. (2015). Incorporates a multiplexed,
+% probabilistic approach for aligning HERMES data (MM: 170609)
+showPlots = 0;
+if nargin <5
+    noAlign = 0;
+    if nargin <4
+        F0 = nan;
+        if nargin < 2
+            echo = 1;    
+        end
+    end
+end
+
+% Check whether data is coil-combined. If not, throw error.
+if ~in.flags.addedrcvrs
+    error('ERROR:  I think it only makes sense to do this after you have combined the channels using op_addrcvrs.  ABORTING!!');
+end
+
+% Check whether data is averaged. If it is, throw warning and return unmodified data.
+if in.dims.averages == 0 || in.averages < 2
+    %DO NOTHING
+    warning('No averages found.  Returning input without modification!');
+    out = in;
+    out.flags.freqcorrected = 1;
+    fs  = 0;
+    phs = 0;
+    return
+end
+
+% Check whether input structure has sub-spectra.
+% If it does, prepare to loop over all sub-spectra. Each sub-spectrum will
+% get aligned separately.
+if in.dims.subSpecs == 0
+    numSubSpecs = 1;
+    % Measure drift pre-alignment
+    driftPre = op_measureDrift(in);
+else
+    numSubSpecs = in.sz(in.dims.subSpecs);
+    % Measure drift pre-alignment
+    for mm = 1:numSubSpecs
+        in_sub = op_takesubspec(in, mm);
+        driftPre{mm} = op_measureDrift(in_sub);
+    end
+end
+
+
+% Pre-allocate memory.
+fids    = zeros(in.sz(in.dims.t),1,numSubSpecs);
+fs      = zeros(in.sz(in.dims.averages),numSubSpecs);
+phs     = zeros(in.sz(in.dims.averages),numSubSpecs);
+zMSE = zeros(numSubSpecs,in.sz(in.dims.averages));
+CorrParsML = zeros(in.sz(in.dims.averages),2);
+parsGuess = [0 0];
+
+% Initialize common variables and loop over sub-spectra
+t                 = in.t;
+input.dwelltime   = in.dwelltime;
+
+% Probability density function and parameter bounds
+Cauchy = @(x,s,l) s./(pi.*(s.^2+(x-l).^2));
+lb = [0 -Inf];
+ub = [Inf Inf];
+
+% Optimization options
+lsqnonlinopts = optimoptions(@lsqnonlin);
+lsqnonlinopts = optimoptions(lsqnonlinopts,'Display','off','Algorithm','levenberg-marquardt');
+% nlinopts = statset('nlinfit');
+% nlinopts = statset(nlinopts,'MaxIter',400,'TolX',1e-8,'TolFun',1e-8);
+mleopts = statset('mlecustom');
+mleopts = statset(mleopts,'MaxIter',400,'MaxFunEvals',800,'TolX',1e-6,'TolFun',1e-6,'TolBnd',1e-6);
+
+% Set dimensions of figures of histograms
+if showPlots == 1
+    d.w = 0.6;
+    d.h = 0.45;
+    d.l = (1-d.w)/2;
+    d.b = (1-d.h)/2;
+end
+
+for mm=1:numSubSpecs
+    clear m;
+    freq        = in.ppm;
+    
+    DataToAlign = in.fids(:,:,mm);
+    
+    
+    % Use first n points of time-domain data, where n is the last point where SNR > 3
+    signal = abs(DataToAlign(:,mm));
+    noise = 2*std(signal(ceil(0.75*size(signal,1)):end,:));
+    SNR = signal ./ repmat(noise, [size(DataToAlign,1) 1]);
+    SNR = abs(diff(mean(SNR,2)));
+    
+    %SNR = abs(diff(mean(SNR,2))); This is how it's done in RobSpecReg
+    
+%     n = find(SNR > 1.5); % This is original Gannet
+%     if isempty(n)
+%         tMax = 100;
+%     else
+%         tMax = n(end);
+%         if tMax < 100
+%             tMax = 100;
+%         end
+%     end
+
+    SNR = SNR(t <= 0.2);
+    tMax = find(SNR > 1.5,1,'last');
+    if isempty(tMax) || tMax < find(t <= 0.1,1,'last')
+        tMax = find(t <= 0.1,1,'last');
+    end
+    
+    % 'Flatten' complex data for use in nlinfit
+    clear flatdata;
+    flatdata(:,1,:) = real(DataToAlign(1:tMax,:));
+    flatdata(:,2,:) = imag(DataToAlign(1:tMax,:));
+    
+    % Reference transient
+    flattarget = median(flatdata,3); % median across transients
+    target = flattarget(:);
+    
+    % Scalar to normalize transients (reduces optimization time)
+    a = max(abs(target));
+    
+    % Pre-allocate memory
+    if mm == 1
+        params = zeros(size(flatdata,3), 2);
+        MSE = zeros(1, size(flatdata,3));
+    end
+    
+     % Frame-by-frame determination of frequency of residual water and Cr (if HERMES or GSH editing)
+    if strcmp(seqType, 'HERMES') || strcmp(seqType, 'HERCULES')
+        F0freqRange = freq - 3.02 >= -0.15 & freq - 3.02 <= 0.15;
+    else
+        F0freqRange = freq - 3.22 >= -0.15 & freq - 3.22 <= 0.15;
+    end
+    [~,FrameMaxPos] = max(abs(real(in.specs(F0freqRange,:))),[],1);
+    F0freqRange = freq(F0freqRange);
+    F0freq = F0freqRange(FrameMaxPos);
+
+    % Starting values for optimization
+    if isnan(F0) 
+        f0 = F0freq * in.txfrq * 1e-6;
+        f0 = f0 - f0(1);
+    else
+        f0 = F0;
+        f0 = f0 - f0(1);
+    end
+    phi0 = zeros(size(f0));
+    x0 = [f0(:) phi0(:)];
+        
+    
+    % Determine frequency and phase offsets by spectral registration
+     if noAlign == 0
+        time = 0:input.dwelltime:(length(target)/2 - 1)*input.dwelltime;
+        reverseStr = '';
+        for jj = 1:size(flatdata,3)
+            if echo
+                msg = sprintf('\nSpectral registration - Fitting transient: %d', corrloop);
+                fprintf([reverseStr, msg]);
+                reverseStr = repmat(sprintf('\b'), 1, length(msg));
+            end
+
+            transient = squeeze(flatdata(:,:,jj));
+
+            fun = @(x) SpecReg(transient(:)/a, target/a, time, x);
+            params(jj,:) = lsqnonlin(fun, x0(jj,:), [], [], lsqnonlinopts);
+
+            f   = params(jj,1);
+            phi = params(jj,2);
+            m_c = complex(flatdata(:,1,jj), flatdata(:,2,jj));
+            m_c = m_c .* exp(1i*pi*(time'*f*2+phi/180));
+            m(:,jj) = [real(m_c); imag(m_c)];
+            resid = target - m(:,jj);
+            MSE(jj) = sum(resid.^2) / (length(resid) - 2);
+
+%             x0(jj,2) = phi; % Update phase value according to the last iteration
+
+    %         [parsFit(corrloop,:), ~, ~, ~, MSE(corrloop)] = nlinfit(input, target, @FreqPhaseShiftNest, parsGuess, nlinopts);
+    %         parsGuess = parsFit(corrloop,:);
+        end
+         % Prepare the frequency and phase corrections to be returned by this
+        % function for further use
+        fs(:,mm)  = params(:,1);
+        phs(:,mm) = params(:,2);
+    end
+    
+    
+    % Probability distribution of frequency offsets (estimated by maximum likelihood)
+    start = [iqr(fs(:,mm))/2, median(fs(:,mm))];
+    [fs_p(:,mm), fs_p_ci(:,:,mm)] = ...
+        mle(fs(:,mm), 'pdf', Cauchy, 'start', start, 'lower', lb, 'upper', ub, 'options', mleopts);
+    fs_fx(:,mm) = ...
+        linspace(1.5*min(fs(:,mm)), 1.5*max(fs(:,mm)), 1e3);
+    fs_pdf(:,mm) = Cauchy(fs_fx(:,mm), fs_p(1,mm), fs_p(2,mm));
+    
+    % Probability distribution of phase offsets (estimated by maximum likelihood)
+    start = [iqr(phs(:,mm))/2, median(phs(:,mm))];
+    [phs_p(:,mm), phs_p_ci(:,:,mm)] = ...
+        mle(phs(:,mm), 'pdf', Cauchy, 'start', start, 'lower', lb, 'upper', ub, 'options', mleopts);
+    phs_fx(:,mm) = ...
+        linspace(1.5*min(phs(:,mm)), 1.5*max(phs(:,mm)), 1e3);
+    phs_pdf(:,mm) = Cauchy(phs_fx(:,mm), phs_p(1,mm), phs_p(2,mm));
+             
+
+    zMSE = zscore(MSE); % standardized MSEs
+    
+    % Apply frequency and phase corrections
+     for jj = 1:size(flatdata,3)
+        % Default correction
+%         fids(:,jj,mm) = in.fids(:,jj,mm) .* ...
+%             exp(1i*params(jj,1)*2*pi*t') * exp(1i*pi/180*params(jj,2));
+
+        % Freq/phase correction + Cauchy pdf location parameter shift
+        fids(:,jj,mm) = in.fids(:,jj,mm) .* ...
+            exp(1i*(params(jj,1)-fs_p(2))*2*pi*t') * exp(1i*pi/180*(params(jj,2)-phs_p(2)));
+     end
+
+           
+        
+end
+
+if echo
+fprintf('\n');
+end
+
+%%% 2. CALCULATE THE DRIFT AFTER CORRECTIONS
+out_temp=in;
+out_temp.fids=fids;
+%re-calculate Specs using fft
+out_temp.specs=fftshift(fft(fids,[],in.dims.t),in.dims.t);
+% Measure drift post-alignment
+if in.dims.subSpecs == 0
+    numSubSpecs = 1;
+    driftPost = op_measureDrift(out_temp);
+else
+    numSubSpecs = in.sz(in.dims.subSpecs);
+    % Measure drift post-alignment
+    for mm = 1:numSubSpecs
+        out_temp_sub = op_takesubspec(out_temp, mm);
+        driftPost{mm} = op_measureDrift(out_temp_sub);
+    end
+end
+
+
+%%% 3. RE-RUN THE RANKING ALGORITHM TO DETERMINE WEIGHTS %%%
+
+% Determine optimal iteration order by calculating a similarity metric (mean squared error)
+w = cell(1,numSubSpecs);
+
+for mm=1:numSubSpecs
+    DataToWeight = fids(:,:,mm);
+    D = zeros(size(DataToWeight,2));
+    for jj = 1:size(DataToWeight,2)
+        for kk = 1:size(DataToWeight,2)
+            tmp = sum((real(DataToWeight(1:tMax,jj)) - real(DataToWeight(1:tMax,kk))).^2) / tMax;
+            if tmp == 0
+                D(jj,kk) = NaN;
+            else
+                D(jj,kk) = tmp;
+            end
+        end
+    end
+    d = nanmedian(D);
+    if isnan(sum(d))        
+        d(isnan(d)) = 1;
+    end
+    w{mm} = 1./d.^2;
+    w{mm} = w{mm}/sum(w{mm});
+    w{mm} = repmat(w{mm}, [size(DataToWeight,1) 1]);
+    
+    
+    
+    % Apply the weighting
+fids_out(:,:,mm) = w{mm} .* DataToWeight;
+    
+end
+
+% Or remove based on zMSE score
+% for mm=1:numSubSpecs
+%     DataToWeight = fids(:,:,mm);
+%     w{mm} = zMSE < 3;
+%     w{mm} = w{mm}/sum(w{mm});
+%     w{mm} = repmat(w{mm}, [size(DataToWeight,1) 1]);
+%         
+%     % Apply the weighting
+%     fids_out(:,:,mm) = w{mm} .* DataToWeight;
+%     
+% end
+
+% Perform weighted averaging
+% No need for 'mean', since the weights vector w is normalized
+fids = squeeze(sum(fids_out, in.dims.averages));
+    
+%%% 3. WRITE THE NEW STRUCTURE %%%
+%re-calculate Specs using fft
+specs=fftshift(fft(fids,[],in.dims.t),in.dims.t);
+
+%change the dims variables.
+if in.dims.t>in.dims.averages
+    dims.t=in.dims.t-1;
+else
+    dims.t=in.dims.t;
+end
+if in.dims.coils>in.dims.averages
+    dims.coils=in.dims.coils-1;
+else
+    dims.coils=in.dims.coils;
+end
+dims.averages=0;
+if in.dims.subSpecs>in.dims.averages
+    dims.subSpecs=in.dims.subSpecs-1;
+else
+    dims.subSpecs=in.dims.subSpecs;
+end
+if in.dims.extras>in.dims.averages
+    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.averages=1;
+    
+%FILLING IN THE FLAGS
+out.flags=in.flags;
+out.flags.writtentostruct=1;
+out.flags.freqcorrected=1;
+out.flags.averaged=1;
+
+
+
+end
+
+function out = SpecReg(data, target, time, x)
+
+f   = x(1);
+phi = x(2);
+
+data = reshape(data, [length(data)/2 2]);
+fid = complex(data(:,1), data(:,2));
+
+fid = fid .* exp(1i*pi*(time'*f*2+phi/180));
+fid = [real(fid); imag(fid)];
+
+out = target - fid;
+
+end
+
diff --git a/libraries/FID-A/processingTools/op_robustSpecReg.m b/libraries/FID-A/processingTools/op_robustSpecReg.m
index 76ca7182..e23fe1e0 100644
--- a/libraries/FID-A/processingTools/op_robustSpecReg.m
+++ b/libraries/FID-A/processingTools/op_robustSpecReg.m
@@ -51,6 +51,9 @@
     out.flags.freqcorrected = 1;
     fs  = 0;
     phs = 0;
+    w = 0;
+    driftPre = 0;
+    driftPost = 0;
     return
 end
 
@@ -321,7 +324,7 @@
 
 % Perform weighted averaging
 % No need for 'mean', since the weights vector w is normalized
-fids = sum(fids_out, in.dims.averages);
+fids = squeeze(sum(fids_out, in.dims.averages));
     
 %%% 3. WRITE THE NEW STRUCTURE %%%
 %re-calculate Specs using fft

From 6724ea9f6a30238b7d6d09c0031360c0e780b74e Mon Sep 17 00:00:00 2001
From: Helge <47348963+HJZollner@users.noreply.github.com>
Date: Wed, 23 Mar 2022 18:00:30 -0400
Subject: [PATCH 2/4] Revert "MEGA-MM"

This reverts commit f52b57bb72836f59390a77476a72533fabd454c9.
---
 fit/osp_fitMEGA.m                             |  34 -----
 fit/osp_fitUnEdited.m                         |  50 +++----
 .../Osprey/fit_OspreyParamsToModel.m          |  16 +--
 .../Osprey/fit_OspreyPrelimStep2MM.m          |   8 +-
 .../Osprey/fit_Osprey_PrelimReducedMM.m       |   8 +-
 .../FID-A/fitTools/fit_createMetabListMM.m    | 123 ++++++------------
 libraries/FID-A/fitTools/fit_makeBasis.m      |   6 +-
 libraries/FID-A/fitTools/fit_runFitMM.m       |  15 +--
 plot/osp_plotFit.m                            |   7 +-
 9 files changed, 88 insertions(+), 179 deletions(-)

diff --git a/fit/osp_fitMEGA.m b/fit/osp_fitMEGA.m
index f188fcac..419f4016 100644
--- a/fit/osp_fitMEGA.m
+++ b/fit/osp_fitMEGA.m
@@ -83,40 +83,6 @@
             % Save back the basis set and fit parameters to MRSCont
             MRSCont.fit.resBasisSet.diff1{kk}           = resBasisSetDiff1;
             MRSCont.fit.results.diff1.fitParams{kk} = fitParamsDiff1;
-            
-            %Modeling MM spectra after the main spectrum 
-            if MRSCont.flags.hasMM == 1            
-                dataToFit_mm   = MRSCont.processed.diff1_mm{kk};
-                dataToFit_mm   = op_ampScale(dataToFit_mm, 1/MRSCont.fit.scale{kk});
-                %add some info from the metabolite fit
-                dataToFit_mm.lineShape  = fitParamsOff.lineShape;
-                dataToFit_mm.refFWHM  = fitParamsOff.refFWHM;
-                % Extract fit options
-                fitOpts_mm    = MRSCont.opts.fit;
-                fitModel_mm    = fitOpts.method;
-                fitOpts_mm.sequence = 'MEGA';
-                %Specify a reduced basis set for MM modeling
-                %basisSet_mm    = MRSCont.fit.basisSet;
-                %Reduce the size of the basis set
-
-                %Adjust basis set
-                % Clear existing basis set
-                MRSCont.fit.basisSet_mm = [];
-                basisSet_mm = basisSetDiff1;
-                % Generate the list of basis functions that are supposed to be included in
-                % the basis set
-                % To do: Interface with interactive user input
-                metabList_mm = fit_createMetabListMM('MEGA');
-                % Create the modified basis set
-                basisSet_mm = fit_selectMetabs(basisSet_mm, metabList_mm, 0);
-                 % Call the fit function
-                [fitParams_mm, resBasisSet_mm] = fit_runFitMM(dataToFit_mm, basisSet_mm, fitModel_mm, fitOpts_mm);
-                 % Save back the basis set and fit parameters to MRSCont
-                MRSCont.fit.basisSet_mm                    = basisSet_mm;
-                MRSCont.fit.resBasisSet.diff1_mm{kk}         = resBasisSet_mm;
-                MRSCont.fit.results.diff1_mm.fitParams{kk}   = fitParams_mm;
-            end  
-        
         end
     end
 
diff --git a/fit/osp_fitUnEdited.m b/fit/osp_fitUnEdited.m
index 233c08c0..091b7087 100644
--- a/fit/osp_fitUnEdited.m
+++ b/fit/osp_fitUnEdited.m
@@ -34,12 +34,12 @@
     progressText = '';
 end
 for kk = 1:MRSCont.nDatasets
-
-
+    
+        
     % ----- Osprey fit pipeline -----
     if strcmpi(MRSCont.opts.fit.method, 'Osprey')
-        [~] = printLog('OspreyFit', kk, MRSCont.nDatasets, progressText, MRSCont.flags.isGUI, MRSCont.flags.isMRSI);
-
+        [~] = printLog('OspreyFit', kk, MRSCont.nDatasets, progressText, MRSCont.flags.isGUI, MRSCont.flags.isMRSI);    
+        
         if ~(MRSCont.flags.didFit == 1 && MRSCont.flags.speedUp && isfield(MRSCont, 'fit') && (kk > length(MRSCont.fit.results.off.fitParams))) || ~strcmp(MRSCont.ver.Osp,MRSCont.ver.CheckOsp)
             % Apply scaling factor to the data
             dataToFit   = MRSCont.processed.A{kk};
@@ -47,7 +47,7 @@
             % Extract fit options
             fitOpts     = MRSCont.opts.fit;
             fitModel    = fitOpts.method;
-
+            
             % If MRSI data load priors
             if MRSCont.flags.isMRSI
                 if isfield(MRSCont.fit, 'results')
@@ -57,11 +57,11 @@
                     fitOpts.MRSIpriors = fitParamsOff;
                 end
             end
-
+            
             % Call the fit function
             basisSet    = MRSCont.fit.basisSet;
             [fitParams, resBasisSet] = fit_runFit(dataToFit, basisSet, fitModel, fitOpts);
-
+            
             % Save back the basis set and fit parameters to MRSCont
             MRSCont.fit.basisSet                    = basisSet;
             MRSCont.fit.resBasisSet.off{kk}         = resBasisSet;
@@ -76,11 +76,10 @@
                 % Extract fit options
                 fitOpts_mm    = MRSCont.opts.fit;
                 fitModel_mm    = fitOpts.method;
-                fitOpts_mm.sequence = 'unedited';
                 %Specify a reduced basis set for MM modeling
                 %basisSet_mm    = MRSCont.fit.basisSet;
                 %Reduce the size of the basis set
-
+                
                 %Adjust basis set
                 % Clear existing basis set
                 MRSCont.fit.basisSet_mm = [];
@@ -90,7 +89,7 @@
                 % Generate the list of basis functions that are supposed to be included in
                 % the basis set
                 % To do: Interface with interactive user input
-                metabList_mm = fit_createMetabListMM('unedited');
+                metabList_mm = fit_createMetabListMM;
                 % Collect MMfit flag from the options determined in the job file
                 fitMM = MRSCont.opts.fit.fitMM;
                 % Create the modified basis set
@@ -104,22 +103,22 @@
                 MRSCont.fit.basisSet_mm;
             end                                         % re_mm
         end
-
+        
         % ----- LCModel wrapper fit pipeline -----
     elseif strcmpi(MRSCont.opts.fit.method, 'LCModel')
         [~] = printLog('OspreyFit', kk, MRSCont.nDatasets, progressText, MRSCont.flags.isGUI, MRSCont.flags.isMRSI,1);
-
-
+    
+        
         % Put the scale to be 1 for now. Might have to adjust.
         MRSCont.fit.scale{kk} = 1;
-
+        
         % Create the sub-folder where the LCModel results will be saved,
         % otherwise LCModel will throw 'FATAL ERROR MAIN 2'.
         if ~exist(fullfile(MRSCont.outputFolder, 'LCMoutput'), 'dir')
             mkdir(fullfile(MRSCont.outputFolder, 'LCMoutput'));
         end
-
-
+        
+        
         % Find the path to LCModel binary
         pathLCModelBinaryStruct = osp_platform('lcmodel');
         pathLCModelBinary = fullfile(MRSCont.ospFolder, 'libraries', 'LCModel',pathLCModelBinaryStruct.os,pathLCModelBinaryStruct.osver,['LCModel_' pathLCModelBinaryStruct.os,'_',pathLCModelBinaryStruct.osver]);
@@ -131,24 +130,24 @@
             case 'win'
                 bin = 'LCModel.exe';
         end
-
+        
         if ~exist([pathLCModelBinary filesep bin],'file') %Unzip binary
             unzip([pathLCModelBinary filesep bin '.zip'],[pathLCModelBinary filesep]);
         end
-
+        
         if ~exist([pathLCModelBinary filesep bin],'file') %Binary is still missing
             error('ERROR: No LCModel binary found.  Aborting!!');
         else
             addpath(which('libraries/LCModel'));
-        end
-
+        end    
+        
         callLCModel(MRSCont, MRSCont.opts.fit.lcmodel.controlfileA{kk},[pathLCModelBinary filesep bin]);
-
+         
         % Save the parameters and information about the basis set
         MRSCont.fit.results.off.fitParams{kk} = readLCMFitParams(MRSCont, 'A', kk);
-
+        
     end
-
+    
 end
 
 time = toc(metFitTime);
@@ -238,7 +237,7 @@ function callLCModel(MRSCont, controlFile,pathLCModelBinary)
 end
 % Store amplitudes regardless whether the are fit or not
 fitParams.ampl        = tab.concentration';
-
+    
 % Read the raw area of the unsuppressed water peak from the .print file
 infoPrint = mrs_readLcmodelPRINT( [MRSCont.opts.fit.lcmodel.(lcmOutputFile){kk} '.print'] );
 if isfield(infoPrint,'h2oarea')
@@ -246,5 +245,6 @@ function callLCModel(MRSCont, controlFile,pathLCModelBinary)
 end
 
 
-
+  
 end
+
diff --git a/libraries/FID-A/fitTools/fitModels/Osprey/fit_OspreyParamsToModel.m b/libraries/FID-A/fitTools/fitModels/Osprey/fit_OspreyParamsToModel.m
index c5fbb4ab..81808d3d 100644
--- a/libraries/FID-A/fitTools/fitModels/Osprey/fit_OspreyParamsToModel.m
+++ b/libraries/FID-A/fitTools/fitModels/Osprey/fit_OspreyParamsToModel.m
@@ -36,14 +36,14 @@
 dataToFit     = inputData.dataToFit;
 dataToFit     = op_zeropad(dataToFit, 2); % zero-fill for LCModel
 basisSet      = inputData.basisSet;
-% if (length(fitParams.ampl) == 2)
-%     basisSet      = inputData.basisSet_mm;
-%     fitParams.freqShift = repmat(fitParams.freqShift,[basisSet.nMets+basisSet.nMM 1]);
-%     fitParams.lorentzLB = repmat(fitParams.lorentzLB,[basisSet.nMets+basisSet.nMM 1]);
-%     %dummy=fitParams.ampl;   
-%     %fitParams.ampl=zeros([basisSet.nMets+basisSet.nMM 1]);
-%     %fitParams.ampl([3 4 13])=dummy;
-% end
+if (length(fitParams.ampl) == 3)
+    basisSet      = inputData.basisSet_mm;
+    fitParams.freqShift = repmat(fitParams.freqShift,[basisSet.nMets+basisSet.nMM 1]);
+    fitParams.lorentzLB = repmat(fitParams.lorentzLB,[basisSet.nMets+basisSet.nMM 1]);
+    %dummy=fitParams.ampl;   
+    %fitParams.ampl=zeros([basisSet.nMets+basisSet.nMM 1]);
+    %fitParams.ampl([3 4 13])=dummy;
+end
 % ... settings:
 fitRangePPM         = inputSettings.fitRangePPM;
 minKnotSpacingPPM   = inputSettings.minKnotSpacingPPM;
diff --git a/libraries/FID-A/fitTools/fitModels/Osprey/fit_OspreyPrelimStep2MM.m b/libraries/FID-A/fitTools/fitModels/Osprey/fit_OspreyPrelimStep2MM.m
index 63ce9331..0e1f696b 100644
--- a/libraries/FID-A/fitTools/fitModels/Osprey/fit_OspreyPrelimStep2MM.m
+++ b/libraries/FID-A/fitTools/fitModels/Osprey/fit_OspreyPrelimStep2MM.m
@@ -164,10 +164,10 @@
 % Set the hard box constraints for the parameters
 nMets   = resBasisSet.nMets;
 nMM     = resBasisSet.nMM;
-lb_ph0              = -7.5; 
-ub_ph0              = +7.5; % Zero order phase shift [deg]
-lb_ph1              = -2.5; 
-ub_ph1              = +2.5; % First order phase shift [deg/ppm]
+lb_ph0              = -45; 
+ub_ph0              = +45; % Zero order phase shift [deg]
+lb_ph1              = -10; 
+ub_ph1              = +10; % First order phase shift [deg/ppm]
 lb_gaussLB          = 0; 
 ub_gaussLB          = sqrt(5000); % Gaussian dampening [Hz]
 lb_lorentzLB_mets   = zeros(nMets, 1); 
diff --git a/libraries/FID-A/fitTools/fitModels/Osprey/fit_Osprey_PrelimReducedMM.m b/libraries/FID-A/fitTools/fitModels/Osprey/fit_Osprey_PrelimReducedMM.m
index fd1e515e..d7eb3478 100644
--- a/libraries/FID-A/fitTools/fitModels/Osprey/fit_Osprey_PrelimReducedMM.m
+++ b/libraries/FID-A/fitTools/fitModels/Osprey/fit_Osprey_PrelimReducedMM.m
@@ -100,10 +100,10 @@
 inputSettings.lorentzLB     = lorentzLB;
 
 % Set the hard box constraints for the parameters
-lb_ph0          = -5; 
-ub_ph0          = +5; % Zero order phase shift [deg]
-lb_ph1          = -2.5; 
-ub_ph1          = +2.5;  % First order phase shift [deg/ppm]
+lb_ph0          = -15; 
+ub_ph0          = +15; % Zero order phase shift [deg]
+lb_ph1          = -10; 
+ub_ph1          = +10;  % First order phase shift [deg/ppm]
 lb_gaussLB      = 0; 
 ub_gaussLB      = sqrt(5000); % Gaussian dampening [Hz]
 lb_freqShift    = -5; 
diff --git a/libraries/FID-A/fitTools/fit_createMetabListMM.m b/libraries/FID-A/fitTools/fit_createMetabListMM.m
index 0118c5f0..2c802dbe 100644
--- a/libraries/FID-A/fitTools/fit_createMetabListMM.m
+++ b/libraries/FID-A/fitTools/fit_createMetabListMM.m
@@ -17,90 +17,47 @@
 % INPUTS:
 % NONE
 
-function metabList = fit_createMetabListMM(sequence)
-switch sequence 
-    case 'unedited'
-        % Select metabolites to include in basis set
-        metabList.Ala      = 0;
-        metabList.Asc      = 0;
-        metabList.Asp      = 0;
-        metabList.bHB      = 0;
-        metabList.bHG      = 0;
-        metabList.Cit      = 0;
-        metabList.Cr       = 1;
-        metabList.CrCH2    = 1;
-        metabList.EtOH     = 0;
-        metabList.GABA     = 0;
-        metabList.GPC      = 0;
-        metabList.GSH      = 0;
-        metabList.Glc      = 0;
-        metabList.Gln      = 0;
-        metabList.Glu      = 0;
-        metabList.Gly      = 0;
-        metabList.H2O      = 1;
-        metabList.Ins      = 0;
-        metabList.Lac      = 0;
-        metabList.NAA      = 1;
-        metabList.NAAG     = 0;
-        metabList.PCh      = 0;
-        metabList.PCr      = 0;
-        metabList.PE       = 0;
-        metabList.Phenyl   = 0;
-        metabList.Scyllo   = 0;
-        metabList.Ser      = 0;
-        metabList.Tau      = 0;
-        metabList.Tyros    = 0;
+function metabList = fit_createMetabList
 
-        % Select MM/lipid basis functions to include
-        metabList.MM09     = 1;
-        metabList.MM12     = 1;
-        metabList.MM14     = 1;
-        metabList.MM17     = 1;
-        metabList.MM20     = 1;
-        metabList.Lip09    = 1;
-        metabList.Lip13    = 1;
-        metabList.Lip20    = 1;
-    case 'MEGA'
-        % Select metabolites to include in basis set
-        metabList.Ala      = 0;
-        metabList.Asc      = 0;
-        metabList.Asp      = 0;
-        metabList.bHB      = 0;
-        metabList.bHG      = 0;
-        metabList.Cit      = 0;
-        metabList.Cr       = 0;
-        metabList.CrCH2    = 0;
-        metabList.EtOH     = 0;
-        metabList.GABA     = 0;
-        metabList.GPC      = 0;
-        metabList.GSH      = 0;
-        metabList.Glc      = 0;
-        metabList.Gln      = 0;
-        metabList.Glu      = 0;
-        metabList.Gly      = 0;
-        metabList.H2O      = 0;
-        metabList.Ins      = 0;
-        metabList.Lac      = 0;
-        metabList.NAA      = 1;
-        metabList.NAAG     = 0;
-        metabList.PCh      = 0;
-        metabList.PCr      = 0;
-        metabList.PE       = 0;
-        metabList.Phenyl   = 0;
-        metabList.Scyllo   = 0;
-        metabList.Ser      = 0;
-        metabList.Tau      = 0;
-        metabList.Tyros    = 0;
+% Select metabolites to include in basis set
+metabList.Ala      = 0;
+metabList.Asc      = 0;
+metabList.Asp      = 0;
+metabList.bHB      = 0;
+metabList.bHG      = 0;
+metabList.Cit      = 0;
+metabList.Cr       = 1;
+metabList.CrCH2    = 1;
+metabList.EtOH     = 0;
+metabList.GABA     = 0;
+metabList.GPC      = 0;
+metabList.GSH      = 0;
+metabList.Glc      = 0;
+metabList.Gln      = 0;
+metabList.Glu      = 0;
+metabList.Gly      = 0;
+metabList.H2O      = 1;
+metabList.Ins      = 0;
+metabList.Lac      = 0;
+metabList.NAA      = 1;
+metabList.NAAG     = 0;
+metabList.PCh      = 0;
+metabList.PCr      = 0;
+metabList.PE       = 0;
+metabList.Phenyl   = 0;
+metabList.Scyllo   = 0;
+metabList.Ser      = 0;
+metabList.Tau      = 0;
+metabList.Tyros    = 0;
 
-        % Select MM/lipid basis functions to include
-        metabList.MM09     = 0;
-        metabList.MM12     = 0;
-        metabList.MM14     = 0;
-        metabList.MM17     = 0;
-        metabList.MM20     = 0;
-        metabList.Lip09    = 0;
-        metabList.Lip13    = 0;
-        metabList.Lip20    = 0;
-end
+% Select MM/lipid basis functions to include
+metabList.MM09     = 1;
+metabList.MM12     = 1;
+metabList.MM14     = 1;
+metabList.MM17     = 1;
+metabList.MM20     = 1;
+metabList.Lip09    = 1;
+metabList.Lip13    = 1;
+metabList.Lip20    = 1;
 
 end
\ No newline at end of file
diff --git a/libraries/FID-A/fitTools/fit_makeBasis.m b/libraries/FID-A/fitTools/fit_makeBasis.m
index 2b2aaac1..a1a4472b 100644
--- a/libraries/FID-A/fitTools/fit_makeBasis.m
+++ b/libraries/FID-A/fitTools/fit_makeBasis.m
@@ -89,11 +89,7 @@
     buffer.n(kk)                = temp.(basisFct{1}).sz(1);
     buffer.linewidth(kk)        = temp.(basisFct{1}).linewidth;
     buffer.Bo(kk)               = temp.(basisFct{1}).Bo;
-    if iscell(temp.(basisFct{1}).seq)
-        buffer.seq{kk}              = temp.(basisFct{1}).seq{1};
-    else
-        buffer.seq{kk}              = temp.(basisFct{1}).seq;
-    end
+    buffer.seq{kk}              = temp.(basisFct{1}).seq;
     if isfield(temp.(basisFct{1}),'name')
         buffer.name{kk}             = temp.(basisFct{1}).name;
     else
diff --git a/libraries/FID-A/fitTools/fit_runFitMM.m b/libraries/FID-A/fitTools/fit_runFitMM.m
index 2f83e8ff..8e181c81 100644
--- a/libraries/FID-A/fitTools/fit_runFitMM.m
+++ b/libraries/FID-A/fitTools/fit_runFitMM.m
@@ -32,17 +32,10 @@
 if nargin<4
     fitOpts = fit_defaultFitOpts(fitModel);
 end
-switch fitOpts.sequence
-    case 'unedited'
-        %Invert the NAA, Cr and CrCH2 peaks
-        basisSet.specs(:,1)=-basisSet.specs(:,1);
-        basisSet.specs(:,2)=-basisSet.specs(:,2);
-        basisSet.specs(:,4)=-basisSet.specs(:,4);
-    case 'MEGA'
-        %Invert the NAA peaks
-        basisSet.specs(:,1)=-basisSet.specs(:,1);
-end
-        
+%Invert the NAA, Cr and CrCH2 peaks
+basisSet.specs(:,1)=-basisSet.specs(:,1);
+basisSet.specs(:,2)=-basisSet.specs(:,2);
+basisSet.specs(:,4)=-basisSet.specs(:,4);
 %%% 1. SELECT THE MODEL %%%
 switch fitModel
     case 'Osprey'
diff --git a/plot/osp_plotFit.m b/plot/osp_plotFit.m
index 65994f06..498a86d8 100755
--- a/plot/osp_plotFit.m
+++ b/plot/osp_plotFit.m
@@ -315,9 +315,6 @@
 if (strcmp(which_spec, 'mm'))
     Met_corr_spectrum  = sum(ModelOutput.indivMets(:,1:4),2);
 end
-if (strcmp(which_spec, 'diff1_mm'))
-   Met_corr_spectrum  = sum(ModelOutput.indivMets(:,1:2),2);
-end
 
 
 %%% 4. SET UP FIGURE LAYOUT %%%
@@ -413,8 +410,8 @@
 plot(ppm, (zeros(1,length(ppm)) + max(dataToPlot +  abs(min(dataToPlot - fit))) + stagData)/maxPlot, 'Color', colorData, 'LineStyle','--', 'LineWidth', 0.5); % Zeroline Residue
 plot(ppm, (zeros(1,length(ppm)) + max(dataToPlot +  abs(min(dataToPlot - fit))) + abs(max(dataToPlot - fit)) + stagData)/maxPlot, 'Color', colorData, 'LineWidth', 1); % Max Residue
 
-if (contains(which_spec, 'mm'))
-   plot(ppm, (dataToPlot + stagData-Met_corr_spectrum)/maxPlot, 'Color',[1 0 0.1]); % Data
+if (strcmp(which_spec, 'mm'))
+    plot(ppm, (dataToPlot + stagData-Met_corr_spectrum)/maxPlot, 'Color',[1 0 0.1]); % Data
 end
 
 text(fitRangePPM(1), (0 + stagData)/maxPlot, '0', 'FontSize', 10,'Color', MRSCont.colormap.Foreground); %Zeroline text

From 59d548d2494c106845c7dee0c3d142bbfe766d20 Mon Sep 17 00:00:00 2001
From: Helge <47348963+HJZollner@users.noreply.github.com>
Date: Wed, 23 Mar 2022 18:01:16 -0400
Subject: [PATCH 3/4] Revert "Update SpecReg"

This reverts commit 4793f7065cb9ca224388bf61580ff8b4f1284c7b.
---
 .../processingTools/op_SpecRegFreqRestrict.m  |  15 +-
 .../FID-A/processingTools/op_probabSpecReg.m  | 357 ------------------
 .../FID-A/processingTools/op_robustSpecReg.m  |   5 +-
 3 files changed, 5 insertions(+), 372 deletions(-)
 delete mode 100644 libraries/FID-A/processingTools/op_probabSpecReg.m

diff --git a/libraries/FID-A/processingTools/op_SpecRegFreqRestrict.m b/libraries/FID-A/processingTools/op_SpecRegFreqRestrict.m
index 731c0a47..003c2f9d 100644
--- a/libraries/FID-A/processingTools/op_SpecRegFreqRestrict.m
+++ b/libraries/FID-A/processingTools/op_SpecRegFreqRestrict.m
@@ -41,12 +41,6 @@
     end
 end
 
-if strcmp(seqType,'HERMES') || strcmp(seqType,'HERMES')
-    [in, ~] = osp_onOffClassifyHERMES(in);
-    ppmMin = [0.5,0.5,1.85,1.85];
-    ppmMax = [3.9,3.8,3.9,3.8];
-end
-
 % Check whether data is coil-combined. If not, throw error.
 if ~in.flags.addedrcvrs
     error('ERROR:  I think it only makes sense to do this after you have combined the channels using op_addrcvrs.  ABORTING!!');
@@ -92,12 +86,11 @@
 % Initialize common variables and loop over sub-spectra
 t                 = in.t;
 input.dwelltime   = in.dwelltime;
-
-for mm=1:numSubSpecs
-
 in_restrict = op_freqrange(in,ppmMin,ppmMax);
 t_restrict                 = in_restrict.t;
-freq        = in_restrict.ppm;       
+freq        = in_restrict.ppm;
+for mm=1:numSubSpecs
+       
     
         DataToAlign = in_restrict.fids(:,:,mm);
         
@@ -268,7 +261,7 @@
 
 % Perform weighted averaging
 % No need for 'mean', since the weights vector w is normalized
-fids = squeeze(sum(fids_out, in.dims.averages));
+fids = sum(fids_out, in.dims.averages);
     
 %%% 3. WRITE THE NEW STRUCTURE %%%
 %re-calculate Specs using fft
diff --git a/libraries/FID-A/processingTools/op_probabSpecReg.m b/libraries/FID-A/processingTools/op_probabSpecReg.m
deleted file mode 100644
index 1329bba1..00000000
--- a/libraries/FID-A/processingTools/op_probabSpecReg.m
+++ /dev/null
@@ -1,357 +0,0 @@
-function [out, fs, phs, w, driftPre, driftPost] = op_probabSpecReg(in, seqType, echo, F0, noAlign)
-% Spectral registration is a time-domain frequency-and-phase correction
-% routine as per Near et al. (2015). Incorporates a multiplexed,
-% probabilistic approach for aligning HERMES data (MM: 170609)
-showPlots = 0;
-if nargin <5
-    noAlign = 0;
-    if nargin <4
-        F0 = nan;
-        if nargin < 2
-            echo = 1;    
-        end
-    end
-end
-
-% Check whether data is coil-combined. If not, throw error.
-if ~in.flags.addedrcvrs
-    error('ERROR:  I think it only makes sense to do this after you have combined the channels using op_addrcvrs.  ABORTING!!');
-end
-
-% Check whether data is averaged. If it is, throw warning and return unmodified data.
-if in.dims.averages == 0 || in.averages < 2
-    %DO NOTHING
-    warning('No averages found.  Returning input without modification!');
-    out = in;
-    out.flags.freqcorrected = 1;
-    fs  = 0;
-    phs = 0;
-    return
-end
-
-% Check whether input structure has sub-spectra.
-% If it does, prepare to loop over all sub-spectra. Each sub-spectrum will
-% get aligned separately.
-if in.dims.subSpecs == 0
-    numSubSpecs = 1;
-    % Measure drift pre-alignment
-    driftPre = op_measureDrift(in);
-else
-    numSubSpecs = in.sz(in.dims.subSpecs);
-    % Measure drift pre-alignment
-    for mm = 1:numSubSpecs
-        in_sub = op_takesubspec(in, mm);
-        driftPre{mm} = op_measureDrift(in_sub);
-    end
-end
-
-
-% Pre-allocate memory.
-fids    = zeros(in.sz(in.dims.t),1,numSubSpecs);
-fs      = zeros(in.sz(in.dims.averages),numSubSpecs);
-phs     = zeros(in.sz(in.dims.averages),numSubSpecs);
-zMSE = zeros(numSubSpecs,in.sz(in.dims.averages));
-CorrParsML = zeros(in.sz(in.dims.averages),2);
-parsGuess = [0 0];
-
-% Initialize common variables and loop over sub-spectra
-t                 = in.t;
-input.dwelltime   = in.dwelltime;
-
-% Probability density function and parameter bounds
-Cauchy = @(x,s,l) s./(pi.*(s.^2+(x-l).^2));
-lb = [0 -Inf];
-ub = [Inf Inf];
-
-% Optimization options
-lsqnonlinopts = optimoptions(@lsqnonlin);
-lsqnonlinopts = optimoptions(lsqnonlinopts,'Display','off','Algorithm','levenberg-marquardt');
-% nlinopts = statset('nlinfit');
-% nlinopts = statset(nlinopts,'MaxIter',400,'TolX',1e-8,'TolFun',1e-8);
-mleopts = statset('mlecustom');
-mleopts = statset(mleopts,'MaxIter',400,'MaxFunEvals',800,'TolX',1e-6,'TolFun',1e-6,'TolBnd',1e-6);
-
-% Set dimensions of figures of histograms
-if showPlots == 1
-    d.w = 0.6;
-    d.h = 0.45;
-    d.l = (1-d.w)/2;
-    d.b = (1-d.h)/2;
-end
-
-for mm=1:numSubSpecs
-    clear m;
-    freq        = in.ppm;
-    
-    DataToAlign = in.fids(:,:,mm);
-    
-    
-    % Use first n points of time-domain data, where n is the last point where SNR > 3
-    signal = abs(DataToAlign(:,mm));
-    noise = 2*std(signal(ceil(0.75*size(signal,1)):end,:));
-    SNR = signal ./ repmat(noise, [size(DataToAlign,1) 1]);
-    SNR = abs(diff(mean(SNR,2)));
-    
-    %SNR = abs(diff(mean(SNR,2))); This is how it's done in RobSpecReg
-    
-%     n = find(SNR > 1.5); % This is original Gannet
-%     if isempty(n)
-%         tMax = 100;
-%     else
-%         tMax = n(end);
-%         if tMax < 100
-%             tMax = 100;
-%         end
-%     end
-
-    SNR = SNR(t <= 0.2);
-    tMax = find(SNR > 1.5,1,'last');
-    if isempty(tMax) || tMax < find(t <= 0.1,1,'last')
-        tMax = find(t <= 0.1,1,'last');
-    end
-    
-    % 'Flatten' complex data for use in nlinfit
-    clear flatdata;
-    flatdata(:,1,:) = real(DataToAlign(1:tMax,:));
-    flatdata(:,2,:) = imag(DataToAlign(1:tMax,:));
-    
-    % Reference transient
-    flattarget = median(flatdata,3); % median across transients
-    target = flattarget(:);
-    
-    % Scalar to normalize transients (reduces optimization time)
-    a = max(abs(target));
-    
-    % Pre-allocate memory
-    if mm == 1
-        params = zeros(size(flatdata,3), 2);
-        MSE = zeros(1, size(flatdata,3));
-    end
-    
-     % Frame-by-frame determination of frequency of residual water and Cr (if HERMES or GSH editing)
-    if strcmp(seqType, 'HERMES') || strcmp(seqType, 'HERCULES')
-        F0freqRange = freq - 3.02 >= -0.15 & freq - 3.02 <= 0.15;
-    else
-        F0freqRange = freq - 3.22 >= -0.15 & freq - 3.22 <= 0.15;
-    end
-    [~,FrameMaxPos] = max(abs(real(in.specs(F0freqRange,:))),[],1);
-    F0freqRange = freq(F0freqRange);
-    F0freq = F0freqRange(FrameMaxPos);
-
-    % Starting values for optimization
-    if isnan(F0) 
-        f0 = F0freq * in.txfrq * 1e-6;
-        f0 = f0 - f0(1);
-    else
-        f0 = F0;
-        f0 = f0 - f0(1);
-    end
-    phi0 = zeros(size(f0));
-    x0 = [f0(:) phi0(:)];
-        
-    
-    % Determine frequency and phase offsets by spectral registration
-     if noAlign == 0
-        time = 0:input.dwelltime:(length(target)/2 - 1)*input.dwelltime;
-        reverseStr = '';
-        for jj = 1:size(flatdata,3)
-            if echo
-                msg = sprintf('\nSpectral registration - Fitting transient: %d', corrloop);
-                fprintf([reverseStr, msg]);
-                reverseStr = repmat(sprintf('\b'), 1, length(msg));
-            end
-
-            transient = squeeze(flatdata(:,:,jj));
-
-            fun = @(x) SpecReg(transient(:)/a, target/a, time, x);
-            params(jj,:) = lsqnonlin(fun, x0(jj,:), [], [], lsqnonlinopts);
-
-            f   = params(jj,1);
-            phi = params(jj,2);
-            m_c = complex(flatdata(:,1,jj), flatdata(:,2,jj));
-            m_c = m_c .* exp(1i*pi*(time'*f*2+phi/180));
-            m(:,jj) = [real(m_c); imag(m_c)];
-            resid = target - m(:,jj);
-            MSE(jj) = sum(resid.^2) / (length(resid) - 2);
-
-%             x0(jj,2) = phi; % Update phase value according to the last iteration
-
-    %         [parsFit(corrloop,:), ~, ~, ~, MSE(corrloop)] = nlinfit(input, target, @FreqPhaseShiftNest, parsGuess, nlinopts);
-    %         parsGuess = parsFit(corrloop,:);
-        end
-         % Prepare the frequency and phase corrections to be returned by this
-        % function for further use
-        fs(:,mm)  = params(:,1);
-        phs(:,mm) = params(:,2);
-    end
-    
-    
-    % Probability distribution of frequency offsets (estimated by maximum likelihood)
-    start = [iqr(fs(:,mm))/2, median(fs(:,mm))];
-    [fs_p(:,mm), fs_p_ci(:,:,mm)] = ...
-        mle(fs(:,mm), 'pdf', Cauchy, 'start', start, 'lower', lb, 'upper', ub, 'options', mleopts);
-    fs_fx(:,mm) = ...
-        linspace(1.5*min(fs(:,mm)), 1.5*max(fs(:,mm)), 1e3);
-    fs_pdf(:,mm) = Cauchy(fs_fx(:,mm), fs_p(1,mm), fs_p(2,mm));
-    
-    % Probability distribution of phase offsets (estimated by maximum likelihood)
-    start = [iqr(phs(:,mm))/2, median(phs(:,mm))];
-    [phs_p(:,mm), phs_p_ci(:,:,mm)] = ...
-        mle(phs(:,mm), 'pdf', Cauchy, 'start', start, 'lower', lb, 'upper', ub, 'options', mleopts);
-    phs_fx(:,mm) = ...
-        linspace(1.5*min(phs(:,mm)), 1.5*max(phs(:,mm)), 1e3);
-    phs_pdf(:,mm) = Cauchy(phs_fx(:,mm), phs_p(1,mm), phs_p(2,mm));
-             
-
-    zMSE = zscore(MSE); % standardized MSEs
-    
-    % Apply frequency and phase corrections
-     for jj = 1:size(flatdata,3)
-        % Default correction
-%         fids(:,jj,mm) = in.fids(:,jj,mm) .* ...
-%             exp(1i*params(jj,1)*2*pi*t') * exp(1i*pi/180*params(jj,2));
-
-        % Freq/phase correction + Cauchy pdf location parameter shift
-        fids(:,jj,mm) = in.fids(:,jj,mm) .* ...
-            exp(1i*(params(jj,1)-fs_p(2))*2*pi*t') * exp(1i*pi/180*(params(jj,2)-phs_p(2)));
-     end
-
-           
-        
-end
-
-if echo
-fprintf('\n');
-end
-
-%%% 2. CALCULATE THE DRIFT AFTER CORRECTIONS
-out_temp=in;
-out_temp.fids=fids;
-%re-calculate Specs using fft
-out_temp.specs=fftshift(fft(fids,[],in.dims.t),in.dims.t);
-% Measure drift post-alignment
-if in.dims.subSpecs == 0
-    numSubSpecs = 1;
-    driftPost = op_measureDrift(out_temp);
-else
-    numSubSpecs = in.sz(in.dims.subSpecs);
-    % Measure drift post-alignment
-    for mm = 1:numSubSpecs
-        out_temp_sub = op_takesubspec(out_temp, mm);
-        driftPost{mm} = op_measureDrift(out_temp_sub);
-    end
-end
-
-
-%%% 3. RE-RUN THE RANKING ALGORITHM TO DETERMINE WEIGHTS %%%
-
-% Determine optimal iteration order by calculating a similarity metric (mean squared error)
-w = cell(1,numSubSpecs);
-
-for mm=1:numSubSpecs
-    DataToWeight = fids(:,:,mm);
-    D = zeros(size(DataToWeight,2));
-    for jj = 1:size(DataToWeight,2)
-        for kk = 1:size(DataToWeight,2)
-            tmp = sum((real(DataToWeight(1:tMax,jj)) - real(DataToWeight(1:tMax,kk))).^2) / tMax;
-            if tmp == 0
-                D(jj,kk) = NaN;
-            else
-                D(jj,kk) = tmp;
-            end
-        end
-    end
-    d = nanmedian(D);
-    if isnan(sum(d))        
-        d(isnan(d)) = 1;
-    end
-    w{mm} = 1./d.^2;
-    w{mm} = w{mm}/sum(w{mm});
-    w{mm} = repmat(w{mm}, [size(DataToWeight,1) 1]);
-    
-    
-    
-    % Apply the weighting
-fids_out(:,:,mm) = w{mm} .* DataToWeight;
-    
-end
-
-% Or remove based on zMSE score
-% for mm=1:numSubSpecs
-%     DataToWeight = fids(:,:,mm);
-%     w{mm} = zMSE < 3;
-%     w{mm} = w{mm}/sum(w{mm});
-%     w{mm} = repmat(w{mm}, [size(DataToWeight,1) 1]);
-%         
-%     % Apply the weighting
-%     fids_out(:,:,mm) = w{mm} .* DataToWeight;
-%     
-% end
-
-% Perform weighted averaging
-% No need for 'mean', since the weights vector w is normalized
-fids = squeeze(sum(fids_out, in.dims.averages));
-    
-%%% 3. WRITE THE NEW STRUCTURE %%%
-%re-calculate Specs using fft
-specs=fftshift(fft(fids,[],in.dims.t),in.dims.t);
-
-%change the dims variables.
-if in.dims.t>in.dims.averages
-    dims.t=in.dims.t-1;
-else
-    dims.t=in.dims.t;
-end
-if in.dims.coils>in.dims.averages
-    dims.coils=in.dims.coils-1;
-else
-    dims.coils=in.dims.coils;
-end
-dims.averages=0;
-if in.dims.subSpecs>in.dims.averages
-    dims.subSpecs=in.dims.subSpecs-1;
-else
-    dims.subSpecs=in.dims.subSpecs;
-end
-if in.dims.extras>in.dims.averages
-    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.averages=1;
-    
-%FILLING IN THE FLAGS
-out.flags=in.flags;
-out.flags.writtentostruct=1;
-out.flags.freqcorrected=1;
-out.flags.averaged=1;
-
-
-
-end
-
-function out = SpecReg(data, target, time, x)
-
-f   = x(1);
-phi = x(2);
-
-data = reshape(data, [length(data)/2 2]);
-fid = complex(data(:,1), data(:,2));
-
-fid = fid .* exp(1i*pi*(time'*f*2+phi/180));
-fid = [real(fid); imag(fid)];
-
-out = target - fid;
-
-end
-
diff --git a/libraries/FID-A/processingTools/op_robustSpecReg.m b/libraries/FID-A/processingTools/op_robustSpecReg.m
index e23fe1e0..76ca7182 100644
--- a/libraries/FID-A/processingTools/op_robustSpecReg.m
+++ b/libraries/FID-A/processingTools/op_robustSpecReg.m
@@ -51,9 +51,6 @@
     out.flags.freqcorrected = 1;
     fs  = 0;
     phs = 0;
-    w = 0;
-    driftPre = 0;
-    driftPost = 0;
     return
 end
 
@@ -324,7 +321,7 @@
 
 % Perform weighted averaging
 % No need for 'mean', since the weights vector w is normalized
-fids = squeeze(sum(fids_out, in.dims.averages));
+fids = sum(fids_out, in.dims.averages);
     
 %%% 3. WRITE THE NEW STRUCTURE %%%
 %re-calculate Specs using fft

From 6defefde3c9b964fbbdde9a79906aa2d99611388 Mon Sep 17 00:00:00 2001
From: Helge <47348963+HJZollner@users.noreply.github.com>
Date: Wed, 23 Mar 2022 18:15:56 -0400
Subject: [PATCH 4/4] Update SpecReg

pdate dimensions in specreg functions
---
 .../processingTools/op_SpecRegFreqRestrict.m  |  30 +-
 .../FID-A/processingTools/op_probabSpecReg.m  | 356 ++++++++++++++++++
 .../FID-A/processingTools/op_robustSpecReg.m  |  44 +--
 3 files changed, 396 insertions(+), 34 deletions(-)
 create mode 100644 libraries/FID-A/processingTools/op_probabSpecReg.m

diff --git a/libraries/FID-A/processingTools/op_SpecRegFreqRestrict.m b/libraries/FID-A/processingTools/op_SpecRegFreqRestrict.m
index 003c2f9d..5dd8463a 100644
--- a/libraries/FID-A/processingTools/op_SpecRegFreqRestrict.m
+++ b/libraries/FID-A/processingTools/op_SpecRegFreqRestrict.m
@@ -35,12 +35,18 @@
         if nargin <4
             F0 = nan;
             if nargin < 2
-                echo = 1;    
+                echo = 1;
             end
         end
     end
 end
 
+if strcmp(seqType,'HERMES') || strcmp(seqType,'HERMES')
+    [in, ~] = osp_onOffClassifyHERMES(in);
+    ppmMin = [0.5,0.5,1.85,1.85];
+    ppmMax = [3.9,3.8,3.9,3.8];
+end
+
 % Check whether data is coil-combined. If not, throw error.
 if ~in.flags.addedrcvrs
     error('ERROR:  I think it only makes sense to do this after you have combined the channels using op_addrcvrs.  ABORTING!!');
@@ -90,10 +96,10 @@
 t_restrict                 = in_restrict.t;
 freq        = in_restrict.ppm;
 for mm=1:numSubSpecs
-       
-    
+
+
         DataToAlign = in_restrict.fids(:,:,mm);
-        
+
         % Use first n points of time-domain data, where n is the last point where abs(diff(mean(SNR))) > 0.5
         signal = abs(DataToAlign);
         noise = 2*std(signal(ceil(0.75*size(signal,1)):end,:));
@@ -149,7 +155,7 @@
         F0freq = F0freqRange(FrameMaxPos);
 
         % Starting values for optimization
-        if isnan(F0) 
+        if isnan(F0)
             f0 = F0freq * in.txfrq * 1e-6;
             f0 = f0(alignOrd);
             f0 = f0 - f0(1);
@@ -199,7 +205,7 @@
         fs(:,mm)  = params(:,1);
         phs(:,mm) = params(:,2);
     end
-    
+
     % Apply frequency and phase corrections to raw data
     for jj = 1:size(flatdata,3)
         fids(:,jj,mm) = in.fids(:,jj,mm) .* ...
@@ -253,16 +259,16 @@
     w{mm} = 1./d.^2;
     w{mm} = w{mm}/sum(w{mm});
     w{mm} = repmat(w{mm}, [size(fids(:,:,mm),1) 1]);
-    
+
     % Apply the weighting
 fids_out(:,:,mm) = w{mm} .* fids(:,:,mm);
-    
+
 end
 
 % Perform weighted averaging
 % No need for 'mean', since the weights vector w is normalized
-fids = sum(fids_out, in.dims.averages);
-    
+fids = squeeze(sum(fids_out, in.dims.averages);
+
 %%% 3. WRITE THE NEW STRUCTURE %%%
 %re-calculate Specs using fft
 specs=fftshift(fft(fids,[],in.dims.t),in.dims.t);
@@ -292,7 +298,7 @@
 
 %re-calculate the sz variable
 sz=size(fids);
-    
+
 %FILLING IN DATA STRUCTURE
 out=in;
 out.fids=fids;
@@ -300,7 +306,7 @@
 out.sz=sz;
 out.dims=dims;
 out.averages=1;
-    
+
 %FILLING IN THE FLAGS
 out.flags=in.flags;
 out.flags.writtentostruct=1;
diff --git a/libraries/FID-A/processingTools/op_probabSpecReg.m b/libraries/FID-A/processingTools/op_probabSpecReg.m
new file mode 100644
index 00000000..a2c5e4eb
--- /dev/null
+++ b/libraries/FID-A/processingTools/op_probabSpecReg.m
@@ -0,0 +1,356 @@
+function [out, fs, phs, w, driftPre, driftPost] = op_probabSpecReg(in, seqType, echo, F0, noAlign)
+% Spectral registration is a time-domain frequency-and-phase correction
+% routine as per Near et al. (2015). Incorporates a multiplexed,
+% probabilistic approach for aligning HERMES data (MM: 170609)
+showPlots = 0;
+if nargin <5
+    noAlign = 0;
+    if nargin <4
+        F0 = nan;
+        if nargin < 2
+            echo = 1;
+        end
+    end
+end
+
+% Check whether data is coil-combined. If not, throw error.
+if ~in.flags.addedrcvrs
+    error('ERROR:  I think it only makes sense to do this after you have combined the channels using op_addrcvrs.  ABORTING!!');
+end
+
+% Check whether data is averaged. If it is, throw warning and return unmodified data.
+if in.dims.averages == 0 || in.averages < 2
+    %DO NOTHING
+    warning('No averages found.  Returning input without modification!');
+    out = in;
+    out.flags.freqcorrected = 1;
+    fs  = 0;
+    phs = 0;
+    return
+end
+
+% Check whether input structure has sub-spectra.
+% If it does, prepare to loop over all sub-spectra. Each sub-spectrum will
+% get aligned separately.
+if in.dims.subSpecs == 0
+    numSubSpecs = 1;
+    % Measure drift pre-alignment
+    driftPre = op_measureDrift(in);
+else
+    numSubSpecs = in.sz(in.dims.subSpecs);
+    % Measure drift pre-alignment
+    for mm = 1:numSubSpecs
+        in_sub = op_takesubspec(in, mm);
+        driftPre{mm} = op_measureDrift(in_sub);
+    end
+end
+
+
+% Pre-allocate memory.
+fids    = zeros(in.sz(in.dims.t),1,numSubSpecs);
+fs      = zeros(in.sz(in.dims.averages),numSubSpecs);
+phs     = zeros(in.sz(in.dims.averages),numSubSpecs);
+zMSE = zeros(numSubSpecs,in.sz(in.dims.averages));
+CorrParsML = zeros(in.sz(in.dims.averages),2);
+parsGuess = [0 0];
+
+% Initialize common variables and loop over sub-spectra
+t                 = in.t;
+input.dwelltime   = in.dwelltime;
+
+% Probability density function and parameter bounds
+Cauchy = @(x,s,l) s./(pi.*(s.^2+(x-l).^2));
+lb = [0 -Inf];
+ub = [Inf Inf];
+
+% Optimization options
+lsqnonlinopts = optimoptions(@lsqnonlin);
+lsqnonlinopts = optimoptions(lsqnonlinopts,'Display','off','Algorithm','levenberg-marquardt');
+% nlinopts = statset('nlinfit');
+% nlinopts = statset(nlinopts,'MaxIter',400,'TolX',1e-8,'TolFun',1e-8);
+mleopts = statset('mlecustom');
+mleopts = statset(mleopts,'MaxIter',400,'MaxFunEvals',800,'TolX',1e-6,'TolFun',1e-6,'TolBnd',1e-6);
+
+% Set dimensions of figures of histograms
+if showPlots == 1
+    d.w = 0.6;
+    d.h = 0.45;
+    d.l = (1-d.w)/2;
+    d.b = (1-d.h)/2;
+end
+
+for mm=1:numSubSpecs
+    clear m;
+    freq        = in.ppm;
+
+    DataToAlign = in.fids(:,:,mm);
+
+
+    % Use first n points of time-domain data, where n is the last point where SNR > 3
+    signal = abs(DataToAlign(:,mm));
+    noise = 2*std(signal(ceil(0.75*size(signal,1)):end,:));
+    SNR = signal ./ repmat(noise, [size(DataToAlign,1) 1]);
+    SNR = abs(diff(mean(SNR,2)));
+
+    %SNR = abs(diff(mean(SNR,2))); This is how it's done in RobSpecReg
+
+%     n = find(SNR > 1.5); % This is original Gannet
+%     if isempty(n)
+%         tMax = 100;
+%     else
+%         tMax = n(end);
+%         if tMax < 100
+%             tMax = 100;
+%         end
+%     end
+
+    SNR = SNR(t <= 0.2);
+    tMax = find(SNR > 1.5,1,'last');
+    if isempty(tMax) || tMax < find(t <= 0.1,1,'last')
+        tMax = find(t <= 0.1,1,'last');
+    end
+
+    % 'Flatten' complex data for use in nlinfit
+    clear flatdata;
+    flatdata(:,1,:) = real(DataToAlign(1:tMax,:));
+    flatdata(:,2,:) = imag(DataToAlign(1:tMax,:));
+
+    % Reference transient
+    flattarget = median(flatdata,3); % median across transients
+    target = flattarget(:);
+
+    % Scalar to normalize transients (reduces optimization time)
+    a = max(abs(target));
+
+    % Pre-allocate memory
+    if mm == 1
+        params = zeros(size(flatdata,3), 2);
+        MSE = zeros(1, size(flatdata,3));
+    end
+
+     % Frame-by-frame determination of frequency of residual water and Cr (if HERMES or GSH editing)
+    if strcmp(seqType, 'HERMES') || strcmp(seqType, 'HERCULES')
+        F0freqRange = freq - 3.02 >= -0.15 & freq - 3.02 <= 0.15;
+    else
+        F0freqRange = freq - 3.22 >= -0.15 & freq - 3.22 <= 0.15;
+    end
+    [~,FrameMaxPos] = max(abs(real(in.specs(F0freqRange,:))),[],1);
+    F0freqRange = freq(F0freqRange);
+    F0freq = F0freqRange(FrameMaxPos);
+
+    % Starting values for optimization
+    if isnan(F0)
+        f0 = F0freq * in.txfrq * 1e-6;
+        f0 = f0 - f0(1);
+    else
+        f0 = F0;
+        f0 = f0 - f0(1);
+    end
+    phi0 = zeros(size(f0));
+    x0 = [f0(:) phi0(:)];
+
+
+    % Determine frequency and phase offsets by spectral registration
+     if noAlign == 0
+        time = 0:input.dwelltime:(length(target)/2 - 1)*input.dwelltime;
+        reverseStr = '';
+        for jj = 1:size(flatdata,3)
+            if echo
+                msg = sprintf('\nSpectral registration - Fitting transient: %d', corrloop);
+                fprintf([reverseStr, msg]);
+                reverseStr = repmat(sprintf('\b'), 1, length(msg));
+            end
+
+            transient = squeeze(flatdata(:,:,jj));
+
+            fun = @(x) SpecReg(transient(:)/a, target/a, time, x);
+            params(jj,:) = lsqnonlin(fun, x0(jj,:), [], [], lsqnonlinopts);
+
+            f   = params(jj,1);
+            phi = params(jj,2);
+            m_c = complex(flatdata(:,1,jj), flatdata(:,2,jj));
+            m_c = m_c .* exp(1i*pi*(time'*f*2+phi/180));
+            m(:,jj) = [real(m_c); imag(m_c)];
+            resid = target - m(:,jj);
+            MSE(jj) = sum(resid.^2) / (length(resid) - 2);
+
+%             x0(jj,2) = phi; % Update phase value according to the last iteration
+
+    %         [parsFit(corrloop,:), ~, ~, ~, MSE(corrloop)] = nlinfit(input, target, @FreqPhaseShiftNest, parsGuess, nlinopts);
+    %         parsGuess = parsFit(corrloop,:);
+        end
+         % Prepare the frequency and phase corrections to be returned by this
+        % function for further use
+        fs(:,mm)  = params(:,1);
+        phs(:,mm) = params(:,2);
+    end
+
+
+    % Probability distribution of frequency offsets (estimated by maximum likelihood)
+    start = [iqr(fs(:,mm))/2, median(fs(:,mm))];
+    [fs_p(:,mm), fs_p_ci(:,:,mm)] = ...
+        mle(fs(:,mm), 'pdf', Cauchy, 'start', start, 'lower', lb, 'upper', ub, 'options', mleopts);
+    fs_fx(:,mm) = ...
+        linspace(1.5*min(fs(:,mm)), 1.5*max(fs(:,mm)), 1e3);
+    fs_pdf(:,mm) = Cauchy(fs_fx(:,mm), fs_p(1,mm), fs_p(2,mm));
+
+    % Probability distribution of phase offsets (estimated by maximum likelihood)
+    start = [iqr(phs(:,mm))/2, median(phs(:,mm))];
+    [phs_p(:,mm), phs_p_ci(:,:,mm)] = ...
+        mle(phs(:,mm), 'pdf', Cauchy, 'start', start, 'lower', lb, 'upper', ub, 'options', mleopts);
+    phs_fx(:,mm) = ...
+        linspace(1.5*min(phs(:,mm)), 1.5*max(phs(:,mm)), 1e3);
+    phs_pdf(:,mm) = Cauchy(phs_fx(:,mm), phs_p(1,mm), phs_p(2,mm));
+
+
+    zMSE = zscore(MSE); % standardized MSEs
+
+    % Apply frequency and phase corrections
+     for jj = 1:size(flatdata,3)
+        % Default correction
+%         fids(:,jj,mm) = in.fids(:,jj,mm) .* ...
+%             exp(1i*params(jj,1)*2*pi*t') * exp(1i*pi/180*params(jj,2));
+
+        % Freq/phase correction + Cauchy pdf location parameter shift
+        fids(:,jj,mm) = in.fids(:,jj,mm) .* ...
+            exp(1i*(params(jj,1)-fs_p(2))*2*pi*t') * exp(1i*pi/180*(params(jj,2)-phs_p(2)));
+     end
+
+
+
+end
+
+if echo
+fprintf('\n');
+end
+
+%%% 2. CALCULATE THE DRIFT AFTER CORRECTIONS
+out_temp=in;
+out_temp.fids=fids;
+%re-calculate Specs using fft
+out_temp.specs=fftshift(fft(fids,[],in.dims.t),in.dims.t);
+% Measure drift post-alignment
+if in.dims.subSpecs == 0
+    numSubSpecs = 1;
+    driftPost = op_measureDrift(out_temp);
+else
+    numSubSpecs = in.sz(in.dims.subSpecs);
+    % Measure drift post-alignment
+    for mm = 1:numSubSpecs
+        out_temp_sub = op_takesubspec(out_temp, mm);
+        driftPost{mm} = op_measureDrift(out_temp_sub);
+    end
+end
+
+
+%%% 3. RE-RUN THE RANKING ALGORITHM TO DETERMINE WEIGHTS %%%
+
+% Determine optimal iteration order by calculating a similarity metric (mean squared error)
+w = cell(1,numSubSpecs);
+
+for mm=1:numSubSpecs
+    DataToWeight = fids(:,:,mm);
+    D = zeros(size(DataToWeight,2));
+    for jj = 1:size(DataToWeight,2)
+        for kk = 1:size(DataToWeight,2)
+            tmp = sum((real(DataToWeight(1:tMax,jj)) - real(DataToWeight(1:tMax,kk))).^2) / tMax;
+            if tmp == 0
+                D(jj,kk) = NaN;
+            else
+                D(jj,kk) = tmp;
+            end
+        end
+    end
+    d = nanmedian(D);
+    if isnan(sum(d))
+        d(isnan(d)) = 1;
+    end
+    w{mm} = 1./d.^2;
+    w{mm} = w{mm}/sum(w{mm});
+    w{mm} = repmat(w{mm}, [size(DataToWeight,1) 1]);
+
+
+
+    % Apply the weighting
+fids_out(:,:,mm) = w{mm} .* DataToWeight;
+
+end
+
+% Or remove based on zMSE score
+% for mm=1:numSubSpecs
+%     DataToWeight = fids(:,:,mm);
+%     w{mm} = zMSE < 3;
+%     w{mm} = w{mm}/sum(w{mm});
+%     w{mm} = repmat(w{mm}, [size(DataToWeight,1) 1]);
+%
+%     % Apply the weighting
+%     fids_out(:,:,mm) = w{mm} .* DataToWeight;
+%
+% end
+
+% Perform weighted averaging
+% No need for 'mean', since the weights vector w is normalized
+fids = sum(fids_out, in.dims.averages);
+
+%%% 3. WRITE THE NEW STRUCTURE %%%
+%re-calculate Specs using fft
+specs=fftshift(fft(fids,[],in.dims.t),in.dims.t);
+
+%change the dims variables.
+if in.dims.t>in.dims.averages
+    dims.t=in.dims.t-1;
+else
+    dims.t=in.dims.t;
+end
+if in.dims.coils>in.dims.averages
+    dims.coils=in.dims.coils-1;
+else
+    dims.coils=in.dims.coils;
+end
+dims.averages=0;
+if in.dims.subSpecs>in.dims.averages
+    dims.subSpecs=in.dims.subSpecs-1;
+else
+    dims.subSpecs=in.dims.subSpecs;
+end
+if in.dims.extras>in.dims.averages
+    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.averages=1;
+
+%FILLING IN THE FLAGS
+out.flags=in.flags;
+out.flags.writtentostruct=1;
+out.flags.freqcorrected=1;
+out.flags.averaged=1;
+
+
+
+end
+
+function out = SpecReg(data, target, time, x)
+
+f   = x(1);
+phi = x(2);
+
+data = reshape(data, [length(data)/2 2]);
+fid = complex(data(:,1), data(:,2));
+
+fid = fid .* exp(1i*pi*(time'*f*2+phi/180));
+fid = [real(fid); imag(fid)];
+
+out = target - fid;
+
+end
diff --git a/libraries/FID-A/processingTools/op_robustSpecReg.m b/libraries/FID-A/processingTools/op_robustSpecReg.m
index 76ca7182..58af9226 100644
--- a/libraries/FID-A/processingTools/op_robustSpecReg.m
+++ b/libraries/FID-A/processingTools/op_robustSpecReg.m
@@ -33,7 +33,7 @@
     if nargin <4
         F0 = nan;
         if nargin < 2
-            echo = 1;    
+            echo = 1;
         end
     end
 end
@@ -85,7 +85,7 @@
 input.dwelltime   = in.dwelltime;
 
 for mm=1:numSubSpecs
-   
+
         %%% Automatic lipid/unstable residual water removal %%%
         freq        = in.ppm;
         waterLim = freq <= 4.68 + 0.25 & freq >= 4.68 - 0.25;
@@ -206,7 +206,7 @@
         F0freq = F0freqRange(FrameMaxPos);
 
         % Starting values for optimization
-        if isnan(F0) 
+        if isnan(F0)
             f0 = F0freq * in.txfrq * 1e-6;
             f0 = f0(alignOrd);
             f0 = f0 - f0(1);
@@ -256,7 +256,7 @@
         fs(:,mm)  = params(:,1);
         phs(:,mm) = params(:,2);
     end
-    
+
     % Apply frequency and phase corrections to raw data
     for jj = 1:size(flatdata,3)
         fids(:,jj,mm) = in.fids(:,jj,mm) .* ...
@@ -305,24 +305,24 @@
         end
     end
     d = nanmedian(D);
-    if isnan(sum(d))        
+    if isnan(sum(d))
         d(isnan(d)) = 1;
     end
     w{mm} = 1./d.^2;
     w{mm} = w{mm}/sum(w{mm});
     w{mm} = repmat(w{mm}, [size(DataToWeight,1) 1]);
-    
-    
-    
+
+
+
     % Apply the weighting
 fids_out(:,:,mm) = w{mm} .* DataToWeight;
-    
+
 end
 
 % Perform weighted averaging
 % No need for 'mean', since the weights vector w is normalized
-fids = sum(fids_out, in.dims.averages);
-    
+fids = squeeze(sum(fids_out, in.dims.averages));
+
 %%% 3. WRITE THE NEW STRUCTURE %%%
 %re-calculate Specs using fft
 specs=fftshift(fft(fids,[],in.dims.t),in.dims.t);
@@ -352,7 +352,7 @@
 
 %re-calculate the sz variable
 sz=size(fids);
-    
+
 %FILLING IN DATA STRUCTURE
 out=in;
 out.fids=fids;
@@ -360,7 +360,7 @@
 out.sz=sz;
 out.dims=dims;
 out.averages=1;
-    
+
 %FILLING IN THE FLAGS
 out.flags=in.flags;
 out.flags.writtentostruct=1;
@@ -392,7 +392,7 @@
 %%% BELOW ARE THE SIGNAL FILTERING FUNCTIONS
 
 function out = SignalFilter(spec, r_flag, q_flag, in)
-% Based on: Cobas, J.C., Bernstein, M.A., Martín-Pastor, M., Tahoces, P.G.,
+% Based on: Cobas, J.C., Bernstein, M.A., Mart�n-Pastor, M., Tahoces, P.G.,
 % 2006. A new general-purpose fully automatic baseline-correction procedure
 % for 1D and 2D NMR data. J. Magn. Reson. 183, 145-51.
 
@@ -401,41 +401,41 @@
 [z.real, z.imag] = BaselineModeling(spec, r_flag, q_flag, in);
 
 if showPlots == 1
-    
+
     freq = in.ppm;
-    
+
     H = figure(77);
     d.w = 0.5;
     d.h = 0.65;
     d.l = (1-d.w)/2;
     d.b = (1-d.h)/2;
     set(H,'Color', 'w', 'Units', 'Normalized', 'OuterPosition', [d.l d.b d.w d.h]);
-    
+
     subplot(2,2,1); cla;
     plot(freq, real(spec), 'k');
     set(gca, 'XDir', 'reverse', 'XLim', [0 6], 'Box', 'off');
     hold on;
     plot(freq, z.real, 'r');
     set(gca, 'XDir', 'reverse', 'XLim', [0 6], 'Box', 'off');
-    
+
     subplot(2,2,2); cla;
     plot(freq, imag(spec), 'k');
     set(gca, 'XDir', 'reverse', 'XLim', [0 6], 'Box', 'off');
     hold on;
     plot(freq, z.imag, 'r');
     set(gca, 'XDir', 'reverse', 'XLim', [0 6], 'Box', 'off');
-    
+
     subplot(2,2,3); cla;
     plot(freq, real(spec) - z.real, 'k');
     set(gca, 'XDir', 'reverse', 'XLim', [0 6], 'Box', 'off');
-    
+
     subplot(2,2,4); cla;
     plot(freq, imag(spec) - z.imag, 'k');
     set(gca, 'XDir', 'reverse', 'XLim', [0 6], 'Box', 'off');
-    
+
     drawnow;
     %pause(0.5);
-    
+
 end
 
 out = ifft(fftshift(complex(real(spec) - z.real, imag(spec) - z.imag),in.dims.t),[],in.dims.t);