Sunday, May 14, 2017

signal processing - How is energy conserved for discrete Fourier transforms of different bin size?


I have gotten myself into a muddle regarding the discrete Fourier transform (DST) and would appreciate if someone could help me become unstuck.



I noticed that the fft of the following are the same:



  • A signal

  • The same signal but with zero padding

  • The same signal but repeated many times over and then normalised by the number of repetitions after fft.


All of these signals in frequency space have different bin sizes (df = fs/N). I understand the DFT as the coefficients of a given signal in the basis set of sinusoids who's periods are factors of the window length. By reducing df and therefore increasing the length of my sum how can energy by conserved if the coefficients of the sinusoids present in shorter sums are the same whilst still adding further sinusoids?


These observations only make sense to me if:



  • Matlab is normalising each bin by its width such that the fft instead of giving coefficients instead gives $\mathrm{Hz}^{-1}$ spectral levels.


  • Or my understanding of the DFT and the terminology 'bin' are entirely wrong and it instead produces samples of a $\mathrm{Hz}^{1}$ value.



Answer



Note: The location and/or order of applying the factor of $N$ in the following will vary between different software languages.



First, be careful with the normalizations. The vary from language-to-language not only in units but also in order of operations. For instance, Matlab adds the $N^{-1}$ normalization only on the inverse FFT whereas Mathematica applies a $N^{-1/2}$ normalization on both the forward and inverse transform. Some languages add the sample period to the normalization to make the units include $Hz^{-1}$ but most do not to allow the user freedom to input data in various time units other than seconds (technically, many FFT routines do not even ask for the time stamps associated with the array of data to be transformed).


Let us backup a little and define $f_{s}$ as the discrete sample rate (i.e., number of samples per second), $\tau$ the sample period (i.e., $f_{s}^{-1}$), $x_{n}$ the $n^{th}$ input array element, $X_{k} = FFT(x_{n})$, and $S_{k}$ the power spectral density, which we define as: $$ S_{k} = \frac{2 \ N}{f_{s}} \lvert X_{k} \rvert^{2} \text{; for k = 1,2,...,N/2 - 1} \tag{1} $$ Note that in most FFT algorithms operating on real data signals, the negative frequencies are mirror reflections of the positive except for the zeroth and N/2 frequencies which are not repeated (well I know of none without this feature but I figured I would be safe and say "most" instead of "all"). Meaning, the amplitude in the $f_{-n}$ bin should match that in the $f_{+n}$ bin, where $f_{n} = n/\tau$ for $n$ = -N/2, -N/2 + 1,...,N/2 - 2, N/2 - 1. Thus, Equation 1 is technically wrong for the zeroth and N/2 frequency bins, which should not include the factor of $2$.



You should be careful here as well because the normalization factor changes if you use zero-padding. For instance, suppose your input signal has $N$ points and you add $M$ zeros to make the input an even power of two (i.e., to optimize the computation time of the FFT). Then the new normalization factor would be the following: $$ S_{k} = \frac{2 \ \left( N + M \right)^{2}}{N \ f_{s}} \lvert X_{k} \rvert^{2} \text{; for k = 1,2,...,N/2} \tag{2} $$




A windowing function (Note that the article refereced below by Fredric Harris is an excellent reference for windowing functions.) is applied to reduce what are known as edge effects (e.g., discontinuities introduced by the implicit assumption of periodicity in the FFT algorithm). There are numerous different types but using one requires an additional normalization of which to be aware. Let us define the $j^{th}$ value of the windowing function as $w_{j}$, then we define the mean square value of the entire function as: $$ W_{ms} = \frac{1}{N} \sum_{j=0}^{N-1} \ w_{j}^{2} \tag{3} $$ Then we see that we need to divide $S_{k}$ by $W_{ms}$ to account for the power lost in $S_{k}$ by applying $w_{j}$. Since $W_{ms}$ is a scalar and scalars can commute, I tend to normalize my windowing functions by $W_{ms}$ prior to applying them to the input time series $x_{n}$.



Energy conservation can be derived from Parseval's relation, which shows that: $$ \sum_{n=0}^{N-1} \ \lvert x_{n} \rvert^{2} = \frac{1}{N} \sum_{k=0}^{N-1} \ \lvert X_{k} \rvert^{2} \tag{4} $$ where $X_{k} = FFT(x_{n})$ and $N$ is the number of elements in the input array.



By reducing df and therefore increasing the length of my sum how can energy by conserved if the coefficients of the sinusoids present in shorter sums are the same whilst still adding further sinusoids?



I am not sure I agree with this. The magnitude of and corresponding $f_{n}$ for the Fourier coefficients change when you change $N$ regardless of whether you use zero-padding or a proper windowing function. The reason zero-padding does not alter the energy is because you are adding zeros to a finite signal (i.e., FFTs expect finite numbers on input), thus zero amplitude signals. So the total power should not change if one properly normalizes $S_{k}$.



Or my understanding of the DFT and the terminology 'bin' are entirely wrong and it instead produces samples of a $Hz^{-1}$ value.




Again, be careful about units. Most FFT routines will not try to normalize by the sample period so there should be no $Hz^{-1}$ in $S_{k}$ until you apply one the normalizations I discussed above.




  • Donnelly, D. and B. Rust "The fast Fourier transform for experimentalists, Part I: Concepts," Comput. Sci. & Eng. 7, pp. 80--88, 2005.

  • Harris, F.J. "On the Use of Windows for Harmonic Analysis with the Discrete Fourier Transform," Proc. IEEE 66(1), pp. 51--83, doi:10.1109/PROC.1978.10837, 1978.

  • Paschmann, G. and P.W. Daly "Analysis Methods for Multi-Spacecraft Data," ISSI Scientific Reports Series 1, 1998.


No comments:

Post a Comment

classical mechanics - Moment of a force about a given axis (Torque) - Scalar or vectorial?

I am studying Statics and saw that: The moment of a force about a given axis (or Torque) is defined by the equation: $M_X = (\vec r \times \...