#+TITLE: Compute Spectral Densities of signals with Matlab
:DRAWER:
#+LANGUAGE: en
#+EMAIL: dehaeze.thomas@gmail.com
#+AUTHOR: Dehaeze Thomas

#+HTML_HEAD: <link rel="stylesheet" type="text/css" href="./css/htmlize.css"/>
#+HTML_HEAD: <link rel="stylesheet" type="text/css" href="./css/readtheorg.css"/>
#+HTML_HEAD: <link rel="stylesheet" type="text/css" href="./css/zenburn.css"/>
#+HTML_HEAD: <script type="text/javascript" src="./js/jquery.min.js"></script>
#+HTML_HEAD: <script type="text/javascript" src="./js/bootstrap.min.js"></script>
#+HTML_HEAD: <script type="text/javascript" src="./js/jquery.stickytableheaders.min.js"></script>
#+HTML_HEAD: <script type="text/javascript" src="./js/readtheorg.js"></script>

#+PROPERTY: header-args:latex  :headers '("\\usepackage{tikz}" "\\usepackage{import}" "\\import{$HOME/MEGA/These/LaTeX/}{config.tex}")
#+PROPERTY: header-args:latex+ :imagemagick t :fit yes
#+PROPERTY: header-args:latex+ :iminoptions -scale 100% -density 150
#+PROPERTY: header-args:latex+ :imoutoptions -quality 100
#+PROPERTY: header-args:latex+ :results raw replace :buffer no
#+PROPERTY: header-args:latex+ :eval no-export
#+PROPERTY: header-args:latex+ :exports both
#+PROPERTY: header-args:latex+ :mkdirp yes
#+PROPERTY: header-args:latex+ :output-dir figs
#+PROPERTY: header-args:latex+ :post pdf2svg(file=*this*, ext="png")

#+PROPERTY: header-args:matlab  :session *MATLAB*
#+PROPERTY: header-args:matlab+ :tangle filters.m
#+PROPERTY: header-args:matlab+ :comments org
#+PROPERTY: header-args:matlab+ :exports both
#+PROPERTY: header-args:matlab+ :results none
#+PROPERTY: header-args:matlab+ :eval no-export
#+PROPERTY: header-args:matlab+ :noweb yes
#+PROPERTY: header-args:matlab+ :mkdirp yes
#+PROPERTY: header-args:matlab+ :output-dir figs
:END:

This document presents the mathematics as well as the matlab scripts to do the spectral analysis of a measured signal.

Typically this signal is coming from an inertial sensor, a force sensor or any other sensor.

We here take the example of a signal coming from a Geophone measurement the vertical velocity of the floor at the ESRF.

* Matlab Init                                               :noexport:ignore:
#+begin_src matlab :tangle no :exports none :results silent :noweb yes :var current_dir=(file-name-directory buffer-file-name)
  <<matlab-dir>>
#+end_src

#+begin_src matlab :exports none :results silent :noweb yes
  <<matlab-init>>
#+end_src

* Sensitivity of the instrumentation
The measured signal $x$ by the ADC is in Volts.
The corresponding real velocity $v$ in m/s.

To obtain the real quantity as measured by the sensor, one have to know the sensitivity of the sensors and electronics used.

#+begin_src latex :file velocity_to_voltage.pdf :post pdf2svg(file=*this*, ext="png") :exports results
  \begin{tikzpicture}
    \node[block] (geophone) at (0, 0) {$G_g(s)$};
    \node[above] at (geophone.north) {Geophone};
    \node[block, right=1 of geophone] (ampli) {$G_a(s)$};
    \node[above] at (ampli.north) {Amplifier};
    \node[ADC, right=1 of ampli] (adc) {ADC};

    \draw[double, <-] (geophone.west) -- node[midway, above]{$v$ [m/s]} ++(-1.4, 0);
    \draw[->] (geophone.east) -- node[midway, above]{[V]} (ampli.west);
    \draw[->] (ampli.east) -- node[midway, above]{[V]} (adc.west);
    \draw[->] (adc.east) -- node[sloped]{$/$}node[midway, above]{$x$ [V]} ++(1.4, 0);
  \end{tikzpicture}
#+end_src

#+NAME: fig:velocity_to_voltage
#+CAPTION: Schematic of the instrumentation used for the measurement
#+RESULTS: fig:velocity_to_voltage
[[file:figs/velocity_to_voltage.png]]

* Convert the time domain from volts to velocity
Let's say, we know that the sensitivity of the geophone used is
\[ G_g(s) = G_0 \frac{\frac{s}{2\pi f_0}}{1 + \frac{s}{2\pi f_0}} \quad \left[\frac{V}{m/s}\right] \]

#+begin_src matlab
  G0 = 88; % Sensitivity [V/(m/s)]
  f0 = 2; % Cut-off frequency [Hz]

  Gg = G0*(s/2/pi/f0)/(1+s/2/pi/f0);
#+end_src

And the gain of the amplifier is 1000: $G_m(s) = 1000$.

#+begin_src matlab
  Gm = 1000;
#+end_src

If ${G_m(s)}^{-1} {G_g(s)}^{-1}$ is proper, we can simulate this dynamical system to go from the voltage to the velocity units (figure [[fig:voltage_to_velocity]]).

#+begin_src matlab
  data = load('mat/data_028.mat', 'data'); data = data.data;

  t = data(:, 3); % [s]
  x = data(:, 1)-mean(data(:, 1)); % The offset if removed (coming from the voltage amplifier) [v]

  dt = t(2)-t(1); Fs = 1/dt;
#+end_src

#+begin_src latex :file voltage_to_velocity.pdf :post pdf2svg(file=*this*, ext="png") :exports results
  \begin{tikzpicture}
    \node[block] (ampli) at (0, 0) {${G_a(s)}^{-1}$};
    \node[above] at (ampli.north) {Amplifier};
    \node[block, right=1 of ampli] (geophone) {${G_g(s)}^{-1}$};
    \node[above] at (geophone.north) {Geophone};

    \draw[<-] (ampli.west) -- node[midway, above]{$x$ [V]}node[sloped]{$/$} ++(-1.4, 0);
    \draw[->] (ampli.east) -- node[midway, above]{[V]}node[sloped]{$/$} (geophone.west);
    \draw[->] (geophone.east) -- node[midway, above]{$v$ [m/s]}node[sloped]{$/$} ++(1.4, 0);
  \end{tikzpicture}
#+end_src

#+NAME: fig:voltage_to_velocity
#+CAPTION: Schematic of the instrumentation used for the measurement
#+RESULTS: fig:voltage_to_velocity
[[file:figs/voltage_to_velocity.png]]

We simulate this system with matlab:
#+begin_src matlab
  v = lsim(inv(Gg*Gm), v, t);
#+end_src

And we plot the obtained velocity
#+begin_src matlab
  figure;
  plot(t, v);
  xlabel("Time [s]"); ylabel("Velocity [m/s]");
#+end_src

#+NAME: fig:velocity_time
#+HEADER: :tangle no :exports results :results value raw replace :noweb yes
#+begin_src matlab :var filepath="figs/velocity_time.pdf" :var figsize="wide-tall" :post pdf2svg(file=*this*, ext="png")
  <<plt-matlab>>
#+end_src

#+NAME: fig:velocity_time
#+CAPTION: Measured Velocity
#+RESULTS: fig:velocity_time
[[file:figs/velocity_time.png]]

* Power Spectral Density and Amplitude Spectral Density
We now have the velocity in the time domain:
\[ v(t)\ [m/s] \]

To compute the Power Spectral Density (PSD):
\[ S_v(f)\ \left[\frac{(m/s)^2}{Hz}\right] \]

To compute that with matlab, we use the =pwelch= function.

We first have to defined a window:
#+begin_src matlab
  win = hanning(ceil(10*Fs)); % 10s window
#+end_src

#+begin_src matlab
  [Sv, f] = pwelch(v, win, [], [], Fs);
#+end_src

#+begin_src matlab
  figure;
  loglog(f, Sv);
  xlabel('Frequency [Hz]');
  ylabel('Power Spectral Density $\left[\frac{(m/s)^2}{Hz}\right]$')
#+end_src

#+NAME: fig:psd_velocity
#+HEADER: :tangle no :exports results :results value raw replace :noweb yes
#+begin_src matlab :var filepath="figs/psd_velocity.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
  <<plt-matlab>>
#+end_src

#+NAME: fig:psd_velocity
#+CAPTION: Power Spectral Density of the measured velocity
#+RESULTS: fig:psd_velocity

The Amplitude Spectral Density (ASD) is the square root of the Power Spectral Density:
\begin{equation}
  \Gamma_{vv}(f) = \sqrt{S_{vv}(f)} \quad \left[ \frac{m/s}{\sqrt{Hz}} \right]
\end{equation}

#+begin_src matlab
  figure;
  loglog(f, sqrt(Sv));
  xlabel('Frequency [Hz]');
  ylabel('Amplitude Spectral Density $\left[\frac{m/s}{\sqrt{Hz}}\right]$')
#+end_src

#+NAME: fig:asd_velocity
#+HEADER: :tangle no :exports results :results value raw replace :noweb yes
#+begin_src matlab :var filepath="figs/asd_velocity.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
  <<plt-matlab>>
#+end_src

#+NAME: fig:asd_velocity
#+CAPTION: Power Spectral Density of the measured velocity
#+RESULTS: fig:asd_velocity

* Modification of a signal's Power Spectral Density when going through an LTI system

#+begin_src latex :file velocity_to_voltage_psd.pdf :post pdf2svg(file=*this*, ext="png") :exports results
  \begin{tikzpicture}
    \node[block] (G) at (0, 0) {$G(s)$};

    \draw[<-] (G.west) -- node[midway, above]{$x$} ++(-1.4, 0);
    \draw[->] (G.east) -- node[midway, above]{$y$} ++(1.4, 0);
  \end{tikzpicture}
#+end_src

#+NAME: fig:velocity_to_voltage_psd
#+CAPTION: Schematic of the instrumentation used for the measurement
#+RESULTS: fig:velocity_to_voltage_psd
[[file:figs/velocity_to_voltage_psd.png]]

We can show that:
\begin{equation}
  S_{yy}(\omega) = \left|G(j\omega)\right|^2 S_{xx}(\omega)
\end{equation}

And we also have:
\begin{equation}
  \Gamma_{yy}(\omega) = \left|G(j\omega)\right| \Gamma_{xx}(\omega)
\end{equation}

* From PSD of the velocity to the PSD of the displacement

#+begin_src latex :file velocity_to_displacement_psd.pdf :post pdf2svg(file=*this*, ext="png") :exports results
  \begin{tikzpicture}
    \node[block] (G) at (0, 0) {$\frac{1}{s}$};

    \draw[<-] (G.west) -- node[midway, above]{$v$} ++(-1, 0);
    \draw[->] (G.east) -- node[midway, above]{$x$} ++(1, 0);
  \end{tikzpicture}
#+end_src

#+NAME: fig:velocity_to_displacement_psd
#+CAPTION: Schematic of the instrumentation used for the measurement
#+RESULTS: fig:velocity_to_displacement_psd
[[file:figs/velocity_to_displacement_psd.png]]

The displacement is the integral of the velocity.

We then have that
\begin{equation}
  S_{xx}(\omega) = \left|\frac{1}{j \omega}\right|^2 S_{vv}(\omega)
\end{equation}

Using a frequency variable in Hz:
\begin{equation}
  S_{xx}(f) = \left| \frac{1}{j 2\pi f} \right|^2 S_{vv}(f)
\end{equation}

For the Amplitude Spectral Density:
\begin{equation}
  \Gamma_{xx}(f) = \frac{1}{2\pi f} \Gamma_{vv}(f)
\end{equation}

#+begin_note
  \begin{equation}
    S_{xx}(\omega = 1) = S_{vv}(\omega = 1)
  \end{equation}
#+end_note




Now if we want to obtain the Power Spectral Density of the Position or Acceleration:
For each frequency:
\[ \left| \frac{d sin(2 \pi f t)}{dt} \right| = | 2 \pi f | \times | \cos(2\pi f t) | \]

\[ \left| \int_0^t sin(2 \pi f \tau) d\tau \right| = \left| \frac{1}{2 \pi f} \right| \times | \cos(2\pi f t) | \]

\[ ASD_x(f) = \frac{1}{2\pi f} ASD_v(f) \ \left[\frac{m}{\sqrt{Hz}}\right] \]

\[ ASD_a(f) = 2\pi f ASD_v(f) \ \left[\frac{m/s^2}{\sqrt{Hz}}\right] \]
And we have
\[ PSD_x(f) = {ASD_x(f)}^2 = \frac{1}{(2 \pi f)^2} {ASD_v(f)}^2 = \frac{1}{(2 \pi f)^2} PSD_v(f) \]

Note here that we always have
\[ PSD_x \left(f = \frac{1}{2\pi}\right) = PSD_v \left(f = \frac{1}{2\pi}\right) = PSD_a \left(f = \frac{1}{2\pi}\right), \quad \frac{1}{2\pi} \approx 0.16 [Hz] \]

If we want to compute the Cumulative Power Spectrum:
\[ CPS_v(f) = \int_0^f PSD_v(\nu) d\nu \quad [(m/s)^2] \]

We can also want to integrate from high frequency to low frequency:
\[ CPS_v(f) = \int_f^\infty PSD_v(\nu) d\nu \quad [(m/s)^2] \]

The Cumulative Amplitude Spectrum is then the square root of the Cumulative Power Spectrum:
\[ CAS_v(f) = \sqrt{CPS_v(f)} = \sqrt{\int_f^\infty PSD_v(\nu) d\nu} \quad [m/s] \]

Then, we can obtain the Root Mean Square value of the velocity:
\[ v_{\text{rms}} = CAS_v(0) \quad [m/s \ \text{rms}] \]

#+begin_src matlab

#+end_src

* Bibliography                                                       :ignore:
bibliographystyle:unsrt
bibliography:ref.bib