The example is single-sided from F=0:Nyquist (Fs/2). Since the FFT returned by Matlab is double-sided, from -Fmax:+Fmax,
f = Fs/2*linspace(0,NFFT/2+1);
as the example uses is the positive frequency vector. NFFT/2+1 is the location of the Nyquist frequency component in the returned FFT; (+)ive frequencies start from the origin of the array df apart positive with higher index to Nyquist, then negative with -Fs at the next element past NFFT/2+1 --> NFFT/2+2, the midpoint plus one being -(Fs/2-df).
The 1/L normalization is the normalization of the output for the length of the input signal; in the example there are L actual data points with power associated with them but NFFT points in the FFT--the normalization to return the amplitude of the PSD to match the input is only for the number of elements in the vector with data; not including the additional NFFT-L points that are just augmented zeros.
To do the two-sided plot, just change the range for f and the subscripting for the FFT data--
I get the following figure for this case--
Again, read the code for the example at FFT very carefully to see the difference between L and NFFT...the former is the length of the actual time series (1000 in the example, 1500 in your case) whereas the latter ( NFFT ) is the length of the transform which was determined using nextpow2 to get the next higher power-of-two length series for the (slight) efficiency in computation that makes internally. Thus NFFT would be 1024 in the example but 2048 for your case as 1500>1024 so must augment to next.
But, for proper normalization, only the length of the original signal L is used as all the additional points are simply introduced zeros and so to divide by NFFT where NFFT > L will reduce the apparent power by that ratio.
If you use the default argument for size in the call to FFT, then the routine will return the same length transform as the input time series and you should set
to get proper normalization.
Oh, the factor of two in the plot is since Matlab FFT returns double-sided, half the input energy is in each the positive and negative frequencies...strictly speaking the factor there should be included only for the one-sided plot but I kept it just for presentation purposes to make the actual peak match the input time series amplitude -- in reality, it should be unity for the two-sided so the sum of the two peaks is the total.