130 lines
4.4 KiB
Markdown
130 lines
4.4 KiB
Markdown
+++
|
|
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](https://research.tdehaeze.xyz/spectral-analysis/).
|
|
|
|
A good article about how to use the `pwelch` function with Matlab <schmid12_how_to_use_fft_matlab>.
|
|
|
|
|
|
## Parseval's Theorem - Linking the Frequency and Time domain {#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) {#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 {#matlab-code-for-computing-the-psd-and-cps}
|
|
|
|
Let's compute the PSD of a signal by "hand".
|
|
The signal is defined below.
|
|
|
|
```matlab
|
|
%% 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.
|
|
|
|
```matlab
|
|
%% 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](#org41c99c6).
|
|
|
|
<a id="org41c99c6"></a>
|
|
|
|
{{< figure src="/ox-hugo/psd_manual_example.png" caption="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.
|
|
|
|
```matlab
|
|
%% 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](#orge7a31a8).
|
|
|
|
<a id="orge7a31a8"></a>
|
|
|
|
{{< figure src="/ox-hugo/psd_comp_pwelch_manual_example.png" caption="Figure 2: Amplitude Spectral Density with manual computation" >}}
|