Again a signal processing post. I'm writing this because I couldn't find a complete and concise explanation of practically computing correct power spectral densities (PSD) from DFT online.

Generating a noise signal

We will consider a complex random process $s(t)$ of variance $\sigma^2$. $$s(t) = \frac{\sigma}{\sqrt{2}}(n_i + jn_q)$$ The two random variables $n_i$ and $n_q$ being sampled from the normal centered distribution $\mathcal{N}(0, \sigma^2)$.

Because we are not just writing equations but solving real world problems, let's write code that can actually generate this signal.

import numpy as np
σ, fs, t = 5, 100, 10
n = t * fs
s = (sigma / np.sqrt(2)) * (np.random.randn(n) + np.random.randn(n))

With the two parameters fs and t being the sampling rate in Hz and duration in seconds.

Discrete Fourier Transform (DFT)

The DFT is defined as the following operation:

$$S_k = \sum_{n=0}^{N-1} s_n \cdot e^{-\frac {i 2\pi}{N}kn}$$

You can use np.fft module or code it your own (it will be slow, don't do that unless you want to have fun). As a different example, I here compute the dot product with complex exponentials.

t_vec = np.arange(0, t, 1./fs)
f_vec = np.linspace(-fs/2, fs/2, n)
S = np.zeros(n, dtype=np.complex128)
for i, f in enumerate(f_vec):
    S[i] = np.dot(s, np.exp(-1j*2*np.pi*f*t_vec))
# or S = np.fft.fft(s)

From our DFT result S, let's move on to the Energy Spectrum Density (ESD) computation.

Energy Spectrum Density

Again, from definition: $$\bar{S}_{xx}(f) \triangleq \left |S(f) \right |^2 $$ Let's take the square modulus of the DFT we just computed.

esd = np.abs(S)**2 

Now, this is supposed to be a density, which whole purpose in life is to be integrated to get some sort of energy. Here is the first little detail that needs attention. From Parseval's theorem we know that:

$$\int_{-\infty}^\infty | x(t) |^2 \ \mathrm{d}t = \frac{1}{2\pi} \int_{-\infty}^\infty | X(\omega) |^2 \ \mathrm{d}\omega$$

But in the discrete world, this looks like: $$\sum_{n=0}^{N-1} | x[n] |^2 = \frac{1}{N} \sum_{k=0}^{N-1} | X[k] |^2$$ Notice the $\frac{1}{N}$ factor ? Don't forget it if you want to get back the total energy of your signal from the spectral estimation. Just to check, I will actually compute it.

esd.sum()/n # don't forget this 1/n factor
> 25375.0

Compare it against energy computed from temporal samples:

np.sum(np.abs(s)**2)
> 25485.9

We are doing fine.

Power Spectrum Density

The real meat here. Always, starting from the definition: $$S_{xx}(f) \triangleq \lim_{T\to \infty} \frac 1 {T} |\hat{x}_{T}(f)|^2$$ This seems to be the Energy Spectrum Density divided by the duration, taken at the limit of time. Practically, let's divide by the number of samples in our signal.

psd = esd / n

Simple enough. But can we retrieve the signal power from it ? For a 0-mean noise, the average power of the signal should be it's variance. If our PSD computation is wrong, we won't get it back.

First, computing the integrated power, not forgetting the $\frac{1}{N}$ factor from our old friend Parseval.

psd.sum()*parseval
> 25.4

Ok, 25. This should be the variance of our signal, which is $\sigma^2$.

σ**2
> 25

Well, good ! We also could have computed empirically the variance from the generated samples.

s.var()
> 25.5

This seems so simple, if you know how and why you do it. Obviously, things often get a little bit more complicated, like when windowing signals.

Windowing signals ?

Windowing signals, yes

Why window signals ? It worked before ! Well, we were already windowing the signal, and in a bad way.

What we were essentially doing when generating $s(t)$ was applying a rectangular window to it, by default. In time domain, multiplying a signal by a rectangle window is the same as convoluting the fourier transform of the rectangle window with the signal spectrum. Because the fourier transform of a rectangle window is an ugly sinc function, a lot of secondary lobes are observed.

$$\operatorname{sinc}x = \frac{\sin x}{x}$$

By applying smooth windows (like hamming, blackman, etc...) with prettier frequency responses (lower lobes), we get a more accurate depiction of the actual signal we are receiving.

Does it change the computation of the PSD ? well yes and no. Windowing reduces amplitudes at the start and end of the signal, effectively reducing the energy of the signal. We need to compensate this loss when computing the energy, or the power.

Let's create an hamming window:

$$w[n] = 0.54\ \left[1 - \cos \left ( \frac{2 \pi n}{N} \right) \right]$$

w = np.hamming(n)

And apply it to our signal.

s *= w

Going back to the ESD computation, we simply (but nonetheless need) to normalize back with the energy ratio loss of the window if we were to integrate it. This, my friend, will annoy if you forget it.

esd = np.abs(S)**2 
esd *= n/np.dot(w, w)

The ratio n/np.dot(w, w) effectively correct the energy lost, n being the rectangular window energy.

Wrapping up

This seemed so simple, but still it had to be done correctly, accounting for the details and understanding them.

I hope this writing was useful for you and, who knows, even for future me ?