Only to be complete

here is the algorithm (matlab) for 2N+1 points

Code:

N=3;
io=find(z==max(z),1); % locate peak
zm=mean(diff(z(io-N:io))); % average negative differences (i=i0-N to io)
zp=mean(diff(z(io:io+N))); % average positive differences (i=io to io+N)
b=1/(2*N)*(zp-zm);
dxo=-N/2*(zp+zm)/(zp-zm);
a=z(io)-b*dxo.^2;
xo=x(io)+dxo;
yo=a;

Having more than 3 points for fitting a quadratic function, i.e. N>1, is useful if data are noisy.

The formula is partially as done by @frankzappa. it only includes the linear interpolation and estimates also the modelled peak of the quadratic.

Edit:

obviously

Code:

zm=mean(diff(z(io-N:io))); % average negative differences (i=i0-N to io)
zp=mean(diff(z(io:io+N))); % average positive differences (i=io to io+N)

is equivalent to

Code:

zm=(z(io)-z(io-N))/N;
zp=(z(io+N)-z(io))/N;

which proves that averaging does not help with noise (in fact it only undoes the differentiation noise)