filters-matlab-bank/index.org

#+TITLE: List of filters - Matlab Implementation
:DRAWER:
#+LANGUAGE: en
#+EMAIL: dehaeze.thomas@gmail.com
#+AUTHOR: Dehaeze Thomas

#+HTML_LINK_HOME: ../index.html
#+HTML_LINK_UP:   ../index.html

#+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/Cloud/tikz/org/}{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:

* 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

* Low Pass
** First Order Low Pass Filter
\[ H(s) = \frac{1}{1 + s/\omega_0} \]

Parameters:
- $\omega_0$: cut-off frequency in [rad/s]

Characteristics:
- Low frequency gain of $1$
- Roll-off equals to -20 dB/dec

Matlab code:
#+begin_src matlab
  w0 = 2*pi; % Cut-off Frequency [rad/s]

  H = 1/(1 + s/w0);
#+end_src

#+begin_src matlab :exports none
  freqs = logspace(-2, 2, 1000);

  figure;
  tiledlayout(2, 1, 'TileSpacing', 'None', 'Padding', 'None');

  % Magnitude
  ax1 = nexttile;
  hold on;
  plot(freqs, abs(squeeze(freqresp(H, freqs, 'Hz'))), 'k-');
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  ylabel('Magnitude'); set(gca, 'XTickLabel',[]);
  axis padded
  hold off;
  xline(1, '--', 'LineWidth', 2);

  % Phase
  ax2 = nexttile;
  hold on;
  plot(freqs, 180/pi*angle(squeeze(freqresp(H, freqs, 'Hz'))), 'k-');
  set(gca,'xscale','log');
  yticks(-90:15:0); ylim([-90 0]);
  xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
  hold off;
  xline(1, '--', {'\omega_0'}, 'LineWidth', 2, 'LabelVerticalAlignment', 'bottom', 'LabelOrientation', 'horizontal');


  linkaxes([ax1,ax2],'x');
  xlim([freqs(1), freqs(end)]);
#+end_src

#+begin_src matlab :tangle no :exports results :results file replace
  exportFig('figs/filter_low_pass_first_order.pdf', 'width', 'wide', 'height', 'tall');
#+end_src

#+name: fig:filter_low_pass_first_order
#+caption:
#+RESULTS:
[[file:figs/filter_low_pass_first_order.png]]


** Second Order
\[ H(s) = \frac{1}{1 + 2 \xi / \omega_0 s + s^2/\omega_0^2} \]

Parameters:
- $\omega_0$:
- $\xi$: Damping ratio

Characteristics:
- Low frequency gain: 1
- High frequency roll off: - 40 dB/dec

Matlab code:
#+begin_src matlab
  w0 = 2*pi; % Cut-off frequency [rad/s]
  xi = 0.3; % Damping Ratio

  H = 1/(1 + 2*xi/w0*s + s^2/w0^2);
#+end_src

#+begin_src matlab :exports none
  xis = [0.1, 0.3, 0.5, 0.7, 0.9];
  Hs = {zeros(length(xis), 1)};

  for xi_i = 1:length(xis)
      Hs(xi_i) = {1/(1 + 2*xis(xi_i)/w0*s + s^2/w0^2)};
  end
#+end_src

#+begin_src matlab :exports none
  freqs = logspace(-2, 2, 1000);

  figure;
  tiledlayout(2, 1, 'TileSpacing', 'None', 'Padding', 'None');

  % Magnitude
  ax1 = nexttile;
  hold on;
  for xi_i = 1:length(xis)
      plot(freqs, abs(squeeze(freqresp(Hs{xi_i}, freqs, 'Hz'))), '-', ...
           'DisplayName', sprintf('$\\xi = %.1f$', xis(xi_i)));
  end
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  ylabel('Magnitude'); set(gca, 'XTickLabel',[]);
  axis padded
  legend('location', 'northeast');
  hold off;
  xline(1, '--', 'LineWidth', 2, 'HandleVisibility', 'off');

  % Phase
  ax2 = nexttile;
  hold on;
  for xi_i = 1:length(xis)
    plot(freqs, 180/pi*angle(squeeze(freqresp(Hs{xi_i}, freqs, 'Hz'))), '-');
  end
  set(gca,'xscale','log');
  yticks(-180:30:0); ylim([-180 0]);
  xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
  hold off;
  xline(1, '--', {'\omega_0'}, 'LineWidth', 2, 'LabelVerticalAlignment', 'bottom', 'LabelOrientation', 'horizontal');


  linkaxes([ax1,ax2],'x');
  xlim([freqs(1), freqs(end)]);
#+end_src

#+begin_src matlab :tangle no :exports results :results file replace
  exportFig('figs/filter_low_pass_second_order.pdf', 'width', 'wide', 'height', 'tall');
#+end_src

#+name: fig:filter_low_pass_second_order
#+caption:
#+RESULTS:
[[file:figs/filter_low_pass_second_order.png]]

** Combine multiple first order filters
\[ H(s) = \left( \frac{1}{1 + s/\omega_0} \right)^n \]

Matlab code:
#+begin_src matlab
  w0 = 2*pi; % Cut-off frequency [rad/s]
  n = 3; % Filter order

  H = (1/(1 + s/w0))^n;
#+end_src

#+begin_src matlab :exports none
  ns = [1, 2, 3, 4];
  Hs = {zeros(length(ns), 1)};

  for n_i = 1:length(ns)
      Hs(n_i) = {1/(1 + s/w0)^ns(n_i)};
  end
#+end_src

#+begin_src matlab :exports none
  freqs = logspace(-2, 2, 1000);

  figure;
  tiledlayout(2, 1, 'TileSpacing', 'None', 'Padding', 'None');

  % Magnitude
  ax1 = nexttile;
  hold on;
  for n_i = 1:length(ns)
      plot(freqs, abs(squeeze(freqresp(Hs{n_i}, freqs, 'Hz'))), '-', ...
           'DisplayName', sprintf('$n = %i$', ns(n_i)));
  end
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  ylabel('Magnitude'); set(gca, 'XTickLabel',[]);
  axis padded
  legend('location', 'northeast');
  hold off;
  xline(1, '--', 'LineWidth', 2, 'HandleVisibility', 'off');

  % Phase
  ax2 = nexttile;
  hold on;
  for n_i = 1:length(ns)
    plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(Hs{n_i}, freqs, 'Hz')))), '-');
  end
  set(gca,'xscale','log');
  yticks(-360:45:0); ylim([-360 0]);
  xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
  hold off;
  xline(1, '--', {'\omega_0'}, 'LineWidth', 2, 'LabelVerticalAlignment', 'bottom', 'LabelOrientation', 'horizontal');


  linkaxes([ax1,ax2],'x');
  xlim([freqs(1), freqs(end)]);
#+end_src

#+begin_src matlab :tangle no :exports results :results file replace
  exportFig('figs/filter_low_pass_first_order_add.pdf', 'width', 'wide', 'height', 'tall');
#+end_src

#+name: fig:filter_low_pass_first_order_add
#+caption:
#+RESULTS:
[[file:figs/filter_low_pass_first_order_add.png]]

* High Pass
** First Order
\[ H(s) = \frac{s/\omega_0}{1 + s/\omega_0} \]

Parameters:
- $\omega_0$: cut-off frequency in [rad/s]

Characteristics:
- High frequency gain of $1$
- Low frequency slow of +20 dB/dec

Matlab code:
#+begin_src matlab
  w0 = 2*pi; % Cut-off frequency [rad/s]

  H = (s/w0)/(1 + s/w0);
#+end_src

#+begin_src matlab :exports none
  freqs = logspace(-2, 2, 1000);

  figure;
  tiledlayout(2, 1, 'TileSpacing', 'None', 'Padding', 'None');

  % Magnitude
  ax1 = nexttile;
  hold on;
  plot(freqs, abs(squeeze(freqresp(H, freqs, 'Hz'))), 'k-');
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  ylabel('Magnitude'); set(gca, 'XTickLabel',[]);
  axis padded
  hold off;
  xline(1, '--', 'LineWidth', 2);

  % Phase
  ax2 = nexttile;
  hold on;
  plot(freqs, 180/pi*angle(squeeze(freqresp(H, freqs, 'Hz'))), 'k-');
  set(gca,'xscale','log');
  yticks(0:15:90); ylim([0, 90]);
  xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
  hold off;
  xline(1, '--', {'\omega_0'}, 'LineWidth', 2, 'LabelVerticalAlignment', 'bottom', 'LabelOrientation', 'horizontal');


  linkaxes([ax1,ax2],'x');
  xlim([freqs(1), freqs(end)]);
#+end_src

#+begin_src matlab :tangle no :exports results :results file replace
  exportFig('figs/filter_high_pass_first_order.pdf', 'width', 'wide', 'height', 'tall');
#+end_src

#+name: fig:filter_high_pass_first_order
#+caption:
#+RESULTS:
[[file:figs/filter_high_pass_first_order.png]]

** Second Order
\[ H(s) = \frac{s^2/\omega_0^2}{1 + 2 \xi / \omega_0 s + s^2/\omega_0^2} \]

Parameters:
- $\omega_0$:
- $\xi$: Damping ratio

Matlab code:
#+begin_src matlab
  w0 = 2*pi; % [rad/s]
  xi = 0.3;

  H = (s^2/w0^2)/(1 + 2*xi/w0*s + s^2/w0^2);
#+end_src

#+begin_src matlab :exports none
  xis = [0.1, 0.3, 0.5, 0.7, 0.9];
  Hs = {zeros(length(xis), 1)};

  for xi_i = 1:length(xis)
      Hs(xi_i) = {(s^2/w0^2)/(1 + 2*xis(xi_i)/w0*s + s^2/w0^2)};
  end

  freqs = logspace(-2, 2, 1000);

  figure;
  tiledlayout(2, 1, 'TileSpacing', 'None', 'Padding', 'None');

  % Magnitude
  ax1 = nexttile;
  hold on;
  for xi_i = 1:length(xis)
      plot(freqs, abs(squeeze(freqresp(Hs{xi_i}, freqs, 'Hz'))), '-', ...
           'DisplayName', sprintf('$\\xi = %.1f$', xis(xi_i)));
  end
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  ylabel('Magnitude'); set(gca, 'XTickLabel',[]);
  axis padded
  legend('location', 'northeast');
  hold off;
  xline(1, '--', 'LineWidth', 2, 'HandleVisibility', 'off');

  % Phase
  ax2 = nexttile;
  hold on;
  for xi_i = 1:length(xis)
    plot(freqs, 180/pi*angle(squeeze(freqresp(Hs{xi_i}, freqs, 'Hz'))), '-');
  end
  set(gca,'xscale','log');
  yticks(0:30:180); ylim([0 180]);
  xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
  hold off;
  xline(1, '--', {'\omega_0'}, 'LineWidth', 2, 'LabelVerticalAlignment', 'bottom', 'LabelOrientation', 'horizontal');


  linkaxes([ax1,ax2],'x');
  xlim([freqs(1), freqs(end)]);
#+end_src

#+begin_src matlab :tangle no :exports results :results file replace
  exportFig('figs/filter_high_pass_second_order.pdf', 'width', 'wide', 'height', 'tall');
#+end_src

#+name: fig:filter_high_pass_second_order
#+caption:
#+RESULTS:
[[file:figs/filter_high_pass_second_order.png]]

** Combine multiple filters
\[ H(s) = \left( \frac{s/\omega_0}{1 + s/\omega_0} \right)^n \]

Matlab code:
#+begin_src matlab
  w0 = 2*pi; % [rad/s]
  n = 3;

  H = ((s/w0)/(1 + s/w0))^n;
#+end_src

#+begin_src matlab :exports none
  ns = [1, 2, 3, 4];
  Hs = {zeros(length(ns), 1)};

  for n_i = 1:length(ns)
      Hs(n_i) = {((s/w0)/(1 + s/w0))^ns(n_i)};
  end
#+end_src

#+begin_src matlab :exports none
  freqs = logspace(-2, 2, 1000);

  figure;
  tiledlayout(2, 1, 'TileSpacing', 'None', 'Padding', 'None');

  % Magnitude
  ax1 = nexttile;
  hold on;
  for n_i = 1:length(ns)
      plot(freqs, abs(squeeze(freqresp(Hs{n_i}, freqs, 'Hz'))), '-', ...
           'DisplayName', sprintf('$n = %i$', ns(n_i)));
  end
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  ylabel('Magnitude'); set(gca, 'XTickLabel',[]);
  axis padded
  legend('location', 'northeast');
  hold off;
  xline(1, '--', 'LineWidth', 2, 'HandleVisibility', 'off');

  % Phase
  ax2 = nexttile;
  hold on;
  for n_i = 1:length(ns)
    plot(freqs, 180/pi*angle(squeeze(freqresp(Hs{n_i}, freqs, 'Hz'))), '-');
  end
  set(gca,'xscale','log');
  yticks(-180:45:180); ylim([-180 180]);
  xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
  hold off;
  xline(1, '--', {'\omega_0'}, 'LineWidth', 2, 'LabelVerticalAlignment', 'bottom', 'LabelOrientation', 'horizontal');

  linkaxes([ax1,ax2],'x');
  xlim([freqs(1), freqs(end)]);
#+end_src

#+begin_src matlab :tangle no :exports results :results file replace
  exportFig('figs/filter_high_pass_first_order_add.pdf', 'width', 'wide', 'height', 'tall');
#+end_src

#+name: fig:filter_high_pass_first_order_add
#+caption:
#+RESULTS:
[[file:figs/filter_high_pass_first_order_add.png]]

* TODO Band Pass
** Second Order

* Notch
** Second Order
\begin{equation}
  \frac{s^2 + 2 g_c \xi \omega_n s + \omega_n^2}{s^2 + 2 \xi \omega_n s + \omega_n^2}
\end{equation}

Parameters:
- $\omega_n$: frequency of the notch
- $g_c$: gain at the notch frequency
- $\xi$: damping ratio (notch width)

Matlab code:
#+begin_src matlab
  gc = 0.02;
  xi = 0.1;
  wn = 2*pi;

  H = (s^2 + 2*gm*xi*wn*s + wn^2)/(s^2 + 2*xi*wn*s + wn^2);
#+end_src

#+begin_src matlab :exports none
  xi = 0.2;
  gcs = [0.02, 0.1, 0.5];
  Hs = {zeros(length(gcs), 1)};

  for g_i = 1:length(gcs)
      Hs(g_i) = {(s^2 + 2*gcs(g_i)*xi*wn*s + wn^2)/(s^2 + 2*xi*wn*s + wn^2)};
  end

  freqs = logspace(-1, 1, 1000);

  figure;
  tiledlayout(2, 1, 'TileSpacing', 'None', 'Padding', 'None');

  % Magnitude
  ax1 = nexttile;
  hold on;
  for g_i = 1:length(gcs)
      plot(freqs, abs(squeeze(freqresp(Hs{g_i}, freqs, 'Hz'))), '-', ...
           'DisplayName', sprintf('$g_c = %.2f$', gm(g_i)));
  end
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  ylabel('Magnitude'); set(gca, 'XTickLabel',[]);
  axis padded
  legend('location', 'southeast');
  hold off;

  % Phase
  ax2 = nexttile;
  hold on;
  for g_i = 1:length(gcs)
    plot(freqs, 180/pi*angle(squeeze(freqresp(Hs{g_i}, freqs, 'Hz'))), '-');
  end
  set(gca,'xscale','log');
  yticks(-90:30:90); ylim([-90 90]);
  xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
  hold off;

  linkaxes([ax1,ax2],'x');
  xlim([freqs(1), freqs(end)]);
#+end_src

#+begin_src matlab :tangle no :exports results :results file replace
  exportFig('figs/filter_notch_xi.pdf', 'width', 'wide', 'height', 'tall');
#+end_src

#+name: fig:filter_notch_xi
#+caption:
#+RESULTS:
[[file:figs/filter_notch_xi.png]]

#+begin_src matlab :exports none
  xis = [0.05, 0.1, 0.5];
  gc = 0.1;
  Hs = {zeros(length(xis), 1)};

  for xi_i = 1:length(xis)
      Hs(xi_i) = {(s^2 + 2*gc*xis(xi_i)*wn*s + wn^2)/(s^2 + 2*xis(xi_i)*wn*s + wn^2)};
  end

  freqs = logspace(-1, 1, 1000);

  figure;
  tiledlayout(2, 1, 'TileSpacing', 'None', 'Padding', 'None');

  % Magnitude
  ax1 = nexttile;
  hold on;
  for xi_i = 1:length(xis)
      plot(freqs, abs(squeeze(freqresp(Hs{xi_i}, freqs, 'Hz'))), '-', ...
           'DisplayName', sprintf('$\\xi = %.2f$', xis(xi_i)));
  end
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  ylabel('Magnitude'); set(gca, 'XTickLabel',[]);
  axis padded
  legend('location', 'southeast');
  hold off;

  % Phase
  ax2 = nexttile;
  hold on;
  for xi_i = 1:length(xis)
    plot(freqs, 180/pi*angle(squeeze(freqresp(Hs{xi_i}, freqs, 'Hz'))), '-');
  end
  set(gca,'xscale','log');
  yticks(-90:30:90); ylim([-90 90]);
  xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
  hold off;

  linkaxes([ax1,ax2],'x');
  xlim([freqs(1), freqs(end)]);
#+end_src

#+begin_src matlab :tangle no :exports results :results file replace
  exportFig('figs/filter_notch_gc.pdf', 'width', 'wide', 'height', 'tall');
#+end_src

#+name: fig:filter_notch_gc
#+caption:
#+RESULTS:
[[file:figs/filter_notch_gc.png]]

* TODO Chebyshev
** Chebyshev Type I

#+begin_src matlab
  n = 4; % Order of the filter
  Rp = 3; % Maximum peak-to-peak ripple [dB]
  Wp = 2*pi; % passband-edge frequency [rad/s]

  [A,B,C,D] = cheby1(n, Rp, Wp, 'high', 's');
  H = ss(A, B, C, D);
#+end_src

* Lead - Lag
** Lead
\begin{equation}
  H(s) = \frac{1 + \frac{s}{w_c/\sqrt{a}}}{1 + \frac{s}{w_c \sqrt{a}}}, \quad a > 1
\end{equation}

Parameters:
- $\omega_c$: frequency at which the phase lead is maximum
- $a$: parameter to adjust the phase lead, also impacts the high frequency gain

Characteristics:
- the low frequency gain is $1$
- the high frequency gain is $a$
- the phase lead at $\omega_c$ is equal to (Figure [[fig:filter_lead_effect_a_phase]]):
  \[ \angle H(j\omega_c) = \tan^{-1}(\sqrt{a}) - \tan^{-1}(1/\sqrt{a}) \]

Matlab code:
#+begin_src matlab
  a  = 0.6;  % Amount of phase lead / width of the phase lead / high frequency gain
  wc = 2*pi; % Frequency with the maximum phase lead [rad/s]

  H = (1 + s/(wc/sqrt(a)))/(1 + s/(wc*sqrt(a)));
#+end_src

#+begin_src matlab :exports none
  wc = 2*pi;
  as = [1.2, 2, 5];
  Hs = {zeros(length(as), 1)};

  for a_i = 1:length(as)
      Hs(a_i) = {(1 + s/(wc/sqrt(as(a_i))))/(1 + s/(wc*sqrt(as(a_i))))};
  end

  freqs = logspace(-2, 2, 1000);

  figure;
  tiledlayout(2, 1, 'TileSpacing', 'None', 'Padding', 'None');

  % Magnitude
  ax1 = nexttile;
  hold on;
  for a_i = 1:length(as)
      plot(freqs, abs(squeeze(freqresp(Hs{a_i}, freqs, 'Hz'))), '-', ...
           'DisplayName', sprintf('$a = %.1f$', as(a_i)));
  end
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  ylabel('Magnitude'); set(gca, 'XTickLabel',[]);
  axis padded
  legend('location', 'northwest');
  hold off;
  xline(1, '--', 'LineWidth', 2, 'HandleVisibility', 'off');

  % Phase
  ax2 = nexttile;
  hold on;
  for a_i = 1:length(as)
    plot(freqs, 180/pi*angle(squeeze(freqresp(Hs{a_i}, freqs, 'Hz'))), '-');
  end
  set(gca,'xscale','log');
  yticks(0:15:60); ylim([0 60]);
  xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
  hold off;
  xline(1, '--', {'\omega_c'}, 'LineWidth', 2, 'LabelVerticalAlignment', 'top', 'LabelOrientation', 'horizontal');


  linkaxes([ax1,ax2],'x');
  xlim([freqs(1), freqs(end)]);
#+end_src

#+begin_src matlab :tangle no :exports results :results file replace
  exportFig('figs/filter_lead.pdf', 'width', 'wide', 'height', 'tall');
#+end_src

#+name: fig:filter_lead
#+caption:
#+RESULTS:
[[file:figs/filter_lead.png]]

#+begin_src matlab :exports none
  as = logspace(0, 1, 100);

  figure;
  plot(as, 180/pi*(atan(sqrt(as)) - atan(1./sqrt(as))))
  xlabel('$a$'); ylabel('Phase added at $\omega_c$ [deg]');
  yticks(0:15:60); ylim([0, 60]);
  xticks(1:1:10); xlim([as(1), as(end)]);
#+end_src

#+begin_src matlab :tangle no :exports results :results file replace
  exportFig('figs/filter_lead_effect_a_phase.pdf', 'width', 'normal', 'height', 'normal');
#+end_src

#+name: fig:filter_lead_effect_a_phase
#+caption:
#+RESULTS:
[[file:figs/filter_lead_effect_a_phase.png]]

** Lag
\begin{equation}
  H(s) = \frac{w_c \sqrt{a} + s}{\frac{w_c}{\sqrt{a}} + s}, \quad a > 1
\end{equation}

Parameters:
- $\omega_c$: frequency at which the phase lag is maximum
- $a$: parameter to adjust the phase lag, also impacts the low frequency gain

Characteristics:
- the low frequency gain is increased by a factor $a$
- the high frequency gain is $1$ (unchanged)
- the phase lag at $\omega_c$ is equal to (Figure [[fig:filter_lag_effect_a_phase]]):
  \[ \angle H(j\omega_c) = \tan^{-1}(1/\sqrt{a}) - \tan^{-1}(\sqrt{a}) \]

Matlab code:
#+begin_src matlab
  a  = 0.6;  % Amount of phase lag / width of the phase lag / high frequency gain
  wc = 2*pi; % Frequency with the maximum phase lag [rad/s]

  H = (wc*sqrt(a) + s)/(wc/sqrt(a) + s);
#+end_src

#+begin_src matlab :exports none
  wc = 2*pi;
  as = [1.2, 2, 5];
  Hs = {zeros(length(as), 1)};

  for a_i = 1:length(as)
      Hs(a_i) = {(wc*sqrt(as(a_i)) + s)/(wc/sqrt(as(a_i)) + s)};
  end

  freqs = logspace(-2, 2, 1000);

  figure;
  tiledlayout(2, 1, 'TileSpacing', 'None', 'Padding', 'None');

  % Magnitude
  ax1 = nexttile;
  hold on;
  for a_i = 1:length(as)
      plot(freqs, abs(squeeze(freqresp(Hs{a_i}, freqs, 'Hz'))), '-', ...
           'DisplayName', sprintf('$a = %.1f$', as(a_i)));
  end
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  ylabel('Magnitude'); set(gca, 'XTickLabel',[]);
  axis padded
  legend('location', 'southwest');
  hold off;
  xline(1, '--', 'LineWidth', 2, 'HandleVisibility', 'off');

  % Phase
  ax2 = nexttile;
  hold on;
  for a_i = 1:length(as)
    plot(freqs, 180/pi*angle(squeeze(freqresp(Hs{a_i}, freqs, 'Hz'))), '-');
  end
  set(gca,'xscale','log');
  yticks(-60:15:0); ylim([-60 0]);
  xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
  hold off;
  xline(1, '--', {'\omega_c'}, 'LineWidth', 2, 'LabelVerticalAlignment', 'bottom', 'LabelOrientation', 'horizontal');


  linkaxes([ax1,ax2],'x');
  xlim([freqs(1), freqs(end)]);
#+end_src

#+begin_src matlab :tangle no :exports results :results file replace
  exportFig('figs/filter_lag.pdf', 'width', 'wide', 'height', 'tall');
#+end_src

#+name: fig:filter_lag
#+caption:
#+RESULTS:
[[file:figs/filter_lag.png]]

#+begin_src matlab :exports none
  as = logspace(0, 1, 100);

  figure;
  plot(as, 180/pi*(atan(1./sqrt(as)) - atan(sqrt(as))))
  xlabel('$a$'); ylabel('Phase added at $\omega_c$ [deg]');
  yticks(-60:15:0); ylim([-60, 0]);
  xticks(1:1:10); xlim([as(1), as(end)]);
#+end_src

#+begin_src matlab :tangle no :exports results :results file replace
  exportFig('figs/filter_lag_effect_a_phase.pdf', 'width', 'normal', 'height', 'normal');
#+end_src

#+name: fig:filter_lag_effect_a_phase
#+caption:
#+RESULTS:
[[file:figs/filter_lag_effect_a_phase.png]]

* TODO Complementary

* TODO Performance Weight
** Nice combination
\begin{equation}
  W(s) = G_c * \left(\frac{\frac{1}{\omega_0}\sqrt{\frac{1 - \left(\frac{G_0}{G_c}\right)^{\frac{2}{n}}}{1 - \left(\frac{G_c}{G_\infty}\right)^{\frac{2}{n}}}} s + \left(\frac{G_0}{G_c}\right)^{\frac{1}{n}}}{\frac{1}{\omega_0} \sqrt{\frac{1 - \left(\frac{G_0}{G_c}\right)^{\frac{2}{n}}}{\left(\frac{G_\infty}{G_c}\right)^{\frac{2}{n}} - 1}} s + 1}\right)^n
\end{equation}


#+begin_src matlab
  n = 2; w0 = 2*pi*11; G0 = 1/10; G1 = 1000; Gc = 1/2;
  wL = Gc*(((G1/Gc)^(1/n)/w0*sqrt((1-(G0/Gc)^(2/n))/((G1/Gc)^(2/n)-1))*s + (G0/Gc)^(1/n))/(1/w0*sqrt((1-(G0/Gc)^(2/n))/((G1/Gc)^(2/n)-1))*s + 1))^n;

  n = 3; w0 = 2*pi*9; G0 = 10000; G1 = 0.1; Gc = 1/2;
  wH = Gc*(((G1/Gc)^(1/n)/w0*sqrt((1-(G0/Gc)^(2/n))/((G1/Gc)^(2/n)-1))*s + (G0/Gc)^(1/n))/(1/w0*sqrt((1-(G0/Gc)^(2/n))/((G1/Gc)^(2/n)-1))*s + 1))^n;
#+end_src


** Alternative
#+begin_src matlab
  w0 = 2*pi; % [rad/s]
  A = 1e-2;
  M = 5;

  H = (s/sqrt(M) + w0)^2/(s + w0*sqrt(A))^2;
#+end_src

#+begin_src matlab :exports none
  freqs = logspace(-2, 2, 1000);
  resp = squeeze(freqresp(inv(H), freqs, 'Hz'));
  Ha = abs(resp);
  Hp = 180/pi*phase(resp);

  T = table(freqs', Ha, Hp, 'VariableNames', {'freqs', 'amplitude', 'phase'});
  writetable(T,'mat/weight_first_order.csv');
#+end_src

#+begin_src latex :file weight_first_order.pdf :tangle figs/weight_first_order.tex :exports results
  \setlength\fwidth{8cm}
  \setlength\fheight{4cm}

  \definecolor{mycolor1}{rgb}{0.00000,0.44700,0.74100}%
  \definecolor{mycolor2}{rgb}{0.85000,0.32500,0.09800}%

  \begin{tikzpicture}
    \begin{axis}[%
      width=\fwidth,
      height=\fheight,
      at={(0,0)},
      scale only axis,
      separate axis lines,
      every outer x axis line/.append style={black},
      every x tick label/.append style={font=\color{black}},
      every x tick/.append style={black},
      xmode=log,
      xmin=0.01,
      xmax=100,
      xminorticks=true,
      xlabel={Frequency [Hz]},
      every outer y axis line/.append style={black},
      every y tick label/.append style={font=\color{black}},
      every y tick/.append style={black},
      ymode=log,
      ymin=5e-3,
      ymax=1e1,
      yminorticks=true,
      ylabel={Amplitude},
      axis background/.style={fill=white},
      xmajorgrids,
      xminorgrids,
      ymajorgrids,
      yminorgrids
      ]
      \addplot[color=black, mark=none]
      table [x=freqs, y=amplitude, col sep=comma] {/home/thomas/Cloud/thesis/matlab/filters/matweight_first_order.csv};

      \draw[dashed] (0.01,1e-2) -- (1,1e-2) node[right]{$A$};
      \draw[dashed] (1,5) node[left]{$M$} -- (100,5);
      \draw[dashed] (1,1) -- (1,0.5) node[below]{$\omega_b^*$};
    \end{axis}
  \end{tikzpicture}
#+end_src

#+RESULTS:
[[file:figs/weight_first_order.png]]

* TODO Combine Filters
** Additive
- [ ] Explain how phase and magnitude combine

** Multiplicative
* TODO Filters representing noise
Let's consider a noise $n$ that is shaped from a white-noise $\tilde{n}$ with unitary PSD ($\Phi_\tilde{n}(\omega) = 1$) using a transfer function $G(s)$.
The PSD of $n$ is then:
\[ \Phi_n(\omega) = |G(j\omega)|^2 \Phi_{\tilde{n}}(\omega) = |G(j\omega)|^2 \]

The PSD $\Phi_n(\omega)$ is expressed in $\text{unit}^2/\text{Hz}$.

And the root mean square (RMS) of $n(t)$ is:
\[ \sigma_n = \sqrt{\int_{0}^{\infty} \Phi_n(\omega) d\omega} \]

** First Order Low Pass Filter
\[ G(s) = \frac{g_0}{1 + \frac{s}{\omega_c}} \]

#+begin_src matlab
  g0 = 1; % Noise Density in unit/sqrt(Hz)
  wc = 1; % Cut-Off frequency [rad/s]

  G = g0/(1 + s/wc);

  % Frequency vector [Hz]
  freqs = logspace(-3, 3, 1000);

  % PSD of n in [unit^2/Hz]
  Phi_n = abs(squeeze(freqresp(G, freqs, 'Hz'))).^2;

  % RMS value of n in [unit, rms]
  sigma_n = sqrt(trapz(freqs, Phi_n))
#+end_src

\[ \sigma = \frac{1}{2} g_0 \sqrt{\omega_c} \]
with:
- $g_0$ the Noise Density of $n$ in $\text{unit}/\sqrt{Hz}$
- $\omega_c$ the bandwidth over which the noise is located, in rad/s
- $\sigma$ the rms noise

If the cut-off frequency is to be expressed in Hz:
\[ \sigma = \frac{1}{2} g_0 \sqrt{2\pi f_c} = \sqrt{\frac{\pi}{2}} g_0 \sqrt{f_c} \]


Thus, if a sensor is said to have a RMS noise of $\sigma = 10 nm\ rms$ over a bandwidth of $\omega_c = 100 rad/s$, we can estimated the noise density of the sensor to be (supposing a first order low pass filter noise shape):
\[ g_0 = \frac{2 \sigma}{\sqrt{\omega_c}} \quad \left[ m/\sqrt{Hz} \right] \]

#+begin_src matlab :results value replace
  2*10e-9/sqrt(100)
#+end_src

#+RESULTS:
: 2e-09

#+begin_src matlab :results value replace
  6*0.5*20e-12*sqrt(2*pi*100)
#+end_src

#+RESULTS:
: 1.504e-09