3
$\begingroup$

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:

  1. Take derivative, because oscillations are overlapped over intense DC component.
  2. 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)];$$
  3. Take FFT: $$FFT = fft(oneSample,N);$$
  4. Calculate PSD spectrum: $$PSD = (abs(FFT)^2)*normalizationCoefficient;$$ Output for a single sample looks like this: (PSD(f) and signal(t))
  5. 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).

enter image description here

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

enter image description here

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

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

enter image description here

but the distributions do not make sense again:

enter image description here

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.

enter image description here

Improve this question
$\endgroup$
  • $\begingroup$ You mentioned your $f_s$ are different in the beginning. If the $N$ is same for all, the the frequency resolution will not be same for each set. They will be $f_{s1}/N$ for set 1, $f_{s2}/N$ for set 2 etc. So you will not be comparing pdf with the same frequency resolution. $\endgroup$ – jithin Apr 3 '20 at 4:11
  • $\begingroup$ I actually adjust N by zero padding in order to have the same bin size, $N = f_{sampling}/{\delta}f$, where I choose ${\delta}f$. So $N$ is different for different sampling rates. $\endgroup$ – muxika Apr 3 '20 at 19:46
2
$\begingroup$

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.

Improve this answer
$\endgroup$
  • $\begingroup$ This is exactly what I was looking for, thanks! $\endgroup$ – muxika Apr 3 '20 at 19:02
1
$\begingroup$

I would recommend a different approach without zero padding.

Choose your DFT frame sizes based a common time duration and do different sized DFTs on your different signals with a 1/N (or 2/N if you prefer) normalization factor. The rationale for this is that the bins in a DFT correspond to frequencies in units of cycles per frame making the same bin index all your frames represents the same frequency. With those normalizations the magnitudes in the bins will also be comparable.

For pure tones that are whole cycles per frame (isolate to a conjugate bin pair in the DFT), the values from your different sampling rates will match exactly. For pure tones between the bins, the "leakage" will be a very close match. So, for a mix of tones (which every signal can be expressed as) the two DFTs will be directly comparable.

The only difference will be that your more densely sampled difference will have DFT higher frequency bins that are beyond the range of the less densely sampled one.

Of course, you will want to low pass filter for improved results as you have already found out. One of my favorite techniques is to use the Differ with a 0.5 coefficient as described here:

This gives you differentiaion and smoothing at the same time, but the end points are a little distorted. Therefore your smoothing frame should be larger and you should center your DFT frame within that. Multiple smoothings will work if needed and you can either keep the smoothed or differed signal for the next step.

How much overlapping you want to do is up to you. It is sort of a sampling rate question of its own.

Improve this answer
$\endgroup$
  • $\begingroup$ Thanks a lot for interesting suggestion! IF I understood you correctly, you suggest to apply exponential smoother to the samples with higher sampling rate, and then take FFT with different bin sizes. For the comparison of datasets, I wonder why not just use frequencies from the Nyquist interval for the dataset with lowest sampling rate $f_{Nyquist} = min(f_{sampling})/2$? Exponential smoother that you suggest introduces frequency dependent amplitude dampening as mentioned in the link which you are referring to. $\endgroup$ – muxika Apr 3 '20 at 19:43
  • $\begingroup$ @muxika The DFT itself does quite a bit of smoothing for you so you should start with that as baseline, then try different techniques and compare. Any pure DC component gets relegated to bin 0, which fixes that problem for you. To make spectograms I recommend using a VonHann window. To make parameter estimations, no window and the formulas in my blogs. Many other techniques exist as well. The key concept is that the bin index corresponds to the frequency in units of cycles per frame. Therefore your comparisons will only be universally valid up to the lowest Nyquist with no aliases. $\endgroup$ – Cedron Dawg Apr 3 '20 at 20:16
1
$\begingroup$

Zero padding a signal in time domain will not increase resolution in the frequency domain. Consider N< M, It may seem that increase in length by zero padding (adding M-N zeros) increases the resolution, but it simply is observing the original N DFT coefficients interpolated by lagrangian interpolation and now observing this interpolation at M equally spaced points. It is not the same as resolution of an M point time domain signal with actual signal values at those M points

Improve this answer
$\endgroup$
  • $\begingroup$ I realize that zero padding does not increase an actual resolution in the frequency domain, this is why I used term "digital" resolution in my question. However, I find it beneficial to use zero padding, because to compare different datasets I am taking frequency with max PSD and plot the distributions. For this, I use Matlab's findpeaks function. And it seems to work better on zero-padded samples. $\endgroup$ – muxika Apr 3 '20 at 19:09

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.