2
$\begingroup$

One thing that has always puzzled me is making the units of an analytically derived Power Spectrum Density (PSD) consistent with the units of an FFT.

Let's say we record the output of a frequency generator, which we just set to output a signal of the form $$v(t) = V_{\rm{pk}} \cos (2 \pi f_0 t)$$ for some finite amount of time, $T$, and take the FFT of this acquired signal. We might then expect to see something like this: enter image description here

Now as I understand it the FFT is a Linear Amplitude Spectrum (LAS) which will have units of Volts, or, a Power Spectrum (PS) which will have units of Volts-squared, depending on what we choose do with the resultant FFT. But it is clear the FFT can only have units of volts because that is what the units of the time-transient signal are.

Now lets say I want to fit this FFT with some realistic function. This is straight forward as I have the transient signal and I can calculate, analytically a PSD from a time truncated Fourier transform: $$\frac{1}{T}\left|\int_{0}^{T} v(t) e^{-2 \pi j f t} \ {\rm{d}}t \right|^{2}$$ which will be in units of Volts-squared per Hertz ($\rm{V^{2}/Hz}$).

So its clear that the units of the FFT and the units of my function I want to fit to it don't match up! Multiplying by some quantity in units of $\rm{Hz}$ fixes the problem for example either the bandwidth of the FFT, $\Delta f = 1/ T$ or the Effective Noise Bandwidth (ENBW).

How can i reconcile the units of my analytic expression and the units of my FFT?


The objective is to simply fit my FFT data, so I given that all this should only really affect the amplitude, I could of course allow all this to be absorbed by an arbitrary constant -- but it would be nicer to have a consistent result.

$\endgroup$
1
+100
$\begingroup$

I'm a big believer in this form of the continuous Fourier Transform and inverse:

$$ X(f) \triangleq \mathscr{F}\Big\{x(t)\Big\} = \int\limits_{-\infty}^{+\infty} x(t) \ e^{-j 2 \pi f t} \ \mathrm{d}t $$

$$ x(t) \triangleq \mathscr{F}^{-1}\Big\{X(f)\Big\} = \int\limits_{-\infty}^{+\infty} X(f) \ e^{+j 2 \pi f t} \ \mathrm{d}f $$

because I like the symmetry between the two reciprocal domains.

Let $x(t)$ be a finite-power signal as opposed to a finite energy signal. The power of $x(t)$ is

$$\begin{align} \overline{x^2} &= \ \lim_{T \to +\infty} \frac{1}{T}\int_{-\frac{T}2}^{\frac{T}2} \Big|x(t)\Big|^2 \ \mathrm{d}t \\ &= \ \lim_{T \to +\infty} \frac{1}{T}\int_{-\infty}^{\infty} \Big|x_T(t)\Big|^2 \ \mathrm{d}t \\ \end{align}$$

where $x_T(t)$ is the finite-energy signal defined as identical to $x(t)$ within a finite segment of time:

$$ x_T(t) \triangleq \begin{cases} x(t) \qquad & |t| < \frac{T}2 \\ \\ 0 \qquad & |t| > \frac{T}2 \\ \end{cases} $$

Now, fix $T$ to be something large and positive. Parseval's theorem tells us that the energy integral has an equivalent in the frequency domain:

$$ \int_{-\infty}^{\infty} \Big|x_T(t)\Big|^2 \ \mathrm{d}t = \int_{-\infty}^{\infty} \Big|X_T(f)\Big|^2 \ \mathrm{d}f$$

where $ X_T(f) \triangleq \mathscr{F}\Big\{x_T(t)\Big\}$.

Now let's pretend that positive frequencies and negative frequencies are different (and they are for the complex exponential, $e^{j2\pi ft}$), then if $x_T(t)$ is passed through and came out an ideal brickwall filter with a skinny bandwidth $B>0$ and centered at frequency $f_0$, then:

$$ X_T(f) \approx \begin{cases} X_T(f_0) \qquad & |f-f_0| < \frac{B}2 \\ \\ 0 \qquad & |f-f_0| > \frac{B}2 \\ \end{cases} $$

and that energy integral would be proportional to bandwidth, $B$:

$$\begin{align} \int_{-\infty}^{\infty} |x_T(t)|^2 \ \mathrm{d}t &= \int_{-\infty}^{\infty} \Big|X_T(f)\Big|^2 \ \mathrm{d}f \\ &\approx \int_{f_0-\frac{B}2}^{f_0+\frac{B}2} \Big|X_T(f_0)\Big|^2 \ \mathrm{d}f \\ &= \Big|X_T(f_0)\Big|^2 B\\ \end{align}$$

Now that is the energy in a segment of frequency, centered at $f_0$ with a bandwidth of $B$. This energy is expended over a time of width $T$, so the mean power over that time is

$$ \tfrac{1}T \Big|X_T(f_0)\Big|^2 B $$

which is proportional to the bandwidth, $B$, so the power per unit frequency around frequency $f_0$ is what multiplies the bandwidth, $B$, which is $\frac{1}T |X_T(f_0)|^2$ in the vicinity of frequency $f_0$.

If $x(t)$ were in volts and $B$ were in Hz, then $\frac{1}T |X_T(f)|^2$ would be "volts² per Hz" in the vicinity of frequency $f$. So to get the power over all frequencies you would add up (or integrate) all of the power components for all frequencies (negative and positive) and have:

$$\begin{align} \frac{1}T \int_{-\infty}^{\infty} \Big|X_T(f)\Big|^2 \ \mathrm{d}f &= \frac{1}T \int_{-\infty}^{\infty} \Big|x_T(t)\Big|^2 \ \mathrm{d}t \\ &= \frac{1}T \int_{-\frac{T}2}^{\frac{T}2} \Big|x(t)\Big|^2 \ \mathrm{d}t \\ \end{align} $$

Now that's for a large, but finite $T$. Note I am going with $-\frac{T}2<t<\frac{T}2$ instead of $0<t<T$.

Now that's the first half (which confirms we need to keep the $\frac{1}T$). The second half of the problem is expressing the integral as a Riemann sum and relating that to the DFT.

Now, if your sample rate is $f_\mathrm{s}$, that means your sampling period is $\frac{1}{f_\mathrm{s}}$ and Nyquist is $\frac{f_\mathrm{s}}2$. If $x_T(t)$ is sampled at rate $f_\mathrm{s}$, there should be no energy in the spectrum $X_T(f)$ at frequencies having magnitude above Nyquist. Now, it turns that that theoretically, $x_T(t)$ cannot be both time-limited and band-limited at the same time, but if we make the limits high enough, it's good enough for illustration.

$$\begin{align} X_T(f) \triangleq \mathscr{F}\Big\{x_T(t)\Big\} &= \int\limits_{-\infty}^{+\infty} x_T(t) \ e^{-j 2 \pi f t} \ \mathrm{d}t \\ X(f) &\approx \int\limits_{-\frac{T}2}^{+\frac{T}2} x(t) \ e^{-j 2 \pi f t} \ \mathrm{d}t \\ \end{align}$$

$$\begin{align} x_T(t) \triangleq \mathscr{F}^{-1}\Big\{X_T(f)\Big\} &= \int\limits_{-\infty}^{+\infty} X_T(f) \ e^{+j 2 \pi f t} \ \mathrm{d}f \\ x(t) &\approx \int\limits_{-\frac{f_\mathrm{s}}{2}}^{+\frac{f_\mathrm{s}}{2}} X(f) \ e^{+j 2 \pi f t} \ \mathrm{d}f \\ \end{align}$$


Now the form of Riemann summation with equal-width rectangles is

$$ \int\limits_a^b f(x) \ \mathrm{d}x = \lim_{N \to \infty} \sum\limits_{n=0}^{N-1} f(a + n \Delta x) \ \Delta x \qquad \qquad \text{where} \quad \Delta x \triangleq \frac{b-a}{N}$$

Now if $N$ is just left as large and finite (and even, just to make our lives easier), then the two integrals above (with finite limits) have approximations that look like:

$$\begin{align} X(f) &\approx \int\limits_{-\frac{T}2}^{+\frac{T}2} x(t) \ e^{-j 2 \pi f t} \ \mathrm{d}t \\ &\approx \sum\limits_{n=0}^{N-1} x(-\tfrac{T}2 + n \Delta t) \ e^{-j 2 \pi f (-\frac{T}2 + n \Delta t)} \ \Delta t \\ &= \sum\limits_{n=-\frac{N}2}^{\frac{N}2-1} x(n \Delta t) \ e^{-j 2 \pi f (n \Delta t)} \ \Delta t \\ \end{align}$$

where $\qquad \Delta t = \frac{T}{N}$.

$$\begin{align} x(t) &\approx \int\limits_{-\frac{f_\mathrm{s}}{2}}^{+\frac{f_\mathrm{s}}{2}} X(f) \ e^{+j 2 \pi f t} \ \mathrm{d}f \\ &\approx \sum\limits_{k=0}^{N-1} X(-\tfrac{f_\mathrm{s}}2 + k \Delta f) \ e^{+j 2 \pi (-\tfrac{f_\mathrm{s}}2 + k \Delta f) t} \ \Delta f \\ &= \sum\limits_{k=-\frac{N}2}^{\frac{N}2-1} X(k \Delta f) \ e^{+j 2 \pi (k \Delta f) t} \ \Delta f \\ \end{align}$$

where $\qquad \Delta f = \frac{f_\mathrm{s}}{N}$.

Here we need to recognize $\Delta t$ as the sampling period, same as $\frac{1}{f_\mathrm{s}}$, which means that

$$\begin{align} \Delta f &= \frac{f_\mathrm{s}}{N} \\ &= \frac{1}{N \ \Delta t} \\ \end{align}$$

or $\qquad N \ \Delta f \ \Delta t = 1 $.


So, to relate this to the DFT, let's define the discrete-time samples as:

$$ x[n] \triangleq x(n \Delta t) $$

When there are square brackets, the argument must be an integer. So "$x[n]$" is exactly like "$x_n$".

The DFT and inverse are

$$ X[k] = \sum\limits_{n=0}^{N-1} x[n] \ e^{-j2\pi nk/N} $$

$$ x[n] = \tfrac{1}N \sum\limits_{k=0}^{N-1} X[k] \ e^{+j2\pi nk/N} $$

Now there are DFT periodicity deniers hanging around here that deny this, but it is simply true that:

$$\begin{align} x[n+N] &= x[n] \qquad &\forall n \in \mathbb{Z} \\ X[k+N] &= X[k] \qquad &\forall k \in \mathbb{Z} \\ \end{align}$$

This means that the DFT and inverse can have the limits in the sum shifted by any integer amount.

$$ X[k] = \sum\limits_{n=n_0}^{n_0+N-1} x[n] \ e^{-j2\pi nk/N} \qquad \forall n_0 \in \mathbb{Z} $$

$$ x[n] = \tfrac{1}N \sum\limits_{k=k_0}^{k_0+N-1} X[k] \ e^{+j2\pi nk/N} \qquad \forall k_0 \in \mathbb{Z} $$

We can pick $n_0=k_0=-\frac{N}{2}$:

$$ X[k] = \sum\limits_{n=-\frac{N}{2}}^{\frac{N}{2}-1} x[n] \ e^{-j2\pi nk/N} $$

$$ x[n] = \tfrac{1}N \sum\limits_{k=-\frac{N}{2}}^{\frac{N}{2}-1} X[k] \ e^{+j2\pi nk/N} $$

So putting it together, we recognize that $\Delta t\Delta f = \frac{1}N $ and we evaluate $X(f)$ at discrete frequencies, $k\Delta f$,

$$\begin{align} X(f) \Big|_{f=k\Delta f} &= \sum\limits_{n=-\frac{N}2}^{\frac{N}2-1} x(n \Delta t) \ e^{-j 2 \pi f (n \Delta t)} \ \Delta t \Big|_{f=k\Delta f} \\ &= \sum\limits_{n=-\frac{N}2}^{\frac{N}2-1} x(n \Delta t) \ e^{-j 2 \pi (k\Delta f) (n \Delta t)} \ \Delta t \\ &= \sum\limits_{n=-\frac{N}2}^{\frac{N}2-1} x[n] \ e^{-j 2 \pi nk/N} \ \Delta t \\ &= X[k] \cdot \Delta t \\ &= X[k] \cdot \frac{1}{f_\mathrm{s}}\\ \end{align}$$

So your FFT output value is $X[k]=X(k\Delta f) \cdot f_\mathrm{s}$ whereas the input value was defined above to be $x[n]=x(n\Delta t)$. Now magnitude squaring we have

$$\begin{align} \Big|X[k]\Big|^2 &= \Big|X(k\Delta f)\Big|^2 \cdot f_\mathrm{s}^2 \\ \\ &= \frac{1}{N \Delta t} \cdot \Big|X(k\Delta f)\Big|^2 \cdot N \ f_\mathrm{s} \\ \\ &= \frac{1}{T} \Big|X(k\Delta f)\Big|^2 \cdot N \ f_\mathrm{s} \\ \end{align}$$

If $x(t)$ (and also $x[n]$) are in volts, then as above $\frac{1}T |X(f)|^2$ would be "volts² per Hz" in the vicinity of frequency $f$. Then at frequency $k \Delta f = \frac{k}{N} f_\mathrm{s}$, the magnitude-square of the corresponding point in the FFT, scaled down by $\frac{1}N$, is

$$ \tfrac{1}N \Big|X[k]\Big|^2 = \tfrac{1}{T} \Big|X(k\Delta f)\Big|^2 \cdot f_\mathrm{s} $$

which would be "volts² per Hz times the sample rate in Hz" or just volts² at frequency $\frac{k}{N} f_\mathrm{s}$.

| improve this answer | |
$\endgroup$
  • 1
    $\begingroup$ Nobody that I am aware of has denied the periodicity of either the extrapolation of the DFT nor of its inverse. What is incorrect is to identify the extrapolation as the signal unless it is specified that it is a periodic signal with the DFT frame covering a period (not necessarily the fundamental). Technical precision is important in mathematical language. $\endgroup$ – Cedron Dawg Sep 7 at 4:29
  • 2
    $\begingroup$ @Q.P. , thank you. while i have been around since 2013, and i have several times spent cred points on a bounty, this is the first time that i have ever received a bounty. $\endgroup$ – robert bristow-johnson Sep 7 at 21:16
  • 1
    $\begingroup$ @robertbristow-johnson. Congrats! :-) $\endgroup$ – Peter K. Sep 7 at 21:21
  • 1
    $\begingroup$ @robertbristow-johnson it seemed only fair! While your reputation clearly doesn't "need" them, you put such an enormous amount of effort into a very well explained answer, I felt it warranted a symbolic gesture. $\endgroup$ – Q.P. Sep 7 at 23:36
  • 1
    $\begingroup$ It's okay to say $$ x[n] = x[n+N] $$ when you are talking about the extension (extrapolation on the set of integers) of the inverse DFT (and again I think DFS should stand for Discrete Fourier Spectrum, not Discrete Fourier Series, but it's fine with me calling the extension/extrapolation that and I don't know what else it would mean). But it is not okay to say $$ x[n] = x[n+N] $$ about the signal. Likewise, it is not correct to say the DFT "assumes it". The DFT is blind to it. $\endgroup$ – Cedron Dawg Sep 8 at 4:22
4
$\begingroup$

You do this by doing two things:

  1. First establish what the relationship of scale exists between your A/D converter and the numerical values that go into your DFT. Find out how many volts corresponds to a value of $1.0$ in the samples that go into the DFT.

  2. The next thing you gotta do is express your integral above without that $\frac{1}{T}$ factor (which is an error) as a Riemann summation having equal-width little rectangles of width equal to the sampling period. That will become identical to the summation of the DFT.

That should spell it out.

| improve this answer | |
$\endgroup$
  • $\begingroup$ I'm really pleased you said point 2. because just "playing around" and removing that $1/T$ makes my final function match what I see when I simulate the DFT. So thanks! Would you say what I am doing, bearing in mind what I want to achieve, makes sense? $\endgroup$ – Q.P. Aug 31 at 13:35
  • 1
    $\begingroup$ yes. just remember that the "$dt$" in the integral has dimension of time (likely units of seconds). so you'll get $\mathrm{V/Hz}$ coming outa the integral. that means you do need to multiply by some quantity of time to make the result $\mathrm{V}^2\mathrm{/Hz}$. i think you have to multiply by the sample rate. $\endgroup$ – robert bristow-johnson Aug 31 at 17:12
  • 1
    $\begingroup$ maybe i'm wrong about leaving off the $\frac1T$. i gotta think about this a little more. $\endgroup$ – robert bristow-johnson Sep 1 at 1:08
  • 1
    $\begingroup$ i'm gonna add another answer and do this formaly. $\endgroup$ – robert bristow-johnson Sep 1 at 17:21
  • 1
    $\begingroup$ @Q.P. Check out dsp.stackexchange.com/questions/69430/… Clearly the $dt$ corresponds to a $1/N$ normalization. $\endgroup$ – Cedron Dawg Sep 2 at 11:53
1
$\begingroup$

According to this, the units for PSD from a DFT should be volts^2/bin:
https://www.mathworks.com/matlabcentral/answers/47633-what-is-the-relation-between-dft-and-psd-of-a-signal

$$ \mathrm{PSD} = [ X[k] \cdot \operatorname{conj}(X[k]) ] / N $$

$$ \mathrm{units} = \mathrm{volts} \cdot \mathrm{volts} /\mathrm{bins} = \mathrm{volt}^2/\mathrm{bin} $$

Which makes sense since when it is applied to a signal measured in seconds, there is a bin width value of so many Hz per bin (given by $f_s/N$). Thus, you have a conversion factor between the DFT PSD units and the volts^2/Hz you were expecting.

| improve this answer | |
$\endgroup$

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.