The OP posted the code causing the confusion in a duplicate question- I believe it was from here with the code duplicated below. A common confusion with this is the different scaling that will result for noise versus single tones. To complete an accurate power spectral density normalized per Hz that has both tones and noise, we need to use different scaling factors for each (for this reason as I explain in the last paragraph, I would not use a PSD to display such signals, but rather a spectrum plot and report the resolution bandwidth of the measurement). The reason is simply that the power measured in a pure single frequency tone will be independent of the resolution bandwidth of our measurement, while the power measured for distributed noise will go down as the resolution band width goes down (as in the latter case the effective filter of the measurement resolution bandwidth will select less of the total noise in the signal). In the DFT, the resolution bandwidth goes down as we increase the total number of bins, $N$.
That said, the scaling suggested in the link would be correct as $|X[k]|^2/(N f_s)$ for the case of noise that is spread across multiple bins with a rectangular windowed DFT, as I will explain in more detail below. (To a first order; the secondary affects of scalloping loss, spectral leakage, windowing loss etc are not detailed here but given in other posts referenced).
Let's start with the code from the link with scaling in question:
Fs = 1000;
t = 0:1/Fs:1-1/Fs;
x = cos(2*pi*100*t) + randn(size(t));
N = length(x);
xdft = fft(x);
xdft = xdft(1:N/2+1);
psdx = (1/(Fs*N)) * abs(xdft).^2;
psdx(2:end-1) = 2*psdx(2:end-1);
freq = 0:Fs/length(x):Fs/2;
plot(freq,10*log10(psdx))
grid on
title('Periodogram Using FFT')
xlabel('Frequency (Hz)')
ylabel('Power/Frequency (dB/Hz)')
Which results in the following plot:
This is the correct scaling for noise densities, but NOT for single tones. It only happens to work out in this case that Mathworks has posted since the length of the data is the same as the sampling rate! (Which may not have been a coincidence to avoid all this additional explanation).
For purposes of computing a PSD with both positive and negative frequencies, the DFT should be scaled by $1/N$ and then this result should be squared (as a complex conjugate product) to provide a normalized power quantity- to put it in units of per Hz as a power spectral density (applicable to noise only!) we divide that result by the resolution bandwidth of each bin which results in the scaling used. So the answer for scaling a DFT result $X[k]$ is determined as follows:
Normalized power per bin (I will use units of Watts for power which assumes the time domain signal is in units of volts across a normalized 1 ohm resistance):
$$(|X[k]|/N)^2 \space\space\space \text{Watts/bin}$$
Where $N$ is the total number of samples.
The resolution bandwidth (bandwidth of each DFT bin) of an unwindowed (rectangular window) DFT is simply the sampling rate divided by the number of bins:
$$\text{rbw} = \frac{f_s}{N}\space\space\space \text{Hz/bin}$$
Where $f_s$ is the sampling rate in Hz.
Thus the Power Spectral Density in units of Watts/Hz would be:
$$PSD = \frac{(|X[k]|/N)^2 \space\space\space \text{Watts/bin}}{\text{rbw}\space\space\space \text{Hz/bin}} = \frac{(|X[k]|/N)^2}{f_s/N}\space\space\space \text{Watts/Hz}$$
$$ = \frac{|X[k]|^2}{N f_s}$$
As given in the original MATLAB code.
The actual power units will be based on the units of the time domain waveform, so further scaling or conversion may be needed for different units of actual power, and for a one-sided spectrum the level is doubled as they do to account for the components in both the positive and negative frequencies (or left as is for a two-sided spectrum that shows both positive and negative frequencies). For relative power the above would be sufficient when converted to dB and compared to a reference such as total power in the waveform. As a caution, some implementations of the FFT may already include the scaling by $1/N$ as explained in further detail below.
In general, the scaling factor is given by what we want the vertical axis to represent. If we scale the DFT by 1/N as given by:
$$DFT = \frac{1}{N}\sum_{n=0}^{N-1}x[n]e^{j2\pi n k/N}$$
Then the vertical axis of the DFT will be scaled by the magnitude of each exponential component that each bin in the DFT result represents. For a specific example consider the sinusoidal function and it's representation as exponentials as given by Euler's identity:
$$\cos(\omega t) = \frac{1}{2}e^{j\omega t}+\frac{1}{2}e^{-j\omega t}$$
If we were to take the DFT of the sampled version of the above (and sampled at an integer multiple of $\omega$ such that we only get two non-zero bins in the DFT result), if we scaled the DFT by $1/N$ as given in the first formula introduced above, then the result for each DFT bin will have a magnitude of $1/2$ as given by the expression for the cosine above. Squaring this as a complex conjugate product provides a power unit per bin, but as explained earlier in the case of tones, no matter how much we reduce the resolution bandwidth we would expect the same power in each tone (since each tone occupies 0 bandwidth effectively). So in this case we would be done by simply computing the power of each tone as $(|X[k]|/N)^2$. For the case of noise signals, such as white noise where the power is spread evenly, then we can normalize the power in each bin to be a power/Hz quantity by dividing by the resolution bandwidth resulting in a scaling $(|X[k]|^2/(N f_s)$.
This is the case for a rectangular window. When using any other window the resolution bandwidth (or equivalent noise bandwidth) is increased as explained in this post. Additionally spectral tones versus noise would have a very different effect since a pure tone at bin center (for example) would occupy only one bin, while white noise (for example) as a density is spread uniformly across bins, and within a bin such as if we were to add more bins, the power in each bin would go down. In contrast with a single tone at bin center, if we add more bins, the power of the tone at bin center should not change (I say "bin center" to avoid getting into a scalloping loss and spectral leakage discussion).
Here is what we get if we repeat the MATLAB code for 100,000 samples with the same waveform using the same sampling rate. The noise comes out to be the same (with wider distribution given more samples but similar average power level/Hz) but the power level of the tone is now significantly higher as shown. This is because the power of the tone as is is given as $2[X[k]/N]^2$ and should not be modified by the resolution bandwidth per bin. If I were to repeat this with that scaling, the tones would be accurately represented but the noise would be significantly reduced since in this case the bandwidth of each bin is a sub-Hz quantity, so we need to scale that up by multiplying by $f_s/N$ as described above to put it in units of power/Hz.
What I typically do for spectrums that do have discrete stronger frequencies as well as noise is NOT normalize to a /Hz density but instead I compute the normalized power spectrum (using $(|X[k]|/N)^2$) and report on my plot the resolution bandwidth (as computed for whatever window was used), consistent with what we would see on any spectrum analyzer. Depending on if my interest is in an accurate measurement of the tones or the noise I will also scale for the window accordingly and report which quantity the plot was scaled for (but in this case the power difference between noise and tones is much less in the overall plot). However when I am creating plots of noise only (such as a phase noise plots) then in those cases I do give them as a density normalized per Hz as we would typically show.
