5.0 KiB
+++ title = "Power Spectral Density" author = ["Dehaeze Thomas"] draft = false +++
- Tags
- [Signal to Noise Ratio]({{< relref "signal_to_noise_ratio.md" >}})
Tutorial about Power Spectral Density is accessible here.
A good article about how to use the pwelch function with Matlab (Schmid 2012).
Parseval's Theorem - Linking the Frequency and Time domain
For non-periodic finite duration signals, the energy in the time domain is described by:
\begin{equation} \text{Energy} = \int_{-\infty}^\infty x(t)^2 dt \end{equation}
Parseval's Theorem states that energy in the time domain equals energy in the frequency domain:
\begin{equation} \text{Energy} = \int_{-\infty}^{\infty} x(t)^2 dt = \int_{-\infty}^{\infty} |X(f)|^2 df \end{equation}
where \(X(f)\) is the Fourier transform of the time signal \(x(t)\):
\begin{equation} X(f) = \int_{-\infty}^{\infty} x(t) e^{-2\pi j f t} dt \end{equation}
Power Spectral Density function (PSD)
The power distribution over frequency of a time signal \(x(t)\) is described by the PSD denoted the \(S_x(f)\). A PSD is a power density function with units \([\text{SI}^2/Hz]\), meaning that the area underneath the PSD curve equals the power (units \([\text{SI}^2]\)) of the signal (SI is the unit of the signal, e.g. \(m/s\)).
Using the definition of signal power \(\bar{x^2}\) and Parseval's theorem, we can link power in the time domain with power in the frequency domain:
\begin{equation} \text{power} = \lim_{T \to \infty} \frac{1}{2T} \int_{-T}^{T} x_T(t)^2 dt = \lim_{T \to \infty} \frac{1}{2T} \int_{-\infty}^{\infty} |X_T(f)|^2 df = \int_{-\infty}^{\infty} \left( \lim_{T \to \infty} \frac{|X_T(f)|^2}{2T} \right) df \end{equation}
where \(X_T(f)\) denotes the Fourier transform of \(x_T(t)\), which equals \(x(t)\) on the interval \(-T \le t \le T\) and is zero outside this interval.
This term is referred to as the two-sided spectral density:
\begin{equation} S_{x,two} (f) = \lim_{T \to \infty} \frac{|X_T(f)|^2}{2T}, \quad -\infty \le f \le \infty \end{equation}
In practice, the one sided PSD is used, which is only defined on the positive frequency axis but also contains all the power. It is defined as:
\begin{equation} S_{x}(f) = \lim_{T \to \infty} \frac{|X_T(f)|^2}{T}, \quad 0 \le f \le \infty \end{equation}
For discrete time signals, the one-sided PSD estimate is defined as:
\begin{equation} \hat{S}(f_k) = \frac{|X_L(f_k)|^2}{L T_s} \end{equation}
where \(L\) equals the number of time samples and \(T_s\) the sample time, \(X_L(f_k)\) is the N-point discrete Fourier Transform of the discrete time signal \(x_L[n]\) containing \(L\) samples:
\begin{equation} X_L(f_k) = \sum_{n = 0}^{N-1} x_L[n] e^{-j 2 \pi k n/N} \end{equation}
Matlab Code for computing the PSD and CPS
Let's compute the PSD of a signal by "hand". The signal is defined below.
%% Signal generation
T_s = 1e-3; % Sampling Time [s]
t = T_s:T_s:100; % Time vector [s]
L = length(t);
x = lsim(1/(1 + s/2/pi/5), randn(1, L), t);
The computation is performed using the fft function.
%% Parameters
T_r = L*T_s; % signal time range
d_f = 1/T_r; % width of frequency grid
F_s = 1/T_s; % sample frequency
F_n = F_s/2; % Nyquist frequency
F = [0:d_f:F_n]; % one sided frequency grid
% Discrete Time Fourier Transform Wxx
Wxx = fft(x - mean(x))/L;
% Two-sided Power Spectrum Pxx [SI^2]
Pxx = Wxx.*conj(Wxx);
% Two-sided Power Spectral Density Sxx_t [SI^2/Hz]
Sxx_t = Pxx/d_f;
% One-sided Power Spectral Density Sxx_o [SI^2/Hz] defined on F
Sxx_o = 2*Sxx_t(1:L/2+1);
The result is shown in Figure 1.
{{< figure src="/ox-hugo/psd_manual_example.png" caption="<span class="figure-number">Figure 1: Amplitude Spectral Density with manual computation" >}}
This can also be done using the pwelch function which integrated a "window" that permits to do some averaging.
%% Computation using pwelch function
[pxx, f] = pwelch(x, hanning(ceil(5/T_s)), [], [], 1/T_s);
The comparison of the two method is shown in Figure 2.
{{< figure src="/ox-hugo/psd_comp_pwelch_manual_example.png" caption="<span class="figure-number">Figure 2: Amplitude Spectral Density with manual computation" >}}