Hi DK,
Since you are after power spectral density, and not the spectrum of the linear quantity fft(x), this means the square of fft(x) is correct. If you think of x as analogous to a voltage, then power is proportional to x^2, and squaring survives in the psd. The only remaining question is normalization, 1/(2*pi*n) in psdx above.
For power per unit bandwidth in Hz, which I think is what is usually wanted, the factor of 2*pi would NOT be there, and neither would it be in the expression for freq.
It pays to take a look at the fft as an approximation to the continuous fourier transform,
y(f) = Integral x(t) e^(-2*pi*i*f*t) dt
In terms of an n-point fft with a corresponding time array of spacing delt, this is
y = fft(x)*delt
where the fft does the summing and delt takes the role of dt.
The average power of an n-point x array is
Pave = (1/n)*sum(x.^2)
and the idea is to come up with
psd = dPave/df such that Integral psd(f) df = Pave
The approximation using the fft and an n-point frequency array of proper spacing delf is
sum(psd)*delf = Pave
The psd can also be done in terms of circular frequency w = 2*pi*f and a circ frequency array of spacing delw so that
psd = dPave/dw sum(psd)*delw = Pave
These two psds are related by a factor of 2*pi, and are psd1 and psd3 below. For more of a signal processing context one can plug in delt = 1. Then by the fft golden rule, delf = 1/n, which is the normalized frequency spacing. That yields psd2 below. delw = 2*pi/n is the normalized circular frequency and yields psd4 below, which differs from pds2 by a factor of 2*pi. Your code at present uses psd4, so the frequency grid would have to be in terms of w.
Your code also tosses out the negative frequencies and doubles the result for the positive frequencies, which (very nearly) doubles the psd and means you sum over n/2 points, but I did not do that since it does not change the 2*pi-or-not argument.
n = 1000;
delt = .04;
delf = 1/(n*delt) % delt*delf = 1/n fft golden rule
delw = 2*pi*delf
x = 2*rand(1,n)-1;
Pave = (1/n)*sum(x.^2);
% by analogy to the continuous case
y1 = fft(x)*delt;
psd1 = abs(y1).^2/(n*delt); % psd, dPave/df
Pave1 = sum(psd1)*delf;
% use delt = 1, normalized f
y2 = fft(x);
psd2 = abs(y2).^2/n;
Pave2 = sum(psd2)*(1/n);
% w instead of f, continuous
y3 = fft(x)*delt;
psd3 = abs(y3).^2/((2*pi)*(n*delt)); % psd, dPave/dw
Pave3 = sum(psd3)*delw;
% delt = 1, normalized w
y4 = fft(x);
psd4 = abs(y4).^2/((2*pi)*n);
Pave4 = sum(psd4)*(2*pi/n);
[Pave Pave1 Pave2 Pave3 Pave4] % they all agree
