function [est_len, quality] = single_interval_estimator(win_sig, win)


            %% Remove mean
            win_sig = bsxfun(@minus,win_sig,mean(win_sig,1));
                   
            %% Compute Scores
            combined = ones(size(win));
            probabilities = zeros(size(win_sig,2), length(win));
            
            for k = 1:size(win_sig,2)

%                xc = x_spectrum_mid_noscale(win_sig(:,k), win(1), 0);
%                ms = max_spectrum_mid(win_sig(:,k), win(1));
%                ad = admf_mid(win_sig(:,k), win(1), 0);
                
                
                
                %disp(win_sig(1:6,k));
                
                [xc, ms, ad] = xc_ms_ad(win_sig(:,k), win(1), 0);
                
                %disp(xc);
                
                xc = xc - min(xc);
                xc = xc/(sum(xc)+eps)*0.9;
                xc = xc + 0.1/length(xc);
                
                ms = ms - min(ms);
                ms = ms/(sum(ms)+eps)*0.9;
                ms = ms + 0.1/length(ms);
                
                ad = 1./(ad+eps);
                ad = ad - min(ad);
                ad = ad/(sum(ad)+eps)*0.9;
                ad = ad + 0.1/length(ad);

                probabilities(k,:) = xc.*ad.*ms;
                
                [~, extrs, ~, ampls] = peakdetector(probabilities(k,:));
                %[ampls, extrs] = findpeaks(probabilities(k,:));

                
                %[ampls, extrs] = max(probabilities(k,:));
                
                [max_ampl, idx] = max(ampls);
                if isempty(max_ampl)
                    max_ampl = 0;
                end
                
                if length(ampls) > 1
                    [ampls_sorted, sort_idx] = sort(ampls);
                    extrs_sorted = extrs(sort_idx);
                    
                    if (ampls_sorted(end)-ampls_sorted(end-1))/ampls_sorted(end) < 0.999 ...
                            && abs((extrs_sorted(end)+win(1)-1)/2 - (extrs_sorted(end-1)+win(1)-1)) < 5
                        max_ampl = ampls_sorted(end-1);
                        idx = find(extrs == extrs_sorted(end-1));
                    end
                end
                
                quality(k) = max_ampl/sum(probabilities(k,:));

                if isempty(extrs)
                    %If no peak exists, return the last valid interval length
                    %This should rarely happen (if ever)
                    est_len(k) = 0;
                else
%                     peak = extrs(idx);
%                     if peak > 2 && peak < length(probabilities(k,:))-1
%                         %If the peak is not right at the edge of the window
%                         %interpolate using a 2nd order polynomial to achieve
%                         %sub-sample accuracy
%                         poly = polyfit(-1:1, probabilities(k,peak-1:peak+1), 2);
%                         offset = -0.5 * poly(2)/poly(1);
%                     else
%                         offset = 0;
%                     end
                    est_len(k) = round(win(1)+extrs(idx)-1);%+offset;
                end
                
                
                combined = combined .* probabilities(k,:);

            end
              
            %Detect peaks in combined score and select largest peak
%             tic
            [~, extrs, ~, ampls] = peakdetector(combined);
%              toc
            
%            [ampls, extrs] = findpeaks(combined);
            
            %[ampls, extrs] = max(combined);
            [max_ampl, idx] = max(ampls);
            if isempty(max_ampl)
                max_ampl = 0;
            end
            
            if length(ampls) > 1
                    [ampls_sorted, sort_idx] = sort(ampls);
                    extrs_sorted = extrs(sort_idx);
                    
                    if (ampls_sorted(end)-ampls_sorted(end-1))/ampls_sorted(end) < 0.9999 ...
                            && abs((extrs_sorted(end)+win(1)-1)/2 - (extrs_sorted(end-1)+win(1)-1)) < 5
                        max_ampl = ampls_sorted(end-1);
                        idx = find(extrs == extrs_sorted(end-1));
                    end
            end
            
            quality(size(win_sig,2)+1) = max_ampl/sum(combined);
            
            if isempty(extrs)
                %If no peak exists, return the last valid interval length
                %This should rarely happen (if ever)
                est_len(size(win_sig,2)+1) = 0;
            else
%                 peak = extrs(idx);
%                 if peak > 2 && peak < length(combined)-1
%                     %If the peak is not right at the edge of the window
%                     %interpolate using a 2nd order polynomial to achieve
%                     %sub-sample accuracy
%                     poly = polyfit(-1:1, combined(peak-1:peak+1), 2);
%                     offset = -0.5 * poly(2)/poly(1);
%                 else
%                     offset = 0;
%                 end

                est_len(size(win_sig,2)+1) = round(win(1)+extrs(idx)-1);%+offset;
            end
            
end


