From ea35194620392845b6746a9fd08f81933a124287 Mon Sep 17 00:00:00 2001 From: Thomas Dehaeze Date: Thu, 29 Apr 2021 15:03:59 +0200 Subject: [PATCH] Add functions to generate weights and CF --- matlab/index.org | 1394 +++++++++++++++++++------------- matlab/matlab/src/generateCF.m | 28 + matlab/matlab/src/generateWF.m | 39 + 3 files changed, 915 insertions(+), 546 deletions(-) create mode 100644 matlab/matlab/src/generateCF.m create mode 100644 matlab/matlab/src/generateWF.m diff --git a/matlab/index.org b/matlab/index.org index 8930375..f8ef0b7 100644 --- a/matlab/index.org +++ b/matlab/index.org @@ -65,15 +65,24 @@ This document is divided into several sections: ** Matlab Init :noexport:ignore: #+begin_src matlab :tangle no :exports none :results silent :noweb yes :var current_dir=(file-name-directory buffer-file-name) - <> +<> #+end_src #+begin_src matlab :exports none :results silent :noweb yes - <> +<> #+end_src #+begin_src matlab - freqs = logspace(-1, 3, 1000); +freqs = logspace(-1, 3, 1000); +#+end_src + +#+begin_src matlab :tangle no +addpath('./matlab'); +addpath('./matlab/src'); +#+end_src + +#+begin_src matlab :exec no +addpath('./src'); #+end_src ** Synthesis Architecture @@ -118,38 +127,38 @@ The parameters permits to specify: The general shape of a weighting function generated using the formula is shown in figure [[fig:weight_formula]]. #+begin_src matlab :exports none - n = 3; w0 = 2*pi*10; G0 = 1e-3; G1 = 1e1; Gc = 2; +n = 3; w0 = 2*pi*10; G0 = 1e-3; G1 = 1e1; Gc = 2; - W = (((1/w0)*sqrt((1-(G0/Gc)^(2/n))/(1-(Gc/G1)^(2/n)))*s + (G0/Gc)^(1/n))/((1/G1)^(1/n)*(1/w0)*sqrt((1-(G0/Gc)^(2/n))/(1-(Gc/G1)^(2/n)))*s + (1/Gc)^(1/n)))^n; +W = (((1/w0)*sqrt((1-(G0/Gc)^(2/n))/(1-(Gc/G1)^(2/n)))*s + (G0/Gc)^(1/n))/((1/G1)^(1/n)*(1/w0)*sqrt((1-(G0/Gc)^(2/n))/(1-(Gc/G1)^(2/n)))*s + (1/Gc)^(1/n)))^n; - figure; - hold on; - plot(freqs, abs(squeeze(freqresp(W, freqs, 'Hz'))), 'k-'); +figure; +hold on; +plot(freqs, abs(squeeze(freqresp(W, freqs, 'Hz'))), 'k-'); - plot([1e-3 1e0], [G0 G0], 'k--', 'LineWidth', 1) - text(1e0, G0, '$\quad G_0$') +plot([1e-3 1e0], [G0 G0], 'k--', 'LineWidth', 1) +text(1e0, G0, '$\quad G_0$') - plot([1e1 1e3], [G1 G1], 'k--', 'LineWidth', 1) - text(1e1,G1,'$G_{\infty}\quad$','HorizontalAlignment', 'right') +plot([1e1 1e3], [G1 G1], 'k--', 'LineWidth', 1) +text(1e1,G1,'$G_{\infty}\quad$','HorizontalAlignment', 'right') - plot([w0/2/pi w0/2/pi], [1 2*Gc], 'k--', 'LineWidth', 1) - text(w0/2/pi,1,'$\omega_c$','VerticalAlignment', 'top', 'HorizontalAlignment', 'center') +plot([w0/2/pi w0/2/pi], [1 2*Gc], 'k--', 'LineWidth', 1) +text(w0/2/pi,1,'$\omega_c$','VerticalAlignment', 'top', 'HorizontalAlignment', 'center') - plot([w0/2/pi/2 2*w0/2/pi], [Gc Gc], 'k--', 'LineWidth', 1) - text(w0/2/pi/2, Gc, '$G_c \quad$','HorizontalAlignment', 'right') +plot([w0/2/pi/2 2*w0/2/pi], [Gc Gc], 'k--', 'LineWidth', 1) +text(w0/2/pi/2, Gc, '$G_c \quad$','HorizontalAlignment', 'right') - text(w0/5/pi/2, abs(evalfr(W, j*w0/5)), 'Slope: $n \quad$', 'HorizontalAlignment', 'right') +text(w0/5/pi/2, abs(evalfr(W, j*w0/5)), 'Slope: $n \quad$', 'HorizontalAlignment', 'right') - text(w0/2/pi, abs(evalfr(W, j*w0)), '$\bullet$', 'HorizontalAlignment', 'center') - set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); - xlabel('Frequency [Hz]'); ylabel('Magnitude'); - hold off; - xlim([freqs(1), freqs(end)]); - ylim([5e-4, 20]); +text(w0/2/pi, abs(evalfr(W, j*w0)), '$\bullet$', 'HorizontalAlignment', 'center') +set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); +xlabel('Frequency [Hz]'); ylabel('Magnitude'); +hold off; +xlim([freqs(1), freqs(end)]); +ylim([5e-4, 20]); #+end_src #+begin_src matlab :tangle no :exports results :results file replace - exportFig('figs/weight_formula.pdf', 'width', 'wide', 'height', 'normal'); +exportFig('figs/weight_formula.pdf', 'width', 'wide', 'height', 'normal'); #+end_src #+name: fig:weight_formula @@ -158,31 +167,31 @@ The general shape of a weighting function generated using the formula is shown i [[file:figs/weight_formula.png]] #+begin_src matlab - n = 2; w0 = 2*pi*11; G0 = 1/10; G1 = 1000; Gc = 1/2; - W1 = (((1/w0)*sqrt((1-(G0/Gc)^(2/n))/(1-(Gc/G1)^(2/n)))*s + (G0/Gc)^(1/n))/((1/G1)^(1/n)*(1/w0)*sqrt((1-(G0/Gc)^(2/n))/(1-(Gc/G1)^(2/n)))*s + (1/Gc)^(1/n)))^n; +n = 2; w0 = 2*pi*11; G0 = 1/10; G1 = 1000; Gc = 1/2; +W1 = (((1/w0)*sqrt((1-(G0/Gc)^(2/n))/(1-(Gc/G1)^(2/n)))*s + (G0/Gc)^(1/n))/((1/G1)^(1/n)*(1/w0)*sqrt((1-(G0/Gc)^(2/n))/(1-(Gc/G1)^(2/n)))*s + (1/Gc)^(1/n)))^n; - n = 3; w0 = 2*pi*10; G0 = 1000; G1 = 0.1; Gc = 1/2; - W2 = (((1/w0)*sqrt((1-(G0/Gc)^(2/n))/(1-(Gc/G1)^(2/n)))*s + (G0/Gc)^(1/n))/((1/G1)^(1/n)*(1/w0)*sqrt((1-(G0/Gc)^(2/n))/(1-(Gc/G1)^(2/n)))*s + (1/Gc)^(1/n)))^n; +n = 3; w0 = 2*pi*10; G0 = 1000; G1 = 0.1; Gc = 1/2; +W2 = (((1/w0)*sqrt((1-(G0/Gc)^(2/n))/(1-(Gc/G1)^(2/n)))*s + (G0/Gc)^(1/n))/((1/G1)^(1/n)*(1/w0)*sqrt((1-(G0/Gc)^(2/n))/(1-(Gc/G1)^(2/n)))*s + (1/Gc)^(1/n)))^n; #+end_src #+begin_src matlab :exports none - figure; - hold on; - set(gca,'ColorOrderIndex',1) - plot(freqs, 1./abs(squeeze(freqresp(W1, freqs, 'Hz'))), '--', 'DisplayName', '$|W_1|^{-1}$'); - set(gca,'ColorOrderIndex',2) - plot(freqs, 1./abs(squeeze(freqresp(W2, freqs, 'Hz'))), '--', 'DisplayName', '$|W_2|^{-1}$'); - set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); - xlabel('Frequency [Hz]'); ylabel('Magnitude'); - hold off; - xlim([freqs(1), freqs(end)]); - ylim([5e-4, 20]); - xticks([0.1, 1, 10, 100, 1000]); - legend('location', 'northeast', 'FontSize', 8); +figure; +hold on; +set(gca,'ColorOrderIndex',1) +plot(freqs, 1./abs(squeeze(freqresp(W1, freqs, 'Hz'))), '--', 'DisplayName', '$|W_1|^{-1}$'); +set(gca,'ColorOrderIndex',2) +plot(freqs, 1./abs(squeeze(freqresp(W2, freqs, 'Hz'))), '--', 'DisplayName', '$|W_2|^{-1}$'); +set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); +xlabel('Frequency [Hz]'); ylabel('Magnitude'); +hold off; +xlim([freqs(1), freqs(end)]); +ylim([5e-4, 20]); +xticks([0.1, 1, 10, 100, 1000]); +legend('location', 'northeast', 'FontSize', 8); #+end_src #+begin_src matlab :tangle no :exports results :results file replace - exportFig('figs/weights_W1_W2.pdf', 'width', 'wide', 'height', 'normal'); +exportFig('figs/weights_W1_W2.pdf', 'width', 'wide', 'height', 'normal'); #+end_src #+name: fig:weights_W1_W2 @@ -193,14 +202,14 @@ The general shape of a weighting function generated using the formula is shown i ** H-Infinity Synthesis We define the generalized plant $P$ on matlab. #+begin_src matlab - P = [W1 -W1; - 0 W2; - 1 0]; +P = [W1 -W1; + 0 W2; + 1 0]; #+end_src And we do the $\mathcal{H}_\infty$ synthesis using the =hinfsyn= command. #+begin_src matlab :results output replace :exports both - [H2, ~, gamma, ~] = hinfsyn(P, 1, 1,'TOLGAM', 0.001, 'METHOD', 'ric', 'DISPLAY', 'on'); +[H2, ~, gamma, ~] = hinfsyn(P, 1, 1,'TOLGAM', 0.001, 'METHOD', 'ric', 'DISPLAY', 'on'); #+end_src #+RESULTS: @@ -240,55 +249,55 @@ Test bounds: 0.1000 < gamma <= 1050.0000 We then define the high pass filter $H_1 = 1 - H_2$. The bode plot of both $H_1$ and $H_2$ is shown on figure [[fig:hinf_filters_results]]. #+begin_src matlab - H1 = 1 - H2; +H1 = 1 - H2; #+end_src ** Obtained Complementary Filters The obtained complementary filters are shown on figure [[fig:hinf_filters_results]]. #+begin_src matlab :exports none - figure; - tiledlayout(3, 1, 'TileSpacing', 'None', 'Padding', 'None'); +figure; +tiledlayout(3, 1, 'TileSpacing', 'None', 'Padding', 'None'); - % Magnitude - ax1 = nexttile([2, 1]); - hold on; - set(gca,'ColorOrderIndex',1) - plot(freqs, 1./abs(squeeze(freqresp(W1, freqs, 'Hz'))), '--', 'DisplayName', '$w_1$'); - set(gca,'ColorOrderIndex',2) - plot(freqs, 1./abs(squeeze(freqresp(W2, freqs, 'Hz'))), '--', 'DisplayName', '$w_2$'); +% Magnitude +ax1 = nexttile([2, 1]); +hold on; +set(gca,'ColorOrderIndex',1) +plot(freqs, 1./abs(squeeze(freqresp(W1, freqs, 'Hz'))), '--', 'DisplayName', '$w_1$'); +set(gca,'ColorOrderIndex',2) +plot(freqs, 1./abs(squeeze(freqresp(W2, freqs, 'Hz'))), '--', 'DisplayName', '$w_2$'); - set(gca,'ColorOrderIndex',1) - plot(freqs, abs(squeeze(freqresp(H1, freqs, 'Hz'))), '-', 'DisplayName', '$H_1$'); - set(gca,'ColorOrderIndex',2) - plot(freqs, abs(squeeze(freqresp(H2, freqs, 'Hz'))), '-', 'DisplayName', '$H_2$'); - set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); - hold off; - set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); - ylabel('Magnitude'); - set(gca, 'XTickLabel',[]); - ylim([1e-4, 20]); - yticks([1e-4, 1e-3, 1e-2, 1e-1, 1, 1e1]); - legend('location', 'southeast', 'FontSize', 8, 'NumColumns', 2); +set(gca,'ColorOrderIndex',1) +plot(freqs, abs(squeeze(freqresp(H1, freqs, 'Hz'))), '-', 'DisplayName', '$H_1$'); +set(gca,'ColorOrderIndex',2) +plot(freqs, abs(squeeze(freqresp(H2, freqs, 'Hz'))), '-', 'DisplayName', '$H_2$'); +set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); +hold off; +set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); +ylabel('Magnitude'); +set(gca, 'XTickLabel',[]); +ylim([1e-4, 20]); +yticks([1e-4, 1e-3, 1e-2, 1e-1, 1, 1e1]); +legend('location', 'southeast', 'FontSize', 8, 'NumColumns', 2); - % Phase - ax2 = nexttile; - hold on; - set(gca,'ColorOrderIndex',1) - plot(freqs, 180/pi*phase(squeeze(freqresp(H1, freqs, 'Hz'))), '-'); - set(gca,'ColorOrderIndex',2) - plot(freqs, 180/pi*phase(squeeze(freqresp(H2, freqs, 'Hz'))), '-'); - hold off; - xlabel('Frequency [Hz]'); ylabel('Phase [deg]'); - set(gca, 'XScale', 'log'); - yticks([-180:90:180]); +% Phase +ax2 = nexttile; +hold on; +set(gca,'ColorOrderIndex',1) +plot(freqs, 180/pi*phase(squeeze(freqresp(H1, freqs, 'Hz'))), '-'); +set(gca,'ColorOrderIndex',2) +plot(freqs, 180/pi*phase(squeeze(freqresp(H2, freqs, 'Hz'))), '-'); +hold off; +xlabel('Frequency [Hz]'); ylabel('Phase [deg]'); +set(gca, 'XScale', 'log'); +yticks([-180:90:180]); - linkaxes([ax1,ax2],'x'); - xlim([freqs(1), freqs(end)]); +linkaxes([ax1,ax2],'x'); +xlim([freqs(1), freqs(end)]); #+end_src #+begin_src matlab :tangle no :exports results :results file replace - exportFig('figs/hinf_filters_results.pdf', 'width', 'wide', 'height', 'tall'); +exportFig('figs/hinf_filters_results.pdf', 'width', 'wide', 'height', 'tall'); #+end_src #+name: fig:hinf_filters_results @@ -310,15 +319,24 @@ The obtained complementary filters are shown on figure [[fig:hinf_filters_result ** Matlab Init :noexport:ignore: #+begin_src matlab :tangle no :exports none :results silent :noweb yes :var current_dir=(file-name-directory buffer-file-name) - <> +<> #+end_src #+begin_src matlab :exports none :results silent :noweb yes - <> +<> #+end_src #+begin_src matlab - freqs = logspace(-2, 4, 1000); +freqs = logspace(-2, 4, 1000); +#+end_src + +#+begin_src matlab :tangle no +addpath('./matlab'); +addpath('./matlab/src'); +#+end_src + +#+begin_src matlab :exec no +addpath('./src'); #+end_src ** Theory @@ -348,33 +366,33 @@ And thus if we choose $H_1 = 1 - H_2 - H_3$ we have solved the problem. ** Weights First we define the weights. #+begin_src matlab - n = 2; w0 = 2*pi*1; G0 = 1/10; G1 = 1000; Gc = 1/2; - W1 = (((1/w0)*sqrt((1-(G0/Gc)^(2/n))/(1-(Gc/G1)^(2/n)))*s + (G0/Gc)^(1/n))/((1/G1)^(1/n)*(1/w0)*sqrt((1-(G0/Gc)^(2/n))/(1-(Gc/G1)^(2/n)))*s + (1/Gc)^(1/n)))^n; +n = 2; w0 = 2*pi*1; G0 = 1/10; G1 = 1000; Gc = 1/2; +W1 = (((1/w0)*sqrt((1-(G0/Gc)^(2/n))/(1-(Gc/G1)^(2/n)))*s + (G0/Gc)^(1/n))/((1/G1)^(1/n)*(1/w0)*sqrt((1-(G0/Gc)^(2/n))/(1-(Gc/G1)^(2/n)))*s + (1/Gc)^(1/n)))^n; - W2 = 0.22*(1 + s/2/pi/1)^2/(sqrt(1e-4) + s/2/pi/1)^2*(1 + s/2/pi/10)^2/(1 + s/2/pi/1000)^2; +W2 = 0.22*(1 + s/2/pi/1)^2/(sqrt(1e-4) + s/2/pi/1)^2*(1 + s/2/pi/10)^2/(1 + s/2/pi/1000)^2; - n = 3; w0 = 2*pi*10; G0 = 1000; G1 = 0.1; Gc = 1/2; - W3 = (((1/w0)*sqrt((1-(G0/Gc)^(2/n))/(1-(Gc/G1)^(2/n)))*s + (G0/Gc)^(1/n))/((1/G1)^(1/n)*(1/w0)*sqrt((1-(G0/Gc)^(2/n))/(1-(Gc/G1)^(2/n)))*s + (1/Gc)^(1/n)))^n; +n = 3; w0 = 2*pi*10; G0 = 1000; G1 = 0.1; Gc = 1/2; +W3 = (((1/w0)*sqrt((1-(G0/Gc)^(2/n))/(1-(Gc/G1)^(2/n)))*s + (G0/Gc)^(1/n))/((1/G1)^(1/n)*(1/w0)*sqrt((1-(G0/Gc)^(2/n))/(1-(Gc/G1)^(2/n)))*s + (1/Gc)^(1/n)))^n; #+end_src #+begin_src matlab :exports none - figure; - hold on; - set(gca,'ColorOrderIndex',1) - plot(freqs, 1./abs(squeeze(freqresp(W1, freqs, 'Hz'))), '--', 'DisplayName', '$|W_1|^{-1}$'); - set(gca,'ColorOrderIndex',2) - plot(freqs, 1./abs(squeeze(freqresp(W2, freqs, 'Hz'))), '--', 'DisplayName', '$|W_2|^{-1}$'); - set(gca,'ColorOrderIndex',3) - plot(freqs, 1./abs(squeeze(freqresp(W3, freqs, 'Hz'))), '--', 'DisplayName', '$|W_3|^{-1}$'); - set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); - xlabel('Frequency [Hz]'); ylabel('Magnitude'); - hold off; - xlim([freqs(1), freqs(end)]); - legend('location', 'northeast', 'FontSize', 8); +figure; +hold on; +set(gca,'ColorOrderIndex',1) +plot(freqs, 1./abs(squeeze(freqresp(W1, freqs, 'Hz'))), '--', 'DisplayName', '$|W_1|^{-1}$'); +set(gca,'ColorOrderIndex',2) +plot(freqs, 1./abs(squeeze(freqresp(W2, freqs, 'Hz'))), '--', 'DisplayName', '$|W_2|^{-1}$'); +set(gca,'ColorOrderIndex',3) +plot(freqs, 1./abs(squeeze(freqresp(W3, freqs, 'Hz'))), '--', 'DisplayName', '$|W_3|^{-1}$'); +set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); +xlabel('Frequency [Hz]'); ylabel('Magnitude'); +hold off; +xlim([freqs(1), freqs(end)]); +legend('location', 'northeast', 'FontSize', 8); #+end_src #+begin_src matlab :tangle no :exports results :results file replace - exportFig('figs/three_weighting_functions.pdf', 'width', 'wide', 'height', 'normal'); +exportFig('figs/three_weighting_functions.pdf', 'width', 'wide', 'height', 'normal'); #+end_src #+name: fig:three_weighting_functions @@ -385,15 +403,15 @@ First we define the weights. ** H-Infinity Synthesis Then we create the generalized plant =P=. #+begin_src matlab - P = [W1 -W1 -W1; - 0 W2 0 ; - 0 0 W3; - 1 0 0]; +P = [W1 -W1 -W1; + 0 W2 0 ; + 0 0 W3; + 1 0 0]; #+end_src And we do the $\mathcal{H}_\infty$ synthesis. #+begin_src matlab :results output replace :exports both - [H, ~, gamma, ~] = hinfsyn(P, 1, 2,'TOLGAM', 0.001, 'METHOD', 'ric', 'DISPLAY', 'on'); +[H, ~, gamma, ~] = hinfsyn(P, 1, 2,'TOLGAM', 0.001, 'METHOD', 'ric', 'DISPLAY', 'on'); #+end_src #+RESULTS: @@ -433,58 +451,58 @@ Test bounds: 0.1000 < gamma <= 1050.0000 ** Obtained Complementary Filters The obtained filters are: #+begin_src matlab - H2 = tf(H(1)); - H3 = tf(H(2)); - H1 = 1 - H2 - H3; +H2 = tf(H(1)); +H3 = tf(H(2)); +H1 = 1 - H2 - H3; #+end_src #+begin_src matlab :exports none - figure; - tiledlayout(3, 1, 'TileSpacing', 'None', 'Padding', 'None'); +figure; +tiledlayout(3, 1, 'TileSpacing', 'None', 'Padding', 'None'); - % Magnitude - ax1 = nexttile([2, 1]); - hold on; - set(gca,'ColorOrderIndex',1) - plot(freqs, 1./abs(squeeze(freqresp(W1, freqs, 'Hz'))), '--', 'DisplayName', '$|W_1|^{-1}$'); - set(gca,'ColorOrderIndex',2) - plot(freqs, 1./abs(squeeze(freqresp(W2, freqs, 'Hz'))), '--', 'DisplayName', '$|W_2|^{-1}$'); - set(gca,'ColorOrderIndex',3) - plot(freqs, 1./abs(squeeze(freqresp(W3, freqs, 'Hz'))), '--', 'DisplayName', '$|W_3|^{-1}$'); - set(gca,'ColorOrderIndex',1) - plot(freqs, abs(squeeze(freqresp(H1, freqs, 'Hz'))), '-', 'DisplayName', '$H_1$'); - set(gca,'ColorOrderIndex',2) - plot(freqs, abs(squeeze(freqresp(H2, freqs, 'Hz'))), '-', 'DisplayName', '$H_2$'); - set(gca,'ColorOrderIndex',3) - plot(freqs, abs(squeeze(freqresp(H3, freqs, 'Hz'))), '-', 'DisplayName', '$H_3$'); - set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); - hold off; - set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); - ylabel('Magnitude'); - set(gca, 'XTickLabel',[]); - ylim([1e-4, 20]); - legend('location', 'northeast', 'FontSize', 8, 'NumColumns', 2); +% Magnitude +ax1 = nexttile([2, 1]); +hold on; +set(gca,'ColorOrderIndex',1) +plot(freqs, 1./abs(squeeze(freqresp(W1, freqs, 'Hz'))), '--', 'DisplayName', '$|W_1|^{-1}$'); +set(gca,'ColorOrderIndex',2) +plot(freqs, 1./abs(squeeze(freqresp(W2, freqs, 'Hz'))), '--', 'DisplayName', '$|W_2|^{-1}$'); +set(gca,'ColorOrderIndex',3) +plot(freqs, 1./abs(squeeze(freqresp(W3, freqs, 'Hz'))), '--', 'DisplayName', '$|W_3|^{-1}$'); +set(gca,'ColorOrderIndex',1) +plot(freqs, abs(squeeze(freqresp(H1, freqs, 'Hz'))), '-', 'DisplayName', '$H_1$'); +set(gca,'ColorOrderIndex',2) +plot(freqs, abs(squeeze(freqresp(H2, freqs, 'Hz'))), '-', 'DisplayName', '$H_2$'); +set(gca,'ColorOrderIndex',3) +plot(freqs, abs(squeeze(freqresp(H3, freqs, 'Hz'))), '-', 'DisplayName', '$H_3$'); +set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); +hold off; +set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); +ylabel('Magnitude'); +set(gca, 'XTickLabel',[]); +ylim([1e-4, 20]); +legend('location', 'northeast', 'FontSize', 8, 'NumColumns', 2); - % Phase - ax2 = nexttile; - hold on; - set(gca,'ColorOrderIndex',1) - plot(freqs, 180/pi*phase(squeeze(freqresp(H1, freqs, 'Hz')))); - set(gca,'ColorOrderIndex',2) - plot(freqs, 180/pi*phase(squeeze(freqresp(H2, freqs, 'Hz')))); - set(gca,'ColorOrderIndex',3) - plot(freqs, 180/pi*phase(squeeze(freqresp(H3, freqs, 'Hz')))); - hold off; - xlabel('Frequency [Hz]'); ylabel('Phase [deg]'); - set(gca, 'XScale', 'log'); - yticks([-360:90:360]); ylim([-270, 270]); +% Phase +ax2 = nexttile; +hold on; +set(gca,'ColorOrderIndex',1) +plot(freqs, 180/pi*phase(squeeze(freqresp(H1, freqs, 'Hz')))); +set(gca,'ColorOrderIndex',2) +plot(freqs, 180/pi*phase(squeeze(freqresp(H2, freqs, 'Hz')))); +set(gca,'ColorOrderIndex',3) +plot(freqs, 180/pi*phase(squeeze(freqresp(H3, freqs, 'Hz')))); +hold off; +xlabel('Frequency [Hz]'); ylabel('Phase [deg]'); +set(gca, 'XScale', 'log'); +yticks([-360:90:360]); ylim([-270, 270]); - linkaxes([ax1,ax2],'x'); - xlim([freqs(1), freqs(end)]); +linkaxes([ax1,ax2],'x'); +xlim([freqs(1), freqs(end)]); #+end_src #+begin_src matlab :tangle no :exports results :results file replace - exportFig('figs/three_complementary_filters_results.pdf', 'width', 'wide', 'height', 'tall'); +exportFig('figs/three_complementary_filters_results.pdf', 'width', 'wide', 'height', 'tall'); #+end_src #+name: fig:three_complementary_filters_results @@ -510,15 +528,24 @@ The FIR complementary filters designed in cite:hua05_low_ligo are of order 512. ** Matlab Init :noexport:ignore: #+begin_src matlab :tangle no :exports none :results silent :noweb yes :var current_dir=(file-name-directory buffer-file-name) - <> +<> #+end_src #+begin_src matlab :exports none :results silent :noweb yes - <> +<> #+end_src #+begin_src matlab - freqs = logspace(-3, 0, 1000); +freqs = logspace(-3, 0, 1000); +#+end_src + +#+begin_src matlab :tangle no +addpath('./matlab'); +addpath('./matlab/src'); +#+end_src + +#+begin_src matlab :exec no +addpath('./src'); #+end_src ** Specifications @@ -531,26 +558,26 @@ The specifications for the filters are: The specifications are translated in upper bounds of the complementary filters are shown on figure [[fig:ligo_specifications]]. #+begin_src matlab :exports none - figure; - hold on; - set(gca,'ColorOrderIndex',1) - plot([0.0001, 0.008], [8e-3, 8e-3], ':', 'DisplayName', 'Spec. on $H_H$'); - set(gca,'ColorOrderIndex',1) - plot([0.008 0.04], [8e-3, 1], ':', 'HandleVisibility', 'off'); - set(gca,'ColorOrderIndex',1) - plot([0.04 0.1], [3, 3], ':', 'HandleVisibility', 'off'); - set(gca,'ColorOrderIndex',2) - plot([0.1, 10], [0.045, 0.045], ':', 'DisplayName', 'Spec. on $H_L$'); - set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); - xlabel('Frequency [Hz]'); ylabel('Magnitude'); - hold off; - xlim([freqs(1), freqs(end)]); - ylim([1e-4, 10]); - legend('location', 'southeast', 'FontSize', 8); +figure; +hold on; +set(gca,'ColorOrderIndex',1) +plot([0.0001, 0.008], [8e-3, 8e-3], ':', 'DisplayName', 'Spec. on $H_H$'); +set(gca,'ColorOrderIndex',1) +plot([0.008 0.04], [8e-3, 1], ':', 'HandleVisibility', 'off'); +set(gca,'ColorOrderIndex',1) +plot([0.04 0.1], [3, 3], ':', 'HandleVisibility', 'off'); +set(gca,'ColorOrderIndex',2) +plot([0.1, 10], [0.045, 0.045], ':', 'DisplayName', 'Spec. on $H_L$'); +set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); +xlabel('Frequency [Hz]'); ylabel('Magnitude'); +hold off; +xlim([freqs(1), freqs(end)]); +ylim([1e-4, 10]); +legend('location', 'southeast', 'FontSize', 8); #+end_src #+begin_src matlab :tangle no :exports results :results file replace - exportFig('figs/ligo_specifications.pdf', 'width', 'wide', 'height', 'normal'); +exportFig('figs/ligo_specifications.pdf', 'width', 'wide', 'height', 'normal'); #+end_src #+name: fig:ligo_specifications @@ -564,63 +591,63 @@ For that, we use the [[http://cvxr.com/cvx/][CVX matlab Toolbox]]. We setup the CVX toolbox and use the =SeDuMi= solver. #+begin_src matlab - cvx_startup; - cvx_solver sedumi; +cvx_startup; +cvx_solver sedumi; #+end_src We define the frequency vectors on which we will constrain the norm of the FIR filter. #+begin_src matlab - w1 = 0:4.06e-4:0.008; - w2 = 0.008:4.06e-4:0.04; - w3 = 0.04:8.12e-4:0.1; - w4 = 0.1:8.12e-4:0.83; +w1 = 0:4.06e-4:0.008; +w2 = 0.008:4.06e-4:0.04; +w3 = 0.04:8.12e-4:0.1; +w4 = 0.1:8.12e-4:0.83; #+end_src We then define the order of the FIR filter. #+begin_src matlab - n = 512; +n = 512; #+end_src #+begin_src matlab - A1 = [ones(length(w1),1), cos(kron(w1'.*(2*pi),[1:n-1]))]; - A2 = [ones(length(w2),1), cos(kron(w2'.*(2*pi),[1:n-1]))]; - A3 = [ones(length(w3),1), cos(kron(w3'.*(2*pi),[1:n-1]))]; - A4 = [ones(length(w4),1), cos(kron(w4'.*(2*pi),[1:n-1]))]; +A1 = [ones(length(w1),1), cos(kron(w1'.*(2*pi),[1:n-1]))]; +A2 = [ones(length(w2),1), cos(kron(w2'.*(2*pi),[1:n-1]))]; +A3 = [ones(length(w3),1), cos(kron(w3'.*(2*pi),[1:n-1]))]; +A4 = [ones(length(w4),1), cos(kron(w4'.*(2*pi),[1:n-1]))]; - B1 = [zeros(length(w1),1), sin(kron(w1'.*(2*pi),[1:n-1]))]; - B2 = [zeros(length(w2),1), sin(kron(w2'.*(2*pi),[1:n-1]))]; - B3 = [zeros(length(w3),1), sin(kron(w3'.*(2*pi),[1:n-1]))]; - B4 = [zeros(length(w4),1), sin(kron(w4'.*(2*pi),[1:n-1]))]; +B1 = [zeros(length(w1),1), sin(kron(w1'.*(2*pi),[1:n-1]))]; +B2 = [zeros(length(w2),1), sin(kron(w2'.*(2*pi),[1:n-1]))]; +B3 = [zeros(length(w3),1), sin(kron(w3'.*(2*pi),[1:n-1]))]; +B4 = [zeros(length(w4),1), sin(kron(w4'.*(2*pi),[1:n-1]))]; #+end_src We run the convex optimization. #+begin_src matlab :results output replace :wrap example - cvx_begin +cvx_begin - variable y(n+1,1) +variable y(n+1,1) - % t - maximize(-y(1)) +% t +maximize(-y(1)) - for i = 1:length(w1) - norm([0 A1(i,:); 0 B1(i,:)]*y) <= 8e-3; - end +for i = 1:length(w1) + norm([0 A1(i,:); 0 B1(i,:)]*y) <= 8e-3; +end - for i = 1:length(w2) - norm([0 A2(i,:); 0 B2(i,:)]*y) <= 8e-3*(2*pi*w2(i)/(0.008*2*pi))^3; - end +for i = 1:length(w2) + norm([0 A2(i,:); 0 B2(i,:)]*y) <= 8e-3*(2*pi*w2(i)/(0.008*2*pi))^3; +end - for i = 1:length(w3) - norm([0 A3(i,:); 0 B3(i,:)]*y) <= 3; - end +for i = 1:length(w3) + norm([0 A3(i,:); 0 B3(i,:)]*y) <= 3; +end - for i = 1:length(w4) - norm([[1 0]'- [0 A4(i,:); 0 B4(i,:)]*y]) <= y(1); - end +for i = 1:length(w4) + norm([[1 0]'- [0 A4(i,:); 0 B4(i,:)]*y]) <= y(1); +end - cvx_end +cvx_end - h = y(2:end); +h = y(2:end); #+end_src #+RESULTS: @@ -703,41 +730,41 @@ h = y(2:end); Finally, we compute the filter response over the frequency vector defined and the result is shown on figure [[fig:fir_filter_ligo]] which is very close to the filters obtain in cite:hua05_low_ligo. #+begin_src matlab - w = [w1 w2 w3 w4]; - H = [exp(-j*kron(w'.*2*pi,[0:n-1]))]*h; +w = [w1 w2 w3 w4]; +H = [exp(-j*kron(w'.*2*pi,[0:n-1]))]*h; #+end_src #+begin_src matlab :exports none - figure; - tiledlayout(3, 1, 'TileSpacing', 'None', 'Padding', 'None'); +figure; +tiledlayout(3, 1, 'TileSpacing', 'None', 'Padding', 'None'); - % Magnitude - ax1 = nexttile([2, 1]); - hold on; - plot(w, abs(H), 'k-'); - plot(w, abs(1-H), 'k--'); - hold off; - set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); - ylabel('Magnitude'); - set(gca, 'XTickLabel',[]); - ylim([5e-3, 5]); +% Magnitude +ax1 = nexttile([2, 1]); +hold on; +plot(w, abs(H), 'k-'); +plot(w, abs(1-H), 'k--'); +hold off; +set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); +ylabel('Magnitude'); +set(gca, 'XTickLabel',[]); +ylim([5e-3, 5]); - % Phase - ax2 = nexttile; - hold on; - plot(w, 180/pi*angle(H), 'k-'); - plot(w, 180/pi*angle(1-H), 'k--'); - hold off; - xlabel('Frequency [Hz]'); ylabel('Phase [deg]'); - set(gca, 'XScale', 'log'); - yticks([-180:90:180]); ylim([-180, 180]); +% Phase +ax2 = nexttile; +hold on; +plot(w, 180/pi*angle(H), 'k-'); +plot(w, 180/pi*angle(1-H), 'k--'); +hold off; +xlabel('Frequency [Hz]'); ylabel('Phase [deg]'); +set(gca, 'XScale', 'log'); +yticks([-180:90:180]); ylim([-180, 180]); - linkaxes([ax1,ax2],'x'); - xlim([1e-3, 1]); +linkaxes([ax1,ax2],'x'); +xlim([1e-3, 1]); #+end_src #+begin_src matlab :tangle no :exports results :results file replace - exportFig('figs/fir_filter_ligo.pdf', 'width', 'wide', 'height', 'normal'); +exportFig('figs/fir_filter_ligo.pdf', 'width', 'wide', 'height', 'normal'); #+end_src #+name: fig:fir_filter_ligo @@ -756,56 +783,56 @@ Here are the requirements on the filters: The bode plot of the weights is shown on figure [[fig:ligo_weights]]. #+begin_src matlab :exports none - w1 = 2*pi*0.008; x1 = 0.35; - w2 = 2*pi*0.04; x2 = 0.5; - w3 = 2*pi*0.05; x3 = 0.5; +w1 = 2*pi*0.008; x1 = 0.35; +w2 = 2*pi*0.04; x2 = 0.5; +w3 = 2*pi*0.05; x3 = 0.5; - % Slope of +3 from w1 - wH = 0.008*(s^2/w1^2 + 2*x1/w1*s + 1)*(s/w1 + 1); - % Little bump from w2 to w3 - wH = wH*(s^2/w2^2 + 2*x2/w2*s + 1)/(s^2/w3^2 + 2*x3/w3*s + 1); - % No Slope at high frequencies - wH = wH/(s^2/w3^2 + 2*x3/w3*s + 1)/(s/w3 + 1); - % Little bump between w2 and w3 - w0 = 2*pi*0.045; xi = 0.1; A = 2; n = 1; - wH = wH*((s^2 + 2*w0*xi*A^(1/n)*s + w0^2)/(s^2 + 2*w0*xi*s + w0^2))^n; +% Slope of +3 from w1 +wH = 0.008*(s^2/w1^2 + 2*x1/w1*s + 1)*(s/w1 + 1); +% Little bump from w2 to w3 +wH = wH*(s^2/w2^2 + 2*x2/w2*s + 1)/(s^2/w3^2 + 2*x3/w3*s + 1); +% No Slope at high frequencies +wH = wH/(s^2/w3^2 + 2*x3/w3*s + 1)/(s/w3 + 1); +% Little bump between w2 and w3 +w0 = 2*pi*0.045; xi = 0.1; A = 2; n = 1; +wH = wH*((s^2 + 2*w0*xi*A^(1/n)*s + w0^2)/(s^2 + 2*w0*xi*s + w0^2))^n; - wH = 1/wH; - wH = minreal(ss(wH)); +wH = 1/wH; +wH = minreal(ss(wH)); #+end_src #+begin_src matlab :exports none - n = 20; Rp = 1; Wp = 2*pi*0.102; - [b,a] = cheby1(n, Rp, Wp, 'high', 's'); - wL = 0.04*tf(a, b); +n = 20; Rp = 1; Wp = 2*pi*0.102; +[b,a] = cheby1(n, Rp, Wp, 'high', 's'); +wL = 0.04*tf(a, b); - wL = 1/wL; - wL = minreal(ss(wL)); +wL = 1/wL; +wL = minreal(ss(wL)); #+end_src #+begin_src matlab :exports none - figure; - hold on; - set(gca,'ColorOrderIndex',1); - plot(freqs, abs(squeeze(freqresp(inv(wH), freqs, 'Hz'))), '-', 'DisplayName', '$|w_H|^{-1}$'); - set(gca,'ColorOrderIndex',2); - plot(freqs, abs(squeeze(freqresp(inv(wL), freqs, 'Hz'))), '-', 'DisplayName', '$|w_L|^{-1}$'); +figure; +hold on; +set(gca,'ColorOrderIndex',1); +plot(freqs, abs(squeeze(freqresp(inv(wH), freqs, 'Hz'))), '-', 'DisplayName', '$|w_H|^{-1}$'); +set(gca,'ColorOrderIndex',2); +plot(freqs, abs(squeeze(freqresp(inv(wL), freqs, 'Hz'))), '-', 'DisplayName', '$|w_L|^{-1}$'); - plot([0.0001, 0.008], [8e-3, 8e-3], 'k--', 'DisplayName', 'Spec.'); - plot([0.008 0.04], [8e-3, 1], 'k--', 'HandleVisibility', 'off'); - plot([0.04 0.1], [3, 3], 'k--', 'HandleVisibility', 'off'); - plot([0.1, 10], [0.045, 0.045], 'k--', 'HandleVisibility', 'off'); +plot([0.0001, 0.008], [8e-3, 8e-3], 'k--', 'DisplayName', 'Spec.'); +plot([0.008 0.04], [8e-3, 1], 'k--', 'HandleVisibility', 'off'); +plot([0.04 0.1], [3, 3], 'k--', 'HandleVisibility', 'off'); +plot([0.1, 10], [0.045, 0.045], 'k--', 'HandleVisibility', 'off'); - set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); - xlabel('Frequency [Hz]'); ylabel('Magnitude'); - hold off; - xlim([freqs(1), freqs(end)]); - ylim([1e-3, 10]); - legend('location', 'southeast', 'FontSize', 8); +set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); +xlabel('Frequency [Hz]'); ylabel('Magnitude'); +hold off; +xlim([freqs(1), freqs(end)]); +ylim([1e-3, 10]); +legend('location', 'southeast', 'FontSize', 8); #+end_src #+begin_src matlab :tangle no :exports results :results file replace - exportFig('figs/ligo_weights.pdf', 'width', 'wide', 'height', 'normal'); +exportFig('figs/ligo_weights.pdf', 'width', 'wide', 'height', 'normal'); #+end_src #+name: fig:ligo_weights @@ -816,14 +843,14 @@ The bode plot of the weights is shown on figure [[fig:ligo_weights]]. ** H-Infinity Synthesis We define the generalized plant as shown on figure [[fig:h_infinity_robst_fusion]]. #+begin_src matlab - P = [0 wL; - wH -wH; - 1 0]; +P = [0 wL; + wH -wH; + 1 0]; #+end_src And we do the $\mathcal{H}_\infty$ synthesis using the =hinfsyn= command. #+begin_src matlab :results output replace :exports both :wrap example - [Hl, ~, gamma, ~] = hinfsyn(P, 1, 1,'TOLGAM', 0.001, 'METHOD', 'ric', 'DISPLAY', 'on'); +[Hl, ~, gamma, ~] = hinfsyn(P, 1, 1,'TOLGAM', 0.001, 'METHOD', 'ric', 'DISPLAY', 'on'); #+end_src #+RESULTS: @@ -852,18 +879,18 @@ Test bounds: 0.3276 < gamma <= 1.8063 The high pass filter is defined as $H_H = 1 - H_L$. #+begin_src matlab - Hh = 1 - Hl; +Hh = 1 - Hl; #+end_src #+begin_src matlab :exports none - Hh = minreal(Hh); - Hl = minreal(Hl); +Hh = minreal(Hh); +Hl = minreal(Hl); #+end_src The size of the filters is shown below. #+begin_src matlab :exports results :results output replace :wrap example - size(Hh), size(Hl) +size(Hh), size(Hl) #+end_src #+RESULTS: @@ -876,33 +903,33 @@ State-space model with 1 outputs, 1 inputs, and 27 states. The bode plot of the obtained filters as shown on figure [[fig:hinf_synthesis_ligo_results]]. #+begin_src matlab :exports none - figure; - hold on; - set(gca,'ColorOrderIndex',1); - plot([0.0001, 0.008], [8e-3, 8e-3], ':', 'DisplayName', 'Spec. on $H_H$'); - set(gca,'ColorOrderIndex',1); - plot([0.008 0.04], [8e-3, 1], ':', 'HandleVisibility', 'off'); - set(gca,'ColorOrderIndex',1); - plot([0.04 0.1], [3, 3], ':', 'HandleVisibility', 'off'); +figure; +hold on; +set(gca,'ColorOrderIndex',1); +plot([0.0001, 0.008], [8e-3, 8e-3], ':', 'DisplayName', 'Spec. on $H_H$'); +set(gca,'ColorOrderIndex',1); +plot([0.008 0.04], [8e-3, 1], ':', 'HandleVisibility', 'off'); +set(gca,'ColorOrderIndex',1); +plot([0.04 0.1], [3, 3], ':', 'HandleVisibility', 'off'); - set(gca,'ColorOrderIndex',2); - plot([0.1, 10], [0.045, 0.045], ':', 'DisplayName', 'Spec. on $H_L$'); +set(gca,'ColorOrderIndex',2); +plot([0.1, 10], [0.045, 0.045], ':', 'DisplayName', 'Spec. on $H_L$'); - set(gca,'ColorOrderIndex',1); - plot(freqs, abs(squeeze(freqresp(Hh, freqs, 'Hz'))), '-', 'DisplayName', '$H_H$'); - set(gca,'ColorOrderIndex',2); - plot(freqs, abs(squeeze(freqresp(Hl, freqs, 'Hz'))), '-', 'DisplayName', '$H_L$'); +set(gca,'ColorOrderIndex',1); +plot(freqs, abs(squeeze(freqresp(Hh, freqs, 'Hz'))), '-', 'DisplayName', '$H_H$'); +set(gca,'ColorOrderIndex',2); +plot(freqs, abs(squeeze(freqresp(Hl, freqs, 'Hz'))), '-', 'DisplayName', '$H_L$'); - set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); - xlabel('Frequency [Hz]'); ylabel('Magnitude'); - hold off; - xlim([freqs(1), freqs(end)]); - ylim([1e-3, 10]); - legend('location', 'southeast', 'FontSize', 8); +set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); +xlabel('Frequency [Hz]'); ylabel('Magnitude'); +hold off; +xlim([freqs(1), freqs(end)]); +ylim([1e-3, 10]); +legend('location', 'southeast', 'FontSize', 8); #+end_src #+begin_src matlab :tangle no :exports results :results file replace - exportFig('figs/hinf_synthesis_ligo_results.pdf', 'width', 'wide', 'height', 'normal'); +exportFig('figs/hinf_synthesis_ligo_results.pdf', 'width', 'wide', 'height', 'normal'); #+end_src #+name: fig:hinf_synthesis_ligo_results @@ -914,56 +941,56 @@ The bode plot of the obtained filters as shown on figure [[fig:hinf_synthesis_li Let's now compare the FIR filters designed in cite:hua05_low_ligo and the one obtained with the $\mathcal{H}_\infty$ synthesis on figure [[fig:comp_fir_ligo_hinf]]. #+begin_src matlab :exports none - figure; - tiledlayout(3, 1, 'TileSpacing', 'None', 'Padding', 'None'); +figure; +tiledlayout(3, 1, 'TileSpacing', 'None', 'Padding', 'None'); - % Magnitude - ax1 = nexttile([2, 1]); - hold on; - set(gca,'ColorOrderIndex',1); - plot(freqs, abs(squeeze(freqresp(Hh, freqs, 'Hz'))), '-', ... - 'DisplayName', '$H_H(s)$ - $\mathcal{H}_\infty$'); - set(gca,'ColorOrderIndex',2); - plot(freqs, abs(squeeze(freqresp(Hl, freqs, 'Hz'))), '-', ... - 'DisplayName', '$H_L(s)$ - $\mathcal{H}_\infty$'); +% Magnitude +ax1 = nexttile([2, 1]); +hold on; +set(gca,'ColorOrderIndex',1); +plot(freqs, abs(squeeze(freqresp(Hh, freqs, 'Hz'))), '-', ... + 'DisplayName', '$H_H(s)$ - $\mathcal{H}_\infty$'); +set(gca,'ColorOrderIndex',2); +plot(freqs, abs(squeeze(freqresp(Hl, freqs, 'Hz'))), '-', ... + 'DisplayName', '$H_L(s)$ - $\mathcal{H}_\infty$'); - set(gca,'ColorOrderIndex',1); - plot(w, abs(H), '--', ... - 'DisplayName', '$H_H(s)$ - FIR'); - set(gca,'ColorOrderIndex',2); - plot(w, abs(1-H), '--', ... - 'DisplayName', '$H_L(s)$ - FIR'); - hold off; - set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); - ylabel('Magnitude'); - set(gca, 'XTickLabel',[]); - ylim([5e-3, 10]); - legend('location', 'southeast', 'FontSize', 8, 'NumColumns', 2); +set(gca,'ColorOrderIndex',1); +plot(w, abs(H), '--', ... + 'DisplayName', '$H_H(s)$ - FIR'); +set(gca,'ColorOrderIndex',2); +plot(w, abs(1-H), '--', ... + 'DisplayName', '$H_L(s)$ - FIR'); +hold off; +set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); +ylabel('Magnitude'); +set(gca, 'XTickLabel',[]); +ylim([5e-3, 10]); +legend('location', 'southeast', 'FontSize', 8, 'NumColumns', 2); - % Phase - ax2 = nexttile; - hold on; - set(gca,'ColorOrderIndex',1); - plot(freqs, 180/pi*angle(squeeze(freqresp(Hh, freqs, 'Hz'))), '-'); - set(gca,'ColorOrderIndex',2); - plot(freqs, 180/pi*angle(squeeze(freqresp(Hl, freqs, 'Hz'))), '-'); +% Phase +ax2 = nexttile; +hold on; +set(gca,'ColorOrderIndex',1); +plot(freqs, 180/pi*angle(squeeze(freqresp(Hh, freqs, 'Hz'))), '-'); +set(gca,'ColorOrderIndex',2); +plot(freqs, 180/pi*angle(squeeze(freqresp(Hl, freqs, 'Hz'))), '-'); - set(gca,'ColorOrderIndex',1); - plot(w, 180/pi*angle(H), '--'); - set(gca,'ColorOrderIndex',2); - plot(w, 180/pi*angle(1-H), '--'); - set(gca, 'XScale', 'log'); - xlabel('Frequency [Hz]'); ylabel('Phase [deg]'); - hold off; - yticks([-180:90:180]); ylim([-180, 180]); +set(gca,'ColorOrderIndex',1); +plot(w, 180/pi*angle(H), '--'); +set(gca,' H2ColorOrderIndex',2); +plot(w, 180/pi*angle(1-H), '--'); +set(gca, 'XScale', 'log'); +xlabel('Frequency [Hz]'); ylabel('Phase [deg]'); +hold off; +yticks([-180:90:180]); ylim([-180, 180]); - linkaxes([ax1,ax2],'x'); - xlim([freqs(1), freqs(end)]); +linkaxes([ax1,ax2],'x'); +xlim([freqs(1), freqs(end)]); #+end_src #+begin_src matlab :tangle no :exports results :results file replace - exportFig('figs/comp_fir_ligo_hinf.pdf', 'width', 'wide', 'height', 'tall'); +exportFig('figs/comp_fir_ligo_hinf.pdf', 'width', 'wide', 'height', 'tall'); #+end_src #+name: fig:comp_fir_ligo_hinf @@ -972,13 +999,23 @@ Let's now compare the FIR filters designed in cite:hua05_low_ligo and the one ob [[file:figs/comp_fir_ligo_hinf.png]] * Alternative Synthesis +** Introduction :ignore: ** Matlab Init :noexport:ignore: #+begin_src matlab :tangle no :exports none :results silent :noweb yes :var current_dir=(file-name-directory buffer-file-name) - <> +<> #+end_src #+begin_src matlab :exports none :results silent :noweb yes - <> +<> +#+end_src + +#+begin_src matlab :tangle no +addpath('./matlab'); +addpath('./matlab/src'); +#+end_src + +#+begin_src matlab :exec no +addpath('./src'); #+end_src ** Two generalized plants @@ -1027,33 +1064,33 @@ In order to synthesize the complementary filter using the proposed method, we ca \end{equation} #+begin_src latex :file h_infinity_arch_2.pdf - \begin{tikzpicture} - \node[block={4.5cm}{4.5cm}, fill=black!20!white, dashed] (P) {}; - \node[above] at (P.north) {$P_2(s)$}; +\begin{tikzpicture} + \node[block={4.5cm}{4.5cm}, fill=black!20!white, dashed] (P) {}; + \node[above] at (P.north) {$P_2(s)$}; - \coordinate[] (input2) at ($(P.south west)!0.85!(P.north west) + (-0.7, 0)$); - \coordinate[] (input1) at ($(P.south west)!0.55!(P.north west) + (-0.7, 0)$); - \coordinate[] (inputu) at ($(P.south west)!0.3!( P.north west) + (-0.7, 0)$); + \coordinate[] (input2) at ($(P.south west)!0.85!(P.north west) + (-0.7, 0)$); + \coordinate[] (input1) at ($(P.south west)!0.55!(P.north west) + (-0.7, 0)$); + \coordinate[] (inputu) at ($(P.south west)!0.3!( P.north west) + (-0.7, 0)$); - \coordinate[] (outputz) at ($(P.south east)!0.3!(P.north east) + (0.7, 0)$); - \coordinate[] (outputv) at ($(P.south east)!0.1!(P.north east) + (0.7, 0)$); + \coordinate[] (outputz) at ($(P.south east)!0.3!(P.north east) + (0.7, 0)$); + \coordinate[] (outputv) at ($(P.south east)!0.1!(P.north east) + (0.7, 0)$); - \node[block, right=1.4 of input2] (W2){$W_2(s)$}; - \node[block, right=1.4 of input1] (W1){$W_1(s)$}; - \node[addb={+}{-}{}{}{}, right=of W1] (sub) {}; - \node[addb, left=2.5 of outputz] (add) {}; + \node[block, right=1.4 of input2] (W2){$W_2(s)$}; + \node[block, right=1.4 of input1] (W1){$W_1(s)$}; + \node[addb={+}{-}{}{}{}, right=of W1] (sub) {}; + \node[addb, left=2.5 of outputz] (add) {}; - \node[block, below=0.3 of P] (H2) {$H_2(s)$}; + \node[block, below=0.3 of P] (H2) {$H_2(s)$}; - \draw[->] (input2) node[above right]{$w_2$} -- (W2.west); - \draw[->] (input1) node[above right]{$w_1$} -- (W1.west); - \draw[->] (W2.east) -| (sub.north); - \draw[->] (W1.east) -- (sub.west); - \draw[->] (W1-|add)node[branch]{} -- (add.north); - \draw[->] (sub.south) |- (outputv) node[above left]{$v$} |- (H2.east); - \draw[->] (H2.west) -| (inputu) node[above right]{$u$} -- (add.west); - \draw[->] (add.east) -- (outputz) node[above left]{$z$}; - \end{tikzpicture} + \draw[->] (input2) node[above right]{$w_2$} -- (W2.west); + \draw[->] (input1) node[above right]{$w_1$} -- (W1.west); + \draw[->] (W2.east) -| (sub.north); + \draw[->] (W1.east) -- (sub.west); + \draw[->] (W1-|add)node[branch]{} -- (add.north); + \draw[->] (sub.south) |- (outputv) node[above left]{$v$} |- (H2.east); + \draw[->] (H2.west) -| (inputu) node[above right]{$u$} -- (add.west); + \draw[->] (add.east) -- (outputz) node[above left]{$z$}; +\end{tikzpicture} #+end_src #+name: fig:h_infinity_arch_2 @@ -1063,21 +1100,21 @@ In order to synthesize the complementary filter using the proposed method, we ca Let's run the $\mathcal{H}_\infty$ synthesis for both generalized plant using the same weights and see if the obtained filters are the same: #+begin_src matlab - n = 2; w0 = 2*pi*11; G0 = 1/10; G1 = 1000; Gc = 1/2; - W1 = (((1/w0)*sqrt((1-(G0/Gc)^(2/n))/(1-(Gc/G1)^(2/n)))*s + (G0/Gc)^(1/n))/((1/G1)^(1/n)*(1/w0)*sqrt((1-(G0/Gc)^(2/n))/(1-(Gc/G1)^(2/n)))*s + (1/Gc)^(1/n)))^n; +n = 2; w0 = 2*pi*11; G0 = 1/10; G1 = 1000; Gc = 1/2; +W1 = (((1/w0)*sqrt((1-(G0/Gc)^(2/n))/(1-(Gc/G1)^(2/n)))*s + (G0/Gc)^(1/n))/((1/G1)^(1/n)*(1/w0)*sqrt((1-(G0/Gc)^(2/n))/(1-(Gc/G1)^(2/n)))*s + (1/Gc)^(1/n)))^n; - n = 3; w0 = 2*pi*10; G0 = 1000; G1 = 0.1; Gc = 1/2; - W2 = (((1/w0)*sqrt((1-(G0/Gc)^(2/n))/(1-(Gc/G1)^(2/n)))*s + (G0/Gc)^(1/n))/((1/G1)^(1/n)*(1/w0)*sqrt((1-(G0/Gc)^(2/n))/(1-(Gc/G1)^(2/n)))*s + (1/Gc)^(1/n)))^n; +n = 3; w0 = 2*pi*10; G0 = 1000; G1 = 0.1; Gc = 1/2; +W2 = (((1/w0)*sqrt((1-(G0/Gc)^(2/n))/(1-(Gc/G1)^(2/n)))*s + (G0/Gc)^(1/n))/((1/G1)^(1/n)*(1/w0)*sqrt((1-(G0/Gc)^(2/n))/(1-(Gc/G1)^(2/n)))*s + (1/Gc)^(1/n)))^n; #+end_src #+begin_src matlab - P1 = [W1 -W1; - 0 W2; - 1 0]; +P1 = [W1 -W1; + 0 W2; + 1 0]; #+end_src #+begin_src matlab :results output replace :exports both - [H2, ~, gamma, ~] = hinfsyn(P1, 1, 1,'TOLGAM', 0.001, 'METHOD', 'ric', 'DISPLAY', 'on'); +[H2, ~, gamma, ~] = hinfsyn(P1, 1, 1,'TOLGAM', 0.001, 'METHOD', 'ric', 'DISPLAY', 'on'); #+end_src #+RESULTS: @@ -1107,12 +1144,12 @@ Let's run the $\mathcal{H}_\infty$ synthesis for both generalized plant using th #+end_example #+begin_src matlab - P2 = [0 W1 1; - W2 -W1 0]; +P2 = [0 W1 1; + W2 -W1 0]; #+end_src #+begin_src matlab :results output replace :exports both - [H2b, ~, gamma, ~] = hinfsyn(P2, 1, 1,'TOLGAM', 0.001, 'METHOD', 'ric', 'DISPLAY', 'on'); +[H2b, ~, gamma, ~] = hinfsyn(P2, 1, 1,'TOLGAM', 0.001, 'METHOD', 'ric', 'DISPLAY', 'on'); #+end_src #+RESULTS: @@ -1143,20 +1180,20 @@ Let's run the $\mathcal{H}_\infty$ synthesis for both generalized plant using th And indeed, we can see that the exact same filters are obtained (Figure [[fig:hinf_comp_P1_P2_syn]]). #+begin_src matlab :exports none - freqs = logspace(-2, 4, 1000); +freqs = logspace(-2, 4, 1000); - figure; - hold on; - plot(freqs, abs(squeeze(freqresp(H2, freqs, 'Hz'))), '-', 'DisplayName', '$H_2$'); - plot(freqs, abs(squeeze(freqresp(H2b, freqs, 'Hz'))), '--', 'DisplayName', '$H_{2b}$'); - hold off; - set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); - xlabel('Frequency [Hz]'); ylabel('Magnitude'); - legend('location', 'southeast', 'FontSize', 8); +figure; +hold on; +plot(freqs, abs(squeeze(freqresp(H2, freqs, 'Hz'))), '-', 'DisplayName', '$H_2$'); +plot(freqs, abs(squeeze(freqresp(H2b, freqs, 'Hz'))), '--', 'DisplayName', '$H_{2b}$'); +hold off; +set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); +xlabel('Frequency [Hz]'); ylabel('Magnitude'); +legend('location', 'southeast', 'FontSize', 8); #+end_src #+begin_src matlab :tangle no :exports results :results file replace - exportFig('figs/hinf_comp_P1_P2_syn.pdf', 'width', 'wide', 'height', 'normal'); +exportFig('figs/hinf_comp_P1_P2_syn.pdf', 'width', 'wide', 'height', 'normal'); #+end_src #+name: fig:hinf_comp_P1_P2_syn @@ -1168,22 +1205,22 @@ And indeed, we can see that the exact same filters are obtained (Figure [[fig:hi Let's see if there is a difference by explicitly shaping $H_1(s)$ or $H_2(s)$. #+begin_src matlab - n = 2; w0 = 2*pi*11; G0 = 1/10; G1 = 1000; Gc = 1/2; - W1 = (((1/w0)*sqrt((1-(G0/Gc)^(2/n))/(1-(Gc/G1)^(2/n)))*s + (G0/Gc)^(1/n))/((1/G1)^(1/n)*(1/w0)*sqrt((1-(G0/Gc)^(2/n))/(1-(Gc/G1)^(2/n)))*s + (1/Gc)^(1/n)))^n; +n = 2; w0 = 2*pi*11; G0 = 1/10; G1 = 1000; Gc = 1/2; +W1 = (((1/w0)*sqrt((1-(G0/Gc)^(2/n))/(1-(Gc/G1)^(2/n)))*s + (G0/Gc)^(1/n))/((1/G1)^(1/n)*(1/w0)*sqrt((1-(G0/Gc)^(2/n))/(1-(Gc/G1)^(2/n)))*s + (1/Gc)^(1/n)))^n; - n = 3; w0 = 2*pi*10; G0 = 1000; G1 = 0.1; Gc = 1/2; - W2 = (((1/w0)*sqrt((1-(G0/Gc)^(2/n))/(1-(Gc/G1)^(2/n)))*s + (G0/Gc)^(1/n))/((1/G1)^(1/n)*(1/w0)*sqrt((1-(G0/Gc)^(2/n))/(1-(Gc/G1)^(2/n)))*s + (1/Gc)^(1/n)))^n; +n = 3; w0 = 2*pi*10; G0 = 1000; G1 = 0.1; Gc = 1/2; +W2 = (((1/w0)*sqrt((1-(G0/Gc)^(2/n))/(1-(Gc/G1)^(2/n)))*s + (G0/Gc)^(1/n))/((1/G1)^(1/n)*(1/w0)*sqrt((1-(G0/Gc)^(2/n))/(1-(Gc/G1)^(2/n)))*s + (1/Gc)^(1/n)))^n; #+end_src Let's first synthesize $H_1(s)$: #+begin_src matlab - P1 = [W2 -W2; - 0 W1; - 1 0]; +P1 = [W2 -W2; + 0 W1; + 1 0]; #+end_src #+begin_src matlab :results output replace :exports both - [H1, ~, gamma, ~] = hinfsyn(P1, 1, 1,'TOLGAM', 0.001, 'METHOD', 'ric', 'DISPLAY', 'on'); +[H1, ~, gamma, ~] = hinfsyn(P1, 1, 1,'TOLGAM', 0.001, 'METHOD', 'ric', 'DISPLAY', 'on'); #+end_src #+RESULTS: @@ -1213,13 +1250,13 @@ Let's first synthesize $H_1(s)$: And now $H_2(s)$: #+begin_src matlab - P2 = [W1 -W1; - 0 W2; - 1 0]; +P2 = [W1 -W1; + 0 W2; + 1 0]; #+end_src #+begin_src matlab :results output replace :exports both - [H2, ~, gamma, ~] = hinfsyn(P2, 1, 1,'TOLGAM', 0.001, 'METHOD', 'ric', 'DISPLAY', 'on'); +[H2, ~, gamma, ~] = hinfsyn(P2, 1, 1,'TOLGAM', 0.001, 'METHOD', 'ric', 'DISPLAY', 'on'); #+end_src #+RESULTS: @@ -1250,22 +1287,22 @@ And now $H_2(s)$: And compare $H_1(s)$ with $1 - H_2(s)$ and $H_2(s)$ with $1 - H_1(s)$ in Figure [[fig:hinf_comp_H1_H2_syn]]. #+begin_src matlab :exports none - freqs = logspace(-2, 4, 1000); +freqs = logspace(-2, 4, 1000); - figure; - hold on; - plot(freqs, abs(squeeze(freqresp(H1, freqs, 'Hz'))), '-', 'DisplayName', '$H_1$'); - plot(freqs, abs(squeeze(freqresp(1-H1, freqs, 'Hz'))), '-', 'DisplayName', '$H_2$'); - plot(freqs, abs(squeeze(freqresp(1-H2, freqs, 'Hz'))), '--', 'DisplayName', '$H_2$'); - plot(freqs, abs(squeeze(freqresp(H2, freqs, 'Hz'))), '--', 'DisplayName', '$H_1$'); - hold off; - set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); - xlabel('Frequency [Hz]'); ylabel('Magnitude'); - legend('location', 'northeast', 'FontSize', 8, 'NumColumns', 2); +figure; +hold on; +plot(freqs, abs(squeeze(freqresp(H1, freqs, 'Hz'))), '-', 'DisplayName', '$H_1$'); +plot(freqs, abs(squeeze(freqresp(1-H1, freqs, 'Hz'))), '-', 'DisplayName', '$H_2$'); +plot(freqs, abs(squeeze(freqresp(1-H2, freqs, 'Hz'))), '--', 'DisplayName', '$H_2$'); +plot(freqs, abs(squeeze(freqresp(H2, freqs, 'Hz'))), '--', 'DisplayName', '$H_1$'); +hold off; +set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); +xlabel('Frequency [Hz]'); ylabel('Magnitude'); +legend('location', 'northeast', 'FontSize', 8, 'NumColumns', 2); #+end_src #+begin_src matlab :tangle no :exports results :results file replace - exportFig('figs/hinf_comp_H1_H2_syn.pdf', 'width', 'wide', 'height', 'normal'); +exportFig('figs/hinf_comp_H1_H2_syn.pdf', 'width', 'wide', 'height', 'normal'); #+end_src #+name: fig:hinf_comp_H1_H2_syn @@ -1273,14 +1310,140 @@ And compare $H_1(s)$ with $1 - H_2(s)$ and $H_2(s)$ with $1 - H_1(s)$ in Figure #+RESULTS: [[file:figs/hinf_comp_H1_H2_syn.png]] +** Using Feedback architecture + +#+begin_src matlab +n = 2; w0 = 2*pi*11; G0 = 1/10; G1 = 1000; Gc = 1/2; +W1 = (((1/w0)*sqrt((1-(G0/Gc)^(2/n))/(1-(Gc/G1)^(2/n)))*s + (G0/Gc)^(1/n))/((1/G1)^(1/n)*(1/w0)*sqrt((1-(G0/Gc)^(2/n))/(1-(Gc/G1)^(2/n)))*s + (1/Gc)^(1/n)))^n; + +n = 3; w0 = 2*pi*10; G0 = 1000; G1 = 0.1; Gc = 1/2; +W2 = (((1/w0)*sqrt((1-(G0/Gc)^(2/n))/(1-(Gc/G1)^(2/n)))*s + (G0/Gc)^(1/n))/((1/G1)^(1/n)*(1/w0)*sqrt((1-(G0/Gc)^(2/n))/(1-(Gc/G1)^(2/n)))*s + (1/Gc)^(1/n)))^n; +#+end_src + +Let's first synthesize $H_1(s)$: +#+begin_src matlab +P = [W1 -W1; + 0 W2; + 1 -1]; +#+end_src + +#+begin_src matlab :results output replace :exports both +[K, ~, gamma, ~] = hinfsyn(P, 1, 1,'TOLGAM', 0.001, 'METHOD', 'lmi', 'DISPLAY', 'on'); +#+end_src + +#+begin_src matlab +H1 = inv(1 + K); +H2 = 1 - H1; +#+end_src + +#+begin_src matlab :exports none +freqs = logspace(-2, 4, 1000); + +figure; +hold on; +plot(freqs, abs(squeeze(freqresp(H1, freqs, 'Hz'))), '-', 'DisplayName', '$H_1$'); +plot(freqs, abs(squeeze(freqresp(H2, freqs, 'Hz'))), '-', 'DisplayName', '$H_2$'); +set(gca,'ColorOrderIndex',1); +plot(freqs, 1./abs(squeeze(freqresp(W1, freqs, 'Hz'))), '--', 'DisplayName', '$|W_1|^{-1}$'); +plot(freqs, 1./abs(squeeze(freqresp(W2, freqs, 'Hz'))), '--', 'DisplayName', '$|W_2|^{-1}$'); +hold off; +set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); +xlabel('Frequency [Hz]'); ylabel('Magnitude'); +legend('location', 'northeast', 'FontSize', 8, 'NumColumns', 2); +#+end_src + +** Adding feature in the filters + +#+begin_src matlab +n = 2; w0 = 2*pi*11; G0 = 1/10; G1 = 1000; Gc = 1/2; +W1 = (((1/w0)*sqrt((1-(G0/Gc)^(2/n))/(1-(Gc/G1)^(2/n)))*s + (G0/Gc)^(1/n))/((1/G1)^(1/n)*(1/w0)*sqrt((1-(G0/Gc)^(2/n))/(1-(Gc/G1)^(2/n)))*s + (1/Gc)^(1/n)))^n; + +n = 3; w0 = 2*pi*10; G0 = 1000; G1 = 0.1; Gc = 1/2; +W2 = (((1/w0)*sqrt((1-(G0/Gc)^(2/n))/(1-(Gc/G1)^(2/n)))*s + (G0/Gc)^(1/n))/((1/G1)^(1/n)*(1/w0)*sqrt((1-(G0/Gc)^(2/n))/(1-(Gc/G1)^(2/n)))*s + (1/Gc)^(1/n)))^n; +#+end_src + +#+begin_src matlab +Wf = (1 + s/2/pi/1)/s; +Wf = s/(1 + s/2/pi/1e2); + +% W2 = W2/Wf/(1 + s/2/pi/1e3); +#+end_src + +#+begin_src matlab +P = [W1 -Wf*W1; + 0 Wf*W2; + 1 -Wf]; +#+end_src + +#+begin_src matlab :results output replace :exports both +[Ka, ~, gamma, ~] = hinfsyn(P, 1, 1,'TOLGAM', 0.001, 'METHOD', 'lmi', 'DISPLAY', 'on'); +#+end_src + +#+begin_src matlab +K = Ka*Wf; +#+end_src + +#+begin_src matlab +H1 = inv(1 + K); +H2 = 1 - H1; +#+end_src + +#+begin_src matlab :exports none +freqs = logspace(-2, 4, 1000); + +figure; +tiledlayout(3, 1, 'TileSpacing', 'None', 'Padding', 'None'); + +% Magnitude +ax1 = nexttile([2, 1]); +hold on; +plot(freqs, abs(squeeze(freqresp(H1, freqs, 'Hz'))), '-', 'DisplayName', '$H_1$'); +plot(freqs, abs(squeeze(freqresp(H2, freqs, 'Hz'))), '-', 'DisplayName', '$H_2$'); +set(gca,'ColorOrderIndex',1); +plot(freqs, 1./abs(squeeze(freqresp(W1, freqs, 'Hz'))), '--', 'DisplayName', '$|W_1|^{-1}$'); +plot(freqs, 1./abs(squeeze(freqresp(W2, freqs, 'Hz'))), '--', 'DisplayName', '$|W_2|^{-1}$'); +hold off; +set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); +xlabel('Frequency [Hz]'); ylabel('Magnitude'); +legend('location', 'northeast', 'FontSize', 8, 'NumColumns', 2); + +% Phase +ax2 = nexttile; +hold on; +set(gca,'ColorOrderIndex',1) +plot(freqs, 180/pi*phase(squeeze(freqresp(H1, freqs, 'Hz'))), '-'); +set(gca,'ColorOrderIndex',2) +plot(freqs, 180/pi*phase(squeeze(freqresp(H2, freqs, 'Hz'))), '-'); +hold off; +xlabel('Frequency [Hz]'); ylabel('Phase [deg]'); +set(gca, 'XScale', 'log'); +yticks([-180:90:180]); +ylim([-180, 180]); + +linkaxes([ax1,ax2],'x'); +xlim([freqs(1), freqs(end)]); +#+end_src + + * Impose a positive slope at DC or a negative slope at infinite frequency +** Introduction :ignore: + ** Matlab Init :noexport:ignore: #+begin_src matlab :tangle no :exports none :results silent :noweb yes :var current_dir=(file-name-directory buffer-file-name) - <> +<> #+end_src #+begin_src matlab :exports none :results silent :noweb yes - <> +<> +#+end_src + +#+begin_src matlab :tangle no +addpath('./matlab'); +addpath('./matlab/src'); +#+end_src + +#+begin_src matlab :exec no +addpath('./src'); #+end_src ** Manually shift zeros to the origin after synthesis @@ -1291,31 +1454,31 @@ We cannot impose that using the weight $W_2(s)$ as it would be improper. However, we may manually shift 2 of the low frequency zeros to the origin. #+begin_src matlab - n = 2; w0 = 2*pi*11; G0 = 1/10; G1 = 1000; Gc = 1/2; - W1 = (((1/w0)*sqrt((1-(G0/Gc)^(2/n))/(1-(Gc/G1)^(2/n)))*s + (G0/Gc)^(1/n))/((1/G1)^(1/n)*(1/w0)*sqrt((1-(G0/Gc)^(2/n))/(1-(Gc/G1)^(2/n)))*s + (1/Gc)^(1/n)))^n; +n = 2; w0 = 2*pi*11; G0 = 1/10; G1 = 1000; Gc = 1/2; +W1 = (((1/w0)*sqrt((1-(G0/Gc)^(2/n))/(1-(Gc/G1)^(2/n)))*s + (G0/Gc)^(1/n))/((1/G1)^(1/n)*(1/w0)*sqrt((1-(G0/Gc)^(2/n))/(1-(Gc/G1)^(2/n)))*s + (1/Gc)^(1/n)))^n; - n = 3; w0 = 2*pi*10; G0 = 1e4; G1 = 0.1; Gc = 1/2; - W2 = (((1/w0)*sqrt((1-(G0/Gc)^(2/n))/(1-(Gc/G1)^(2/n)))*s + (G0/Gc)^(1/n))/((1/G1)^(1/n)*(1/w0)*sqrt((1-(G0/Gc)^(2/n))/(1-(Gc/G1)^(2/n)))*s + (1/Gc)^(1/n)))^n; +n = 3; w0 = 2*pi*10; G0 = 1e4; G1 = 0.1; Gc = 1/2; +W2 = (((1/w0)*sqrt((1-(G0/Gc)^(2/n))/(1-(Gc/G1)^(2/n)))*s + (G0/Gc)^(1/n))/((1/G1)^(1/n)*(1/w0)*sqrt((1-(G0/Gc)^(2/n))/(1-(Gc/G1)^(2/n)))*s + (1/Gc)^(1/n)))^n; #+end_src #+begin_src matlab - P = [W1 -W1; - 0 W2; - 1 0]; +P = [W1 -W1; + 0 W2; + 1 0]; #+end_src And we do the $\mathcal{H}_\infty$ synthesis using the =hinfsyn= command: #+begin_src matlab - [H2, ~, gamma, ~] = hinfsyn(P, 1, 1,'TOLGAM', 0.001, 'METHOD', 'lmi', 'DISPLAY', 'on'); +[H2, ~, gamma, ~] = hinfsyn(P, 1, 1,'TOLGAM', 0.001, 'METHOD', 'lmi', 'DISPLAY', 'on'); #+end_src #+begin_src matlab - [z,p,k] = zpkdata(H2) +[z,p,k] = zpkdata(H2) #+end_src Looking at the zeros, we see two low frequency complex conjugate zeros. #+begin_src matlab :results output replace :exports results - z{1} +z{1} #+end_src #+RESULTS: @@ -1331,65 +1494,65 @@ ans = We manually put these zeros at the origin: #+begin_src matlab - z{1}([3,4]) = 0; +z{1}([3,4]) = 0; #+end_src And we create a modified filter $H_{2z}(s)$: #+begin_src matlab - H2z = zpk(z,p,k); +H2z = zpk(z,p,k); #+end_src And as usual, $H_{1z}(s)$ is defined as the complementary of $H_{2z}(s)$: #+begin_src matlab - H1z = 1 - H2z; +H1z = 1 - H2z; #+end_src The bode plots of $H_1(s)$, $H_2(s)$, $H_{1z}(s)$ and $H_{2z}(s)$ are shown in Figure [[fig:comp_filters_shift_zero]]. And we see that $H_{1z}(s)$ is slightly modified when setting the zeros at the origin for $H_{2z}(s)$. #+begin_src matlab :exports none - freqs = logspace(-2, 4, 1000); +freqs = logspace(-2, 4, 1000); - figure; - tiledlayout(3, 1, 'TileSpacing', 'None', 'Padding', 'None'); +figure; +tiledlayout(3, 1, 'TileSpacing', 'None', 'Padding', 'None'); - % Magnitude - ax1 = nexttile([2, 1]); - hold on; - set(gca,'ColorOrderIndex',1) - plot(freqs, abs(squeeze(freqresp(1-H2, freqs, 'Hz'))), '-', 'DisplayName', '$H_1$'); - plot(freqs, abs(squeeze(freqresp(H2, freqs, 'Hz'))), '-', 'DisplayName', '$H_2$'); +% Magnitude +ax1 = nexttile([2, 1]); +hold on; +set(gca,'ColorOrderIndex',1) +plot(freqs, abs(squeeze(freqresp(1-H2, freqs, 'Hz'))), '-', 'DisplayName', '$H_1$'); +plot(freqs, abs(squeeze(freqresp(H2, freqs, 'Hz'))), '-', 'DisplayName', '$H_2$'); - set(gca,'ColorOrderIndex',1) - plot(freqs, abs(squeeze(freqresp(1-H2z, freqs, 'Hz'))), '--', 'DisplayName', '$H_{1z}$'); - plot(freqs, abs(squeeze(freqresp(H2z, freqs, 'Hz'))), '--', 'DisplayName', '$H_{2z}$'); - set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); - hold off; - set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); - ylabel('Magnitude'); - set(gca, 'XTickLabel',[]); - legend('location', 'southeast', 'FontSize', 8, 'NumColumns', 2); +set(gca,'ColorOrderIndex',1) +plot(freqs, abs(squeeze(freqresp(1-H2z, freqs, 'Hz'))), '--', 'DisplayName', '$H_{1z}$'); +plot(freqs, abs(squeeze(freqresp(H2z, freqs, 'Hz'))), '--', 'DisplayName', '$H_{2z}$'); +set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); +hold off; +set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); +ylabel('Magnitude'); +set(gca, 'XTickLabel',[]); +legend('location', 'southeast', 'FontSize', 8, 'NumColumns', 2); - % Phase - ax2 = nexttile; - hold on; - set(gca,'ColorOrderIndex',1) - plot(freqs, 180/pi*phase(squeeze(freqresp(1-H2, freqs, 'Hz'))), '-'); - plot(freqs, 180/pi*phase(squeeze(freqresp(H2, freqs, 'Hz'))), '-'); - set(gca,'ColorOrderIndex',1) - plot(freqs, 180/pi*phase(squeeze(freqresp(H1z, freqs, 'Hz'))), '--'); - plot(freqs, 360+180/pi*phase(squeeze(freqresp(H2z, freqs, 'Hz'))), '--'); - hold off; - xlabel('Frequency [Hz]'); ylabel('Phase [deg]'); - set(gca, 'XScale', 'log'); - yticks([-180:90:180]); +% Phase +ax2 = nexttile; +hold on; +set(gca,'ColorOrderIndex',1) +plot(freqs, 180/pi*phase(squeeze(freqresp(1-H2, freqs, 'Hz'))), '-'); +plot(freqs, 180/pi*phase(squeeze(freqresp(H2, freqs, 'Hz'))), '-'); +set(gca,'ColorOrderIndex',1) +plot(freqs, 180/pi*phase(squeeze(freqresp(H1z, freqs, 'Hz'))), '--'); +plot(freqs, 360+180/pi*phase(squeeze(freqresp(H2z, freqs, 'Hz'))), '--'); +hold off; +xlabel('Frequency [Hz]'); ylabel('Phase [deg]'); +set(gca, 'XScale', 'log'); +yticks([-180:90:180]); - linkaxes([ax1,ax2],'x'); - xlim([freqs(1), freqs(end)]); +linkaxes([ax1,ax2],'x'); +xlim([freqs(1), freqs(end)]); #+end_src #+begin_src matlab :tangle no :exports results :results file replace - exportFig('figs/comp_filters_shift_zero.pdf', 'width', 'wide', 'height', 'tall'); +exportFig('figs/comp_filters_shift_zero.pdf', 'width', 'wide', 'height', 'tall'); #+end_src #+name: fig:comp_filters_shift_zero @@ -1447,87 +1610,87 @@ For $H_1(s)$ nothing is changed: $H_1(s) = 1 - H_2(s)$. The weighting functions are defined as usual: #+begin_src matlab - n = 2; w0 = 2*pi*11; G0 = 1/10; G1 = 1000; Gc = 1/2; - W1 = (((1/w0)*sqrt((1-(G0/Gc)^(2/n))/(1-(Gc/G1)^(2/n)))*s + (G0/Gc)^(1/n))/((1/G1)^(1/n)*(1/w0)*sqrt((1-(G0/Gc)^(2/n))/(1-(Gc/G1)^(2/n)))*s + (1/Gc)^(1/n)))^n; +n = 2; w0 = 2*pi*11; G0 = 1/10; G1 = 1000; Gc = 1/2; +W1 = (((1/w0)*sqrt((1-(G0/Gc)^(2/n))/(1-(Gc/G1)^(2/n)))*s + (G0/Gc)^(1/n))/((1/G1)^(1/n)*(1/w0)*sqrt((1-(G0/Gc)^(2/n))/(1-(Gc/G1)^(2/n)))*s + (1/Gc)^(1/n)))^n; - n = 3; w0 = 2*pi*10; G0 = 1e4; G1 = 0.1; Gc = 1/2; - W2 = (((1/w0)*sqrt((1-(G0/Gc)^(2/n))/(1-(Gc/G1)^(2/n)))*s + (G0/Gc)^(1/n))/((1/G1)^(1/n)*(1/w0)*sqrt((1-(G0/Gc)^(2/n))/(1-(Gc/G1)^(2/n)))*s + (1/Gc)^(1/n)))^n; +n = 3; w0 = 2*pi*10; G0 = 1e4; G1 = 0.1; Gc = 1/2; +W2 = (((1/w0)*sqrt((1-(G0/Gc)^(2/n))/(1-(Gc/G1)^(2/n)))*s + (G0/Gc)^(1/n))/((1/G1)^(1/n)*(1/w0)*sqrt((1-(G0/Gc)^(2/n))/(1-(Gc/G1)^(2/n)))*s + (1/Gc)^(1/n)))^n; #+end_src The wanted feature here is a +2 slope at low frequency. For that, we use an high pass filter with a slope of +2 at low frequency. #+begin_src matlab - w0 = 2*pi*50; - H2w = (s/w0/(s/w0+1))^2; +w0 = 2*pi*50; +H2w = (s/w0/(s/w0+1))^2; #+end_src We define the generalized plant as shown in Figure [[fig:h_infinity_arch_H2_feature]]. #+begin_src matlab - P = [W1 -W1; - 0 W2; - H2w 0]; +P = [W1 -W1; + 0 W2; + H2w 0]; #+end_src And we do the $\mathcal{H}_\infty$ synthesis using the =hinfsyn= command. #+begin_src matlab :results output replace :exports both - [H2p, ~, gamma, ~] = hinfsyn(P, 1, 1,'TOLGAM', 0.001, 'METHOD', 'lmi', 'DISPLAY', 'on'); +[H2p, ~, gamma, ~] = hinfsyn(P, 1, 1,'TOLGAM', 0.001, 'METHOD', 'lmi', 'DISPLAY', 'on'); #+end_src Finally, we define $H_2(s)$ as the product of the synthesized filter and the wanted "feature": #+begin_src matlab - H2 = H2p*H2w; +H2 = H2p*H2w; #+end_src And we define $H_1(s)$ to be the complementary of $H_2(s)$: #+begin_src matlab - H1 = 1 - H2; +H1 = 1 - H2; #+end_src The obtained complementary filters are shown in Figure [[fig:comp_filters_H2_feature]]. #+begin_src matlab :exports none - freqs = logspace(-3, 3, 1000); +freqs = logspace(-3, 3, 1000); - figure; - tiledlayout(3, 1, 'TileSpacing', 'None', 'Padding', 'None'); +figure; +tiledlayout(3, 1, 'TileSpacing', 'None', 'Padding', 'None'); - % Magnitude - ax1 = nexttile([2, 1]); - hold on; - set(gca,'ColorOrderIndex',1) - plot(freqs, 1./abs(squeeze(freqresp(W1, freqs, 'Hz'))), '--', 'DisplayName', '$w_1$'); - set(gca,'ColorOrderIndex',2) - plot(freqs, 1./abs(squeeze(freqresp(W2, freqs, 'Hz'))), '--', 'DisplayName', '$w_2$'); +% Magnitude +ax1 = nexttile([2, 1]); +hold on; +set(gca,'ColorOrderIndex',1) +plot(freqs, 1./abs(squeeze(freqresp(W1, freqs, 'Hz'))), '--', 'DisplayName', '$w_1$'); +set(gca,'ColorOrderIndex',2) +plot(freqs, 1./abs(squeeze(freqresp(W2, freqs, 'Hz'))), '--', 'DisplayName', '$w_2$'); - set(gca,'ColorOrderIndex',1) - plot(freqs, abs(squeeze(freqresp(H1, freqs, 'Hz'))), '-', 'DisplayName', '$H_1$'); - set(gca,'ColorOrderIndex',2) - plot(freqs, abs(squeeze(freqresp(H2, freqs, 'Hz'))), '-', 'DisplayName', '$H_2$'); - set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); - hold off; - set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); - ylabel('Magnitude'); - set(gca, 'XTickLabel',[]); - legend('location', 'southeast', 'FontSize', 8, 'NumColumns', 2); +set(gca,'ColorOrderIndex',1) +plot(freqs, abs(squeeze(freqresp(H1, freqs, 'Hz'))), '-', 'DisplayName', '$H_1$'); +set(gca,'ColorOrderIndex',2) +plot(freqs, abs(squeeze(freqresp(H2, freqs, 'Hz'))), '-', 'DisplayName', '$H_2$'); +set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); +hold off; +set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); +ylabel('Magnitude'); +set(gca, 'XTickLabel',[]); +legend('location', 'southeast', 'FontSize', 8, 'NumColumns', 2); - % Phase - ax2 = nexttile; - hold on; - set(gca,'ColorOrderIndex',1) - plot(freqs, 180/pi*phase(squeeze(freqresp(H1, freqs, 'Hz'))), '-'); - set(gca,'ColorOrderIndex',2) - plot(freqs, 360+180/pi*phase(squeeze(freqresp(H2, freqs, 'Hz'))), '-'); - hold off; - xlabel('Frequency [Hz]'); ylabel('Phase [deg]'); - set(gca, 'XScale', 'log'); - yticks([-180:90:180]); +% Phase +ax2 = nexttile; +hold on; +set(gca,'ColorOrderIndex',1) +plot(freqs, 180/pi*phase(squeeze(freqresp(H1, freqs, 'Hz'))), '-'); +set(gca,'ColorOrderIndex',2) +plot(freqs, 360+180/pi*phase(squeeze(freqresp(H2, freqs, 'Hz'))), '-'); +hold off; +xlabel('Frequency [Hz]'); ylabel('Phase [deg]'); +set(gca, 'XScale', 'log'); +yticks([-180:90:180]); - linkaxes([ax1,ax2],'x'); - xlim([freqs(1), freqs(end)]); +linkaxes([ax1,ax2],'x'); +xlim([freqs(1), freqs(end)]); #+end_src #+begin_src matlab :tangle no :exports results :results file replace - exportFig('figs/comp_filters_H2_feature.pdf', 'width', 'wide', 'height', 'tall'); +exportFig('figs/comp_filters_H2_feature.pdf', 'width', 'wide', 'height', 'tall'); #+end_src #+name: fig:comp_filters_H2_feature @@ -1542,85 +1705,85 @@ The used technique is the same as in the previous section, and the generalized p The weights are defined as usual. #+begin_src matlab - n = 3; w0 = 2*pi*10; G0 = 1000; G1 = 0.1; Gc = 1/2; - W1 = (((1/w0)*sqrt((1-(G0/Gc)^(2/n))/(1-(Gc/G1)^(2/n)))*s + (G0/Gc)^(1/n))/((1/G1)^(1/n)*(1/w0)*sqrt((1-(G0/Gc)^(2/n))/(1-(Gc/G1)^(2/n)))*s + (1/Gc)^(1/n)))^n; +n = 3; w0 = 2*pi*10; G0 = 1000; G1 = 0.1; Gc = 1/2; +W1 = (((1/w0)*sqrt((1-(G0/Gc)^(2/n))/(1-(Gc/G1)^(2/n)))*s + (G0/Gc)^(1/n))/((1/G1)^(1/n)*(1/w0)*sqrt((1-(G0/Gc)^(2/n))/(1-(Gc/G1)^(2/n)))*s + (1/Gc)^(1/n)))^n; - n = 2; w0 = 2*pi*11; G0 = 1/10; G1 = 1000; Gc = 1/2; - W2 = (((1/w0)*sqrt((1-(G0/Gc)^(2/n))/(1-(Gc/G1)^(2/n)))*s + (G0/Gc)^(1/n))/((1/G1)^(1/n)*(1/w0)*sqrt((1-(G0/Gc)^(2/n))/(1-(Gc/G1)^(2/n)))*s + (1/Gc)^(1/n)))^n; +n = 2; w0 = 2*pi*11; G0 = 1/10; G1 = 1000; Gc = 1/2; +W2 = (((1/w0)*sqrt((1-(G0/Gc)^(2/n))/(1-(Gc/G1)^(2/n)))*s + (G0/Gc)^(1/n))/((1/G1)^(1/n)*(1/w0)*sqrt((1-(G0/Gc)^(2/n))/(1-(Gc/G1)^(2/n)))*s + (1/Gc)^(1/n)))^n; #+end_src This time, the feature is a low pass filter with a slope of -2 at high frequency. #+begin_src matlab - H2w = 1/(s/(2*pi*10) + 1)^2; +H2w = 1/(s/(2*pi*10) + 1)^2; #+end_src The generalized plant is defined: #+begin_src matlab - P = [W1 -W1; - 0 W2; - H2w 0]; +P = [W1 -W1; + 0 W2; + H2w 0]; #+end_src And we do the $\mathcal{H}_\infty$ synthesis using the =hinfsyn= command. #+begin_src matlab :results output replace :exports both - [H2p, ~, gamma, ~] = hinfsyn(P, 1, 1,'TOLGAM', 0.001, 'METHOD', 'lmi', 'DISPLAY', 'on'); +[H2p, ~, gamma, ~] = hinfsyn(P, 1, 1,'TOLGAM', 0.001, 'METHOD', 'lmi', 'DISPLAY', 'on'); #+end_src The feature is added to the synthesized filter: #+begin_src matlab - H2 = H2p*H2w; +H2 = H2p*H2w; #+end_src And $H_1(s)$ is defined as follows: #+begin_src matlab - H1 = 1 - H2; +H1 = 1 - H2; #+end_src The obtained complementary filters are shown in Figure [[fig:comp_filters_H2_feature_neg_slope]]. #+begin_src matlab :exports none - freqs = logspace(-1, 4, 1000); +freqs = logspace(-1, 4, 1000); - figure; - tiledlayout(3, 1, 'TileSpacing', 'None', 'Padding', 'None'); +figure; +tiledlayout(3, 1, 'TileSpacing', 'None', 'Padding', 'None'); - % Magnitude - ax1 = nexttile([2, 1]); - hold on; - set(gca,'ColorOrderIndex',1) - plot(freqs, 1./abs(squeeze(freqresp(W1, freqs, 'Hz'))), '--', 'DisplayName', '$w_1$'); - set(gca,'ColorOrderIndex',2) - plot(freqs, 1./abs(squeeze(freqresp(W2, freqs, 'Hz'))), '--', 'DisplayName', '$w_2$'); +% Magnitude +ax1 = nexttile([2, 1]); +hold on; +set(gca,'ColorOrderIndex',1) +plot(freqs, 1./abs(squeeze(freqresp(W1, freqs, 'Hz'))), '--', 'DisplayName', '$w_1$'); +set(gca,'ColorOrderIndex',2) +plot(freqs, 1./abs(squeeze(freqresp(W2, freqs, 'Hz'))), '--', 'DisplayName', '$w_2$'); - set(gca,'ColorOrderIndex',1) - plot(freqs, abs(squeeze(freqresp(H1, freqs, 'Hz'))), '-', 'DisplayName', '$H_1$'); - set(gca,'ColorOrderIndex',2) - plot(freqs, abs(squeeze(freqresp(H2, freqs, 'Hz'))), '-', 'DisplayName', '$H_2$'); - set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); - hold off; - set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); - ylabel('Magnitude'); - set(gca, 'XTickLabel',[]); - legend('location', 'southeast', 'FontSize', 8, 'NumColumns', 2); +set(gca,'ColorOrderIndex',1) +plot(freqs, abs(squeeze(freqresp(H1, freqs, 'Hz'))), '-', 'DisplayName', '$H_1$'); +set(gca,'ColorOrderIndex',2) +plot(freqs, abs(squeeze(freqresp(H2, freqs, 'Hz'))), '-', 'DisplayName', '$H_2$'); +set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); +hold off; +set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log'); +ylabel('Magnitude'); +set(gca, 'XTickLabel',[]); +legend('location', 'southeast', 'FontSize', 8, 'NumColumns', 2); - % Phase - ax2 = nexttile; - hold on; - set(gca,'ColorOrderIndex',1) - plot(freqs, 180/pi*phase(squeeze(freqresp(H1, freqs, 'Hz'))), '-'); - set(gca,'ColorOrderIndex',2) - plot(freqs, 180/pi*phase(squeeze(freqresp(H2, freqs, 'Hz'))), '-'); - hold off; - xlabel('Frequency [Hz]'); ylabel('Phase [deg]'); - set(gca, 'XScale', 'log'); - yticks([-180:90:180]); +% Phase +ax2 = nexttile; +hold on; +set(gca,'ColorOrderIndex',1) +plot(freqs, 180/pi*phase(squeeze(freqresp(H1, freqs, 'Hz'))), '-'); +set(gca,'ColorOrderIndex',2) +plot(freqs, 180/pi*phase(squeeze(freqresp(H2, freqs, 'Hz'))), '-'); +hold off; +xlabel('Frequency [Hz]'); ylabel('Phase [deg]'); +set(gca, 'XScale', 'log'); +yticks([-180:90:180]); - linkaxes([ax1,ax2],'x'); - xlim([freqs(1), freqs(end)]); +linkaxes([ax1,ax2],'x'); +xlim([freqs(1), freqs(end)]); #+end_src #+begin_src matlab :tangle no :exports results :results file replace - exportFig('figs/comp_filters_H2_feature_neg_slope.pdf', 'width', 'wide', 'height', 'tall'); +exportFig('figs/comp_filters_H2_feature_neg_slope.pdf', 'width', 'wide', 'height', 'tall'); #+end_src #+name: fig:comp_filters_H2_feature_neg_slope @@ -1632,3 +1795,142 @@ The obtained complementary filters are shown in Figure [[fig:comp_filters_H2_fea * Bibliography :ignore: bibliographystyle:unsrt bibliography:ref.bib + +* Functions +** =generateWF=: Generate Weighting Functions +:PROPERTIES: +:header-args:matlab+: :tangle matlab/src/generateWF.m +:header-args:matlab+: :comments none :mkdirp yes :eval no +:END: +<> + +This Matlab function is accessible [[file:matlab/src/generateWF.m][here]]. + +*** Function description +:PROPERTIES: +:UNNUMBERED: t +:END: + +#+begin_src matlab +function [W] = generateWF(args) +% createWeight - +% +% Syntax: [W] = generateWeight(args) +% +% Inputs: +% - n - Weight Order (integer) +% - G0 - Low frequency Gain +% - G1 - High frequency Gain +% - Gc - Gain of the weight at frequency w0 +% - w0 - Frequency at which |W(j w0)| = Gc [rad/s] +% +% Outputs: +% - W - Generated Weight +#+end_src + +*** Optional Parameters +:PROPERTIES: +:UNNUMBERED: t +:END: +#+begin_src matlab +arguments + args.n (1,1) double {mustBeInteger, mustBePositive} = 1 + args.G0 (1,1) double {mustBeNumeric, mustBePositive} = 0.1 + args.Ginf (1,1) double {mustBeNumeric, mustBePositive} = 10 + args.Gc (1,1) double {mustBeNumeric, mustBePositive} = 1 + args.w0 (1,1) double {mustBeNumeric, mustBePositive} = 1 +end + +mustBeBetween(args.G0, args.Gc, args.Ginf); +#+end_src + +*** Generate the Weighting function +:PROPERTIES: +:UNNUMBERED: t +:END: +#+begin_src matlab +s = zpk('s'); +#+end_src + +#+begin_src matlab +W = (((1/args.w0)*sqrt((1-(args.G0/args.Gc)^(2/args.n))/(1-(args.Gc/args.Ginf)^(2/args.n)))*s + ... + (args.G0/args.Gc)^(1/args.n))/... + ((1/args.Ginf)^(1/args.n)*(1/args.w0)*sqrt((1-(args.G0/args.Gc)^(2/args.n))/(1-(args.Gc/args.Ginf)^(2/args.n)))*s + ... + (1/args.Gc)^(1/args.n)))^args.n; +#+end_src + + +*** Verification of the $G_0$, $G_c$ and $G_\infty$ gains +:PROPERTIES: +:UNNUMBERED: t +:END: +#+begin_src matlab +% Custom validation function +function mustBeBetween(a,b,c) + if ~((a > b && b > c) || (c > b && b > a)) + eid = 'createWeight:inputError'; + msg = 'Gc should be between G0 and Ginf.'; + throwAsCaller(MException(eid,msg)) + end +#+end_src + +** =generateCF=: Generate Complementary Filters +:PROPERTIES: +:header-args:matlab+: :tangle matlab/src/generateCF.m +:header-args:matlab+: :comments none :mkdirp yes :eval no +:END: +<> + +This Matlab function is accessible [[file:matlab/src/generateCF.m][here]]. + +*** Function description +:PROPERTIES: +:UNNUMBERED: t +:END: + +#+begin_src matlab +function [H1, H2] = generateCF(W1, W2, args) +% createWeight - +% +% Syntax: [W] = generateCF(args) +% +% Inputs: +% - W1 - Weighting Function for H1 +% - W2 - Weighting Function for H2 +% - args - +% +% Outputs: +% - H1 - Generated H1 Filter +% - H2 - Generated H2 Filter +#+end_src + +*** Optional Parameters +:PROPERTIES: +:UNNUMBERED: t +:END: +#+begin_src matlab +arguments + W1 + W2 + args.method char {mustBeMember(args.method,{'lmi', 'ric'})} = 'ric' + args.display char {mustBeMember(args.display,{'on', 'off'})} = 'on' +end +#+end_src + +*** H-Infinity Synthesis +:PROPERTIES: +:UNNUMBERED: t +:END: +#+begin_src matlab +P = [W1 -W1; + 0 W2; + 1 0]; +#+end_src + +#+begin_src matlab :results output replace :exports both +[H2, ~, gamma, ~] = hinfsyn(P, 1, 1,'TOLGAM', 0.001, 'METHOD', args.method, 'DISPLAY', args.display); +#+end_src + +#+begin_src matlab +H1 = 1 - H2; +#+end_src diff --git a/matlab/matlab/src/generateCF.m b/matlab/matlab/src/generateCF.m new file mode 100644 index 0000000..1ccecfd --- /dev/null +++ b/matlab/matlab/src/generateCF.m @@ -0,0 +1,28 @@ +function [H1, H2] = generateCF(W1, W2, args) +% createWeight - +% +% Syntax: [W] = generateCF(args) +% +% Inputs: +% - W1 - Weighting Function for H1 +% - W2 - Weighting Function for H2 +% - args - +% +% Outputs: +% - H1 - Generated H1 Filter +% - H2 - Generated H2 Filter + +arguments + W1 + W2 + args.method char {mustBeMember(args.method,{'lmi', 'ric'})} = 'ric' + args.display char {mustBeMember(args.display,{'on', 'off'})} = 'on' +end + +P = [W1 -W1; + 0 W2; + 1 0]; + +[H2, ~, gamma, ~] = hinfsyn(P, 1, 1,'TOLGAM', 0.001, 'METHOD', args.method, 'DISPLAY', args.display); + +H1 = 1 - H2; diff --git a/matlab/matlab/src/generateWF.m b/matlab/matlab/src/generateWF.m new file mode 100644 index 0000000..7276dd1 --- /dev/null +++ b/matlab/matlab/src/generateWF.m @@ -0,0 +1,39 @@ +function [W] = generateWF(args) +% createWeight - +% +% Syntax: [W] = generateWeight(args) +% +% Inputs: +% - n - Weight Order (integer) +% - G0 - Low frequency Gain +% - G1 - High frequency Gain +% - Gc - Gain of the weight at frequency w0 +% - w0 - Frequency at which |W(j w0)| = Gc [rad/s] +% +% Outputs: +% - W - Generated Weight + +arguments + args.n (1,1) double {mustBeInteger, mustBePositive} = 1 + args.G0 (1,1) double {mustBeNumeric, mustBePositive} = 0.1 + args.Ginf (1,1) double {mustBeNumeric, mustBePositive} = 10 + args.Gc (1,1) double {mustBeNumeric, mustBePositive} = 1 + args.w0 (1,1) double {mustBeNumeric, mustBePositive} = 1 +end + +mustBeBetween(args.G0, args.Gc, args.Ginf); + +s = zpk('s'); + +W = (((1/args.w0)*sqrt((1-(args.G0/args.Gc)^(2/args.n))/(1-(args.Gc/args.Ginf)^(2/args.n)))*s + ... + (args.G0/args.Gc)^(1/args.n))/... + ((1/args.Ginf)^(1/args.n)*(1/args.w0)*sqrt((1-(args.G0/args.Gc)^(2/args.n))/(1-(args.Gc/args.Ginf)^(2/args.n)))*s + ... + (1/args.Gc)^(1/args.n)))^args.n; + +% Custom validation function +function mustBeBetween(a,b,c) + if ~((a > b && b > c) || (c > b && b > a)) + eid = 'createWeight:inputError'; + msg = 'Gc should be between G0 and Ginf.'; + throwAsCaller(MException(eid,msg)) + end