digital-brain/content/zettels/power_spectral_density.md

+++
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" >}}