Other than the time plot, the power spectrum and the spectrogram (a plot that shows the signal in time and frequency domain) are probably the most important tools for digital signal analysis. Because of its importance, people invented so many algorithms to estimate and present it in various ways. In Matlab, you can use periodogram, pwelch, pspectrum, pburg, pmusic, and many more of them. It is so overwhelming for me at the beginning (and even now), but from time to time, I gradually get used to it and make sense of it. In this blog, I am going to talk about a few power spectrum-related short topics that confused me in the early days.
Why power spectrum
So why bother power spectrum, why can’t we apply the Fourier Transform directly to the signal and get the frequency spectrum? This may be the first question many people will ask.
If the signal is periodic or more specifically, deterministic, we can use the frequency spectrum to present the signal. However, most signals we are dealing with are stochastic (random, non-deterministic) signals, or with a stochastic component. In addition, these signals are usually infinite in length, and we can’t apply an infinite length of Fourier Transform. If we truncate the stochastic signal into segments and apply the Fourier Transform to each segment, we can find the frequency spectra are random in both magnitude and phase at a given frequency.
The frequency spectrum analysis is only applicable to deterministic signals. As for stochastic signals, they can be analyzed statistically, or in our case, be averaged. Because of the random phase mentioned before, the mean of frequency spectra is probably zero. Although the frequency spectra are complex-valued, the power spectra are real-valued and can be averaged. The periodogram method is essentially the averaged power spectra alone multiply signal segments.
What is the power spectral density
When we talk about the power spectrum, what we’re actually talking about is the power spectral density (PSD). Mathematically, we have to normalize the power in the frequency domain by the sampling rate. In another word, the value of PSD is the power per frequency unit.
If we take a look at the plot from Matlab’s periodogram or the psd from Python’s matplotlib, the unit of the y-axis is dB/Hz when the sampling rate is given. But why do we want to normalize the power by the sampling rate? Let’s have a Python demo to find it out.
We first generate two sine waves with
Different frequencies (for visualization purposes).
The sampling rates are 100 Hz and 200 Hz respectively.
The same amplitude and same length in time.
So they have the same energy, but the number of samples is different.
1
2
3
4
5
6
7
8
9
10
11
12
import numpy as np
import matplotlib.pyplot as plt
f_sig1 = 15
fs1 = 100
t1 = np.arange(0, 512/fs1, 1/fs1)
sig1 = np.sin(2 * np.pi * f_sig1 * t1)
f_sig2 = 25
fs2 = 200
t2 = np.arange(0, 1024/fs2, 1/fs2)
sig2 = np.sin(2 * np.pi * f_sig2 * t2)
If you are trying to calculate the energy of the two signals directly, you may find the energy of the second signal is twice the first one. That is because you forgot the time step during the integration. This will happen in the frequency domain too.
Next, let’s find the power spectrum of the two sine waves. The numbers of samples for the two sine waves are deliberately set to 512 and 1024, so we can use the FFT size of 512 and 1024 directly. Note that, according to the Parseval’s theorem:
\[\sum_{n=0}^{N-1} | x[n] |^2 = \frac{1}{N} \sum_{k=0}^{N-1} | X[k] |^2\]To preserve the energy, the power in the frequency domain needs to be divided by N, and in our case, the FFT size.
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
nfft1 = 512
ft1 = np.fft.fft(sig1, nfft1)
p1 = (ft1 * ft1.conj()).real/nfft1
f1 = np.fft.fftfreq(nfft1, 1/fs1)
nfft2 = 1024
ft2 = np.fft.fft(sig2, nfft2)
p2 = (ft2 * ft2.conj()).real /nfft2
f2 = np.fft.fftfreq(nfft2, 1/fs2)
plt.plot(f1, p1)
plt.plot(f2, p2)
plt.xlim((10,30))
plt.xlabel('frequency (Hz)')
plt.ylabel('power?')
plt.show()
We can find that the second sine wave shows twice the power of the first one, which is obviously wrong. To fix it, we can normalize the power in the frequency domain by the sampling rate, then the two power peaks (in the PSD plot) will have the same height.
When preparing for this post, I found an interesting example in this Matlab document. If you take a close look at the noise level of the same signal in the power spectrum and the persistence spectrum (we will talk about it in a later post), you can find the former is 6 dB lower than the latter. This is also because they are normalized differently. Can you figure out why the difference is 6 dB? 
Matlab didn’t tell you everything
to be continued
Comments powered by Disqus.