DFT for signals of different length and sampling rate I would like to compare oscillations of biological object. I have multiple datasets which are taken with different sampling rate and have different length. 

My current workflow:


Take derivative, because oscillations are overlapped over intense DC component.
Do zero padding to deal with different sample lengths and achieve the same bin size and, consequently, the same "digital" resolution in frequency. For this, I calculate the pad size: $$ pad = N - signalLength$$ ($N$ is padded signal length, $N=f_{sampling}/{\delta}f$, ${\delta}f$ is desired resolution in FFT frequency, $f_{sampling}$ is sampling frequency). Then I add vector of zeros of pad size to the end of signal: $$oneSample = [oneSample;zeros(pad,1)];$$
Take FFT: $$FFT = fft(oneSample,N);$$ 
Calculate PSD spectrum: $$PSD = (abs(FFT)^2)*normalizationCoefficient;$$ Output for a single sample looks like this:

Lastly, I look for maximum PSD value and corresponding FFT frequency for each sample and plot distributions.


Zero padding seems to give a reasonable results for frequencies, however I am still struggling with finding proper normalization coefficient $normalizationCoefficient$ for PDS spectrum. I started by simply dividing by sample length before zero padding: $$normalizationCoefficient = 1/signalLength$$ but PSD distributions are nonsense (sample, for which I am sure it has to have higher energy, has lower average maximum PSD, purple curve).



Then I tried $$normalizationCoefficient = f_{sampling}/signalLength$$ and results seem to make sense but I cannot find anywhere the justification for this.

 

So, my question is: what normalization coefficient for PSD should I use in my case? 

Wikipedia suggests $$normalizationCoefficient = 1/(f_{sampling}^2*signalLength)$$  



but the distributions do not make sense again:



Could anyone point me to the right direction, please?

UPD: on use of zero padding. It does not help to improve the actual resolution in frequency domain, however I find it benificial for the latest stage where I use Matlab's findpeaks function to determine the peak with max amplitude.
Fft normalization is often an overlooked topic and many articles dealing with FFT just say to divide the fft ouput by N (divide by either N or N/2 depending on the actual
FFT algorithm used). 
Real FFT  usually discards the upper part of the spectrum (due to FFT redundancy) so the normalization coefficient here is N/2, rather than N.


  So, my question is: what normalization coefficient for PSD should I
  use in my case?


I prefer the approach taken by Heinzel et al, and according to them, the normalized/scaled PSD  can be computed as

PSD = ( 2 * (yi * yi) ) / ( fs * S2 ); // i = 0…N/2  (1)

S2 = sum(wi * wi); //i = 0...N (2)

(yi is the complex result of the fft at bin i and wi is the value of the window function at point i )


The sum S2 is used instead of N in order to better estimate the scaling
factor, given a particular window (rectangular, Hann etc). Using no windows at all corresponds to a rectangular window and in the case of a rectangular window S2 is equal to N.

Heinzel et al also define a metric called ENBW (equivalent noise bandwidth), which is at the heart of the here defined scaling technique and for those interested, complete derivations (including 1 and 2) can be found here https://holometer.fnal.gov/GH_FFT.pdf

Note also that for a rectangular window, ENBW is equal to fsampling / N, which may explain why the results started to make sense for you once this normalization factor was applied.