1
$\begingroup$

I'm trying to acquire a better understanding of the heuristics behind power spectral density (PSD). In particular, for a very simple function, is it possible to determine the peak PSD values from graphical inspection?

Here's a simple Python example, where my signal is a linear combination of two sine functions having frequencies 150 and 375 with corresponding weights 5 and 20. The signal length is .5 seconds and the sampling rate is 1000 Hz.

import numpy as np
from matplotlib import pyplot as plt
from scipy.signal import periodogram

t = np.arange(0, .5, 0.001);
dt = t[1] - t[0]
fs=1/dt

x = 5*np.sin(2*np.pi*150*t)+20*np.sin(2*np.pi*375*t)

f,Pxx=periodogram(x,nfft=len(x))
fig, ax = plt.subplots()
ax.plot(f,Pxx)
plt.show()

The resulting graph looks as follows:

enter image description here

The frequency vector has length 251 and my sampling frequency is fs=1000, so fs/2=500. Rescaling my horizontal axis by fs, I can see the peaks correspond to the frequencies 150 and 375.

But what about the peak PSD values, i.e. heights? Is there a way to determine or estimate their values based upon the weights of 5 and 20 used above. For example, the maximum PSD is approximately 40,875.

$\endgroup$
1
  • $\begingroup$ "heuristics" is not the right word for a mathematical relationship :) Welcome here! $\endgroup$ Mar 29, 2022 at 18:48

1 Answer 1

2
$\begingroup$

TL; DR: Yes, there is. Apply the DFT to your signal, it doesn't get any easier than this. If your signal is composed of discrete harmonic oscillations, the peaks are at their frequency and their power is proportional to the square of their amplitude.


the heights, given you'd be estimating the PSD exactly, would be 5² and 20² scaled by the vector length (often divided by two).

That's just how the DFT, which your periodogram function uses, is defined.

Other DFT implementations use a scaling of the inverse square root of the vector length or just 1.

To understand: the DFT of a vector $\mathbf v$ of length $N$ is just another vector of the same length, its formula being

$$\text{DFT}(\mathbf v) = \left(\sum_{k=0}^{N-1} \mathbf v_k e^{-j2\pi\frac nNk}\right)_{n=0,\ldots,N-1}.$$

So, if you put in a single sine of amplitude $a$ and frequency $f$ (because we're looking at a discrete transform, let's constrain ourselves to $f$ that are multiple of $\frac1N$), $$ \mathbf v_k= a\sin(2\pi fk)=a\frac 1{2i}\left(e^{j2\pi fk} - e^{-j2\pi fk}\right),$$

you get

$$\text{DFT}(\mathbf v) = \left(\frac a{2i} \sum_{k=0}^{N-1} \left( e^{j2\pi fk} - e^{-j2\pi fk}\right) e^{-j2\pi\frac nNk}\right)_{n=0,\ldots,N-1},$$

and if you squint really hard, you can multiply out all the complex exponentials and get new complex exponentials; they all do one or multiple "rounds" around the unit circle for $k$ going from $0$ to $N-1$. Which means their sum just cancels. Aside from the one, where the $\frac nN$ happens to be the same as $f$ – and then, you get a sum full of $N$ individual $e^0=1$. Hence, the DFT at thatfrequency has the vale $N\cdot \frac a{2i}$; and since your DFT-based PSD estimate takes the magnitude square, that's

$$\left\lvert N\cdot \frac a{2i}\right\rvert^2 =\frac{N^2}{2^2}\left\lvert \frac a{i}\right\rvert^2=\frac{N^2}{2^2}\left\lvert a\right\rvert^2.$$

Typically, since your display is one-sided, you only see the positive frequency part, and the display scales it with two, to signify that the energy at the negative frequencies does alos come into effect.

$\endgroup$
2
  • $\begingroup$ Thanks for taking the time to provide such an in-depth response. A couple issues I'm still having a hard time grasping. Suppose x(t)=20*sin(375*t) and t is as above (.5 seconds with fs=1000) so that N=500 and a=20. Then N^2*|a|^2/4= 25,000,000. On the other hand, using f,Pxx=periodogram(x,nfft=len(x)) whose maximum is about 40784 $\endgroup$
    – fishbacp
    Mar 30, 2022 at 15:49
  • $\begingroup$ To get the answers to agree, the N^2/4*|a)^2 must be divided by the number of frequency bins. $\endgroup$
    – fishbacp
    Mar 30, 2022 at 19:26

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.