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
index cd93d93f..c8656aa6 100644
--- a/libraries/FID-A/processingTools/op_probabSpecReg.m
+++ b/libraries/FID-A/processingTools/op_probabSpecReg.m
@@ -82,18 +82,13 @@
 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;
@@ -109,25 +104,24 @@
     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;
@@ -139,7 +133,8 @@
     F0freq = F0freqRange(FrameMaxPos);
 
     % Starting values for optimization
-    if isnan(F0) 
+    if isnan(F0)
+
         f0 = F0freq * in.txfrq * 1e-6;
         f0 = f0 - f0(1);
     else
@@ -148,8 +143,7 @@
     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;
@@ -184,8 +178,7 @@
         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)] = ...
@@ -193,7 +186,7 @@
     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)] = ...
@@ -201,10 +194,11 @@
     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
@@ -215,9 +209,6 @@
         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
@@ -262,18 +253,18 @@
         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
 
 % Or remove based on zMSE score
@@ -285,13 +276,12 @@
 %         
 %     % 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);
@@ -321,7 +311,7 @@
 
 %re-calculate the sz variable
 sz=size(fids);
-    
+
 %FILLING IN DATA STRUCTURE
 out=in;
 out.fids=fids;
@@ -329,7 +319,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_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);