Merge pull request #381 from schorschinho/develop

Update SpecReg functions

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,15 +298,15 @@
 
 %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;

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,27 +178,27 @@
         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
@@ -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,15 +311,15 @@
 
 %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;