I got a time series data. t = array (of regular interval starting from 0) length of t = n = 2080 dt = temporal spacing = 2.e-10 I have a set of data at t values E(t) = array

Now, since I want to covert it to fourier space,

f = sample frequency = array = (np.fft.rfftfreq(n2,dt))

E_f = E(f) = (np.fft.rfft(E))/n

so, here my concern is do I at all need to divide by "n" as mention in above equation ?

Now, I want to estimate total power in between 30 MHz to 100 MHz

how shall I proceed ?


In most implementations the FFT returns the following for the DFT:

$$X[k] = \sum_{n=0}^{N-1}x[n]e^{-j 2\pi nk/N}$$

Which would result in the total output being scaled by $N$. The FFT implementation may or may not be scaled by N so it is always good to check the documentation directly for the fft that is being used. The total power is the sum square of all the FFT bins for a full FFT of length N where the time series data is of length N, once proper scaling is accounted for an assuming no windowing beyond the implied rectangular window is used. (I refer to another post that covers how to properly handle the windowing for power estimation). So for the result as given above, the total signal power would be (and assuming units of voltage or current in the time domain, with a normalized resistance of 1 then the units of the output would be watts):

$$P = \sum_k\bigg(\frac{X[k]}{N}\bigg)^2$$

Besides reading the documentation, a quick test is to do fft([1, 1, 1, 1]) where fft() refers to whatever fft function is used; if the function needs to be scaled by $N$ before computing the total power, the result for above would be [4,0,0,0]. If the result was [1,0,0,0] then the DFT would have already been scaled by $N$ or another option is the FFT is scaled by $1/\sqrt{N}$ so that the fft and ifft are symmetric operations and the result would be [2,0,0,0].

Tne numpy fft function the OP is using rfft only returns the positive half spectrum, which would be the complex conjugate of the negative half spectrum for real signals and handles odd and even $N$ differently in terms of how many samples are included. For a consistent result, I recommend using fft which will return $N$ samples scaled by $N$, in which case to estimate the power you would first divide the fft results by $N$ and then sum the square of all the bins (consistent with Parseval's Theorem).

Also for Python it is recommended to use the fft in scipy.fftpack and not numpy (which is now in both for reverse compatibility reasons but is expected to eventually be removed from numpy; numpy is to be used for the lower level vector processing and broadcasting operations that scipy depends on but not the higher level numerical processing algorithms that are now all in scipy).

If windowing is applied, then the power computation needs to consider the coherent and non-coherent gain of the window. This is straightforward for a single tone whose power spectral density only occupies one bin; it will be reduced by the coherent gain of the window which is given by the sum of the window coefficients and scaled by $N$ (which is always a loss). However for signals with spectral content that spans multiple bins, we must be careful to consider that the bandwidth of a windowed DFT bin will always span multiple bins, so after the decrease due to coherent gain, there will be an increase due to the double counting as all the bins are summed in the power computation previously given. The net result is a processing gain that is applicable to signals with spectral densities spanning multiple bins and described further in the links below.

Find the Equivalent Noise Bandwidth

Effect of windowing on noise

  • $\begingroup$ So, in case of my time series data, is the way I build my frequency array "f" from dt, is it the correct way of doing ?( Provided I have used scipy.fftpack) and just performing fft(E) will it provide E_f corresponding to the particular frequency inhabitants of that "freq" array ? Like in case of a function I have implemented fft quite precisely but got doubtful for this set of some discrete time series data. So my ultimate concern is whether the way I created "f" and "E_f" data is correct ? Such that we can claim that, "E_f[i]" is the fourier transform value at "f[i] ? $\endgroup$ Apr 19 '21 at 19:39
  • $\begingroup$ @SubhadipSaha is it important to you to use rfft or would you be fine in using fft? For fft the frequency axis for n = 0 to N-1 is n*fs/N $\endgroup$ Apr 19 '21 at 21:04
  • $\begingroup$ Yes, for me rfft() shall be appropriate. And for that, the way frequency array is constructed from dt, is that correct ? I understand it would be from 0 to (N+1)/2. But, will E_f[i] be get assigned to f[i] automatically ? Now, is "f" frequency or angular frequency ? Since Fourier space of time is angular frequency space. Now to get power, sum over (E_f)^2(df) over the desired bandwidth is sufficient? Or any further factor need to be divided to normalise power ? Since dt is 2e-10 sec. $\endgroup$ Apr 20 '21 at 11:15
  • $\begingroup$ I would recommend using fft (simpler as you don’t have to treat even and odd datasets differently). If dt is in units of seconds then the frequency axis would be in units of Hz. (The N fft bins correspond to N being the sampling rate and go over 0 to N-1 thus my last comment. Starting from that, the rest of my answer should make sense. $\endgroup$ Apr 20 '21 at 11:51

Your Answer

