function [ve, Npeak, width, pos, contrast, width_err, pos_err, contrast_err] = esrfit_xN(f,data)

if find(isnan(data))
    ve = [0,0,0,0,0,0,0];
    Npeak = 0;
    width = 0;
    pos = 0;
    contrast = 0;
    return;
end

%smooting the curve
n = length(data)/25; %averging over n neighbors
datasmth = smooth(data,n)';

%find peaks (position, contrast) using addon "peakfinder"
[peakLoc,peakMag] = peakfinder(datasmth(1,:),'', 1.5*mean(datasmth(1,:)), -1, false);

%Defining init parameters for fit
Npeak = length(peakLoc);
bkg = max(data); %background if normalized close to one
pos = f(peakLoc); %peakLoc*(f(end)-f(1))/length(f)+f(1); %comes from peakfinder
amp = bkg - peakMag; %comes from peakfinder
width = 0.1*10^-3*ones(1,Npeak); %(MHz) not a very elegant way, need to think about sth more clever

%compile init parameters for Lorentz fitting
if Npeak>1
    initPar = vertcat(amp,pos,width); %concatenate
    initPar = initPar(1:end); %make line vector
    initPar = [bkg,initPar];  %add bkg
else
    initPar = vertcat(amp,pos,width); %concatenate
    initPar = initPar(1:end); %make line vector
    initPar = [bkg,initPar'];  %add bkg
end

%fit
vs = initPar; 
fct = @(v,x) Lorentz_xN_Func(Npeak, v, x);
[ve,r,J,COVB,mse] = nlinfit(f,data,fct,vs);

%extract error bars
errfitres_abs = nlparci(ve,r,'jacobian',J); %absolute value boundaries
errfitres_diff = abs((errfitres_abs(:,1) - errfitres_abs(:,2))/2); %difference values --> error bar 

bkg = ve(1);
bkg_err = errfitres_diff(1);

for i=1:Npeak
   width(i) = 1000*ve(4+(i-1)*3); %width in MHz
   width_err(i) = 1000*errfitres_diff(4+(i-1)*3); %error bar in MHz
   
   pos(i) = ve(3+(i-1)*3);  %peak position in GHz
   pos_err(i) = errfitres_diff(3+(i-1)*3); %error bar in MHz
   
   contrast(i) = 100*ve(2+(i-1)*3)/bkg; %contrast in %
   contrast_err(i) = 100*(errfitres_diff(2+(i-1)*3)/bkg + ve(2+(i-1)*3)/bkg^2*bkg_err); %error bar in MHz
end

end