function [fp,fp_back,mina,fp_max] = fun_extraction(raw,rg,peak) II = reshape(raw,1,length(raw)); mina = floor(min(II)); maxa = ceil(max(II)); muhat = zeros(1,3); a = (mina:1:maxa); if isequal(rg,'rfp') gap1 = 10; gap2 = 20; else gap1 = 5; gap2 = 10; end [b,centers] = hist(II,a); b = b./(length(raw)); count_end = 1; count_end2 = 1; for i1= 1:length(b) if sum(b(1:i1))<0.9999 count_end = count_end+1; else break; end end for i1= 1:length(b) if sum(b(1:i1))<0.99 count_end2 = count_end2+1; else break; end end if count_end>length(centers) count_end = length(centers); end fp_max = centers(count_end); b(count_end) = b(count_end) + sum(b(count_end+1:end)); b_smooth = smooth(b,0.01)'; if peak==0 [peaks,locs] = findpeaks(b_smooth,'MinPeakHeight',0.001,'MinPeakDistance',10); else locs1 = find(centers==peak); if locs1>10 peak1 = max(b_smooth(locs1-10:locs1+10)); else peak1 = max(b_smooth(1:locs1+10)); end locs = find(b_smooth==peak1(1,1)); end if isempty(locs) [peaks,locs] = findpeaks(b_smooth,'MinPeakHeight',0.0005,'MinPeakDistance',10); end if isempty(locs) [peaks,locs] = findpeaks(b_smooth,'MinPeakHeight',0.0001,'MinPeakDistance',10); end if isempty(locs) [peaks,locs] = findpeaks(b_smooth); end if length(locs)>1 && locs(1,2)-locs(1,1)<10 locs_peaks = round((locs(1,1)+locs(1,2))./2); else locs_peaks = locs(1,1); end if peak == 0 mumin = centers(1,locs_peaks)-50; mumax = centers(1,locs_peaks)+50; else mumin = centers(1,locs_peaks); mumax = centers(1,locs_peaks); end xdata = centers(1:count_end2-1); ydata = b(1:count_end2-1); [xfit,resnorm] = lsqcurvefit(@fun_norm,[centers(1,locs_peaks),6.5,0.8],xdata,ydata,[mumin,0.1,0],[mumax,50,1]); muhat(1,1) = xfit(1,1); muhat(1,2) = 0; muhat(1,3) = 0; y =xfit(1,3).*normpdf(a,xfit(1,1),xfit(1,2)); y1 = xfit(1,3).*normpdf(a,xfit(1,1),xfit(1,2)); y2 = a.*0; y3 = a.*0; [xfit,resnorm] = lsqcurvefit(@fun_norm_2,[muhat(1,1),6.5,0.8,muhat(1,1)+gap1,6.5,0.2],xdata,ydata,[muhat(1,1)-10,0.1,0.01,muhat(1,1)+gap1-10,0.1,0.01],[muhat(1,1)+10,50,0.99,muhat(1,1)+500,50,0.99]); muhat(1,1) = xfit(1,1); muhat(1,2) = xfit(1,4); muhat(1,3) = 0; y =xfit(1,3).*normpdf(a,xfit(1,1),xfit(1,2))+xfit(1,6).*normpdf(a,xfit(1,4),xfit(1,5)); y1 = xfit(1,3).*normpdf(a,xfit(1,1),xfit(1,2)); y2 = xfit(1,6).*normpdf(a,xfit(1,4),xfit(1,5)); y3 = a.*0; if sum(b(1:count_end-1)-y(1:count_end-1))>0.01 [xfit,resnorm] = lsqcurvefit(@fun_norm_3,[muhat(1,1),6.5,0.8,muhat(1,2),6.5,0.2,muhat(1,1)+gap2,6.5,0.1],xdata,ydata,[muhat(1,1)-10,0.1,0.01,muhat(1,2)+gap1-10,0.1,0.01,muhat(1,1)+gap2-10,0.01,0],[muhat(1,1)+10,50,0.99,muhat(1,2)+10,50,0.99,mumax+1000,50,0.99]); muhat(1,1) = xfit(1,1); muhat(1,2) = xfit(1,4); muhat(1,3) = xfit(1,7); sigmahat = xfit(1,2); y =xfit(1,3).*normpdf(a,xfit(1,1),xfit(1,2))+xfit(1,6).*normpdf(a,xfit(1,4),xfit(1,5))+xfit(1,9).*normpdf(a,xfit(1,7),xfit(1,8)); y1 =xfit(1,3).*normpdf(a,xfit(1,1),xfit(1,2)); y2 =xfit(1,6).*normpdf(a,xfit(1,4),xfit(1,5)); y3 =xfit(1,9).*normpdf(a,xfit(1,7),xfit(1,8)); locs_peaks = round(muhat(1,1))-mina+1; x = centers-muhat(1,1)-9; end locs_peaks = round(muhat(1,1))-mina+1; fp(1,1) =sum(b(1,1:count_end).*(1:count_end))-(muhat(1,1)-mina+1).*xfit(1,3); fp_back = muhat(1,1)-mina+1; end