At this point I believe that the documentation is misleading in that it appears to imply Power/unit_frequency but the result is in Power/radian, smaller by a factor of 1/2pi (but doubled for a one-sided spectrum). At least, that seems the most likely possibility. pwelch does not seem to be accessable code, so here is an example for which I constructed a function called pwelch1 to look at the scaling.
Here the normalized f runs from -1/2 to 1/2. Integrating dS/df over all frequencies should give the average power S, which it does. Both are right around 1/3, the theoretical value. The function does contain the Hamming function correction factor U which is approximately 0.4, depending on the number of points.
pwelch1 always averages 8 windows with 50% overlap and there is no attempt to use ffts with length of 2^n since I think that effort is usually a waste of time.
In the Matlab code, w runs from 0 to pi, which certainly looks like normaliazed circular frequency. The code outputs positive frequencies only and multiplies the spectrum by A = 2 to make up for that. For Power/radian instead of Power/unit_frequecy, pwelch is smaller than pwlech1 by a factor of A/(2*pi) = 1/pi as you can see on the plot.
It would be interesting to hear what Mathworks might say on this toplic.
function [Z f] = pwelch1(y)
U = (1/nwin)*trapz(h.^2);
indstart(2*k-1) = 1+(k-1)*nwin;
indstart(10-2*k) = N+1-k*nwin;
ind = indstart(k)+(0:(nwin-1));
Z = Z + abs(fft(y(ind).*h)).^2;
Z = fftshift(Z)/(nwin*8*U);
f = (ceil(-nwin/2):ceil(nwin/2)-1)/nwin;