Making the units of a analytically calculated PSD consistent with the units of an FFT

Question

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:

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.

Answer 1

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

Answer 2

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.

Answer 3

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.