PSD estimate using FFT - scaling and units

Question

I know that similar questions have already been asked – I have read the answers, yet I am still not sure if I understand the topic properly at some points.

I want to use FFT to calculate PSD estimate. FFT gives me a set of complex numbers that I want to transform into the PSD values. I am aware of the fact that only the first $N/2 + 1$ values are useful, and that the other half are complex conjugates of the first half.

1. I know that I will need to take the magnitude squared of each (useful) complex number the FFT gave me in order to find the PSD values. I also found out that I should incorporate some scaling/normalizing of the result. It is at this point that I am unsure. In some of the replies to similar questions it is said that the magnitude squared should be multiplied by $1/N$ only, while in others it is said that it should be multiplied by $1 / \left(N\cdot F_s\right)$. I have not been able to figure out which of these two should be used when. Can anyone explain?

2. The sampling frequency of my signal is $128\textrm{ Hz}$. The length of the signal to be transformed is $256$ samples. In this case I found out that the distance between the resultant values of the FFT (or PSD values) should be $128/256 = 0.5\textrm{ Hz}$. Is that correct?

3. The original signal is in $\mu V$. In what values will the result of the PSD be? Is it $\mu V^2$ or $\mu V^2 /\textrm{ Hz}$? (There might be a connection here with what I am asking in 1), am I right?)

Answer 1

Have a look at this paper, section 9 goes over how to scale or normalize the FFT output. For the PSD you need to divide by the effective noise bandwidth (ENBW) which is

$$ ENBW = f_s\frac{S_2}{S_1^{2}}\\ S_1 = \sum_{j=0}^{N-1}w_j\\ S_2 = \sum_{j=0}^{N-1}w_j^{2}\\ where~w_j~is~the~jth~sample~of~the~window~function $$

Answer 2

1. FFT (DFT) normalization

Your math software may not normalize the FFT output, indeed. In such a case you should either refer to the docs or check manualy which normalization is appropriate.

The following is applicable to FFT (DFT) with rectangular window as were noted by @jojek in the comments.

One way to check whether your normalization is correct is to calculate amplitude (not power) spectrum of say a A*sin(1:N) series. If you see a peak of height A, it's all right.

In Scilab I use the following function to get amplitude and power spectrum of a real signal:

function [SD] = spect(y, pwr)
    // Spectrum of a real signal

    if ~exists('pwr', 'l') then
        pwr = %f
    end

    if ~isreal(y) then
        error('Real signal expected.')
    end

    Y = fft(y) ./ length(y)
    SD = abs(Y(1 : $ / 2 + 1))
    SD(2 : $) = 2 .* SD(2 : $)
    if pwr then
        SD = (SD ./ sqrt(2)) .^ 2
    end
endfunction

2. FFT (DFT) frequency resolution

Resolution of FFT is $\frac{F_s}{2}$, so you're right here.

3. PSD Y-axis units

The Y-axis of the PSD plot is in power units, $\mu V^2$ in your case.