4
$\begingroup$

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?)

$\endgroup$
  • $\begingroup$ To me this is the best post that explained it all $\endgroup$ – Andre Feb 19 at 13:55
3
$\begingroup$

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 $$

$\endgroup$
1
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ I's like to point out that method of normalisation you proposed is correct only for the rectangular window. In case of Hamming window you need to divide by the sum of its samples: Y = fft(y) / length(win) (element-wise division is not necessary here). Please also note that for rectangular window it is also correct - sum(win) == length(y). $\endgroup$ – jojek Apr 14 '14 at 14:08
  • $\begingroup$ @jojek, surely. Thanks for pointing that out. $\endgroup$ – werediver Apr 14 '14 at 14:11

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.