This uses fminspleas, downloadable from,
It requires an initial guess only for f0. The fit appears to be sensitive to the initial guess for r, whose final value is about 50% higher than your initial guess.
load data
flist={@(f0,f) 1./(f0^2 + f.^2)};
[f0,A0]=fminspleas(flist,200, f(4:end),pxx(4:end)) %first pass solve using fminspleas
r0=sqrt(A0*pi^2*b/(k*T)) %Recover r0 from A0
ft = fittype(@(A,f0,f) A./ (f0^2 + f.^2),...
'independent', {'f'}, 'dependent', {'pxx'});
result = fit(f(4:end),pxx(4:end),ft,'StartPoint', [A0,f0]) %now do formal fit
r=sqrt(result.A*pi^2*b/(k*T)) %Recover r from A
%plot(result,f(4:end),pxx(4:end), 'fit')
Yresults = result(f);
loglog(f,pxx, ".")
hold on
loglog(f,Yresults)
hold off
