Because you have a real-valued signal, the power spectral density is an even function of frequency. Therefore, there is no need to keep all 251 values in the PSD estimate. For an odd length input (251) if you keep the first round(251/2)+1 you have PSD estimates from 0 frequency (the first value) up to almost the Nyquist frequency.
Personally, I think a better decision would be to have kept floor(251/2)+1
If the length of the input is even (imagine here it was 250), then if you kept from 1 to 250/2+1, you would have PSD estimates from 0 frequency to the Nyquist.
The 1000 is the sampling frequency. The spacing between the DFT bins is Fs/N where Fs is the sampling frequency and N is the length. In this case 1000/251
Again, I would have done:
df = 1000/251; freq = 0:df:500;
So I would have done:
t = 0:.001:.25;
Fs = 1000;
x = sin(2*pi*50*t) + sin(2*pi*120*t);
y = x + 2*randn(size(t));
Y = fft(y,251);
Y = Y(1:floor(251/2)+1);
Y = 1/(251*Fs)*(Y.*conj(Y));
df = 1000/251;
freq = 0:df:500;
plot(freq,Y);
Here is an example for your situation:
t = 1:168;
x = cos((2*pi)/12*t)+randn(size(t));
% if you have the signal processing toolbox
[Pxx,F] = periodogram(x,rectwin(length(x)),length(x),1);
plot(F,10*log10(Pxx)); xlabel('Cycles/hour');
ylabel('dB/(Cycles/hour');
% if not
xdft = 1/168*fft(x);
xper = abs(xdft(1:length(x)/2+1)).^2;
df = 1/168;
freq = 0:df:1/2;
plot(freq,10*log10(xper))
xlabel('Cycles/hour');
ylabel('dB/Cycles/hour');
Note that MATLAB uses a convention of additionally scaling one-sided PSD estimates by 2 at frequencies beside 0 and the Fs/2, which I did not do with my fft() example.
