137 lines
5.0 KiB
Markdown
137 lines
5.0 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 (<a href="#citeproc_bib_item_1">Schmid 2012</a>).
|
||
|
||
|
||
## 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](#figure--fig:psd-manual-example).
|
||
|
||
<a id="figure--fig:psd-manual-example"></a>
|
||
|
||
{{< figure src="/ox-hugo/psd_manual_example.png" caption="<span class=\"figure-number\">Figure 1: </span>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](#figure--fig:psd-comp-pwelch-manual-example).
|
||
|
||
<a id="figure--fig:psd-comp-pwelch-manual-example"></a>
|
||
|
||
{{< figure src="/ox-hugo/psd_comp_pwelch_manual_example.png" caption="<span class=\"figure-number\">Figure 2: </span>Amplitude Spectral Density with manual computation" >}}
|
||
|
||
|
||
## Bibliography {#bibliography}
|
||
|
||
<style>.csl-entry{text-indent: -1.5em; margin-left: 1.5em;}</style><div class="csl-bib-body">
|
||
<div class="csl-entry"><a id="citeproc_bib_item_1"></a>Schmid, Hanspeter. 2012. “How to Use the Fft and Matlab’s Pwelch Function for Signal and Noise Simulations and Measurements.” <i>Institute of Microelectronics</i>.</div>
|
||
</div>
|