PSD looks to clean/continuous I generated an signal with a base frequency of $50\textrm{ Hz}$ and its harmonics and quantization noise in it (as an array in MATLAB). It's sampled fast enough to make sure there is no aliasing etc.

I calculate the PSD with a rectangular window in MATLAB. Due to the fact that the signal is synthesized in MATLAB I know where my harmonics are, so I chose length of my FFT to have a frequency resolution of $5\textrm{ Hz}$. I did this to be able to see my hamonics without the spectral leakage effect.

My problem is, that the PSD estimates differ from my expectations. I expected to see spikes caused by my harmonics and a kind of random noise floor. But what I get is kind of a continuous looking PSD with spikes as you can see in the image below:



The $x$-axis is the frequency in $\textrm{Hz}$ and the $y$-axis is in $\textrm{dBc/Hz}$ normalised to my base frequency power.

My Matlab code is:

NFFT = length(bitstream)-1; %df = fs / NFFT = 5 Hz
fVals = fs*(0:NFFT/2)/NFFT;

% calc fft
X_f = abs(fft(bitstream,NFFT));
% calculate periodogram (one sided)
pxx = ((X_f(1:NFFT/2+1))/NFFT).^2;
% plot
plot(fVals,10*log10(pxx/max(pxx)),'black');


The MATLAB function periodogram() as shown below (with rectangular window and same FFT length) produces the same output, so I thought my code is correct.

pxx = periodogram(bitstream,rectwin(length(bitstream)),NFFT);



Can you tell my why I see that?
Is the periodogram method not suitable for that case?
Is it caused by the rectangular window and/or the frequency resolution? If I use the built-in non-parametric estimation functions like pwelch(), periodogram(), pmtm() etc. without any further input it looks more noisy (as expected).
Or is it even better than expected because one of the frequencies in my simulation is the greatest common divider of the sampling frequency (using e.g. a signal with a different frequency like $49\textrm{ Hz}$ also leads to more noisy results)?
I believe this is due to your choice of using a rectangular window, which has the best frequency resolution (1 bin) but worst spectral leakage (dynamic range). The spectral leakage is of the form of a $\textrm{sinc}$ function (with nulls at $1/T$ where $T$ is your window length in time), which has an envelope that rolls off very slowly ($1/f$) in frequency. What I believe you see above is the cumulative effect of the $\textrm{sinc}$ function in frequency due to your window convolving with each tone. 

My favorite window choice is the Kaiser window kaiser(N,beta), as you can control the precision and dynamic range trade directly through the use of beta. Beta factors I typically use range from 4 to 12 but you can experiment with it and see what works best for you.  

Try this to see the window performance (MATLAB / Octave) in the frequency domain (showing the main lobe width (precision) and sidelobe level (dynamic range).

wind = kaiser(N, beta)
freqz(wind)