tdehaeze/dehaeze21_desig_compl_filte
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

Rework Matlab figure size + Start section 4

Browse Source
This commit is contained in:
Thomas Dehaeze
2021-05-20 16:27:51 +02:00
parent c114cdbc58
commit 1fda0d3383
26 changed files with 1305 additions and 887 deletions
Download Patch File Download Diff File Expand all files Collapse all files

106
matlab/index.org
View File

@@ -292,7 +292,7 @@ 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', 600);
#+end_src
#+name: fig:hinf_filters_results
@@ -492,14 +492,14 @@ 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]);
yticks([-180:90:180]); ylim([-220, 220]);
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', 600);
#+end_src
#+name: fig:three_complementary_filters_results
@@ -739,30 +739,39 @@ 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,'ColorOrderIndex',1)
plot(w, abs(1-H), '-', 'DisplayName', '$L_1$');
plot([0.1, 10], [0.045, 0.045], 'k:', 'DisplayName', 'Spec. on $L_1$');
set(gca,'ColorOrderIndex',2)
plot(w, abs(H), '-', 'DisplayName', '$H_1$');
plot([0.0001, 0.008], [8e-3, 8e-3], 'k--', 'DisplayName', 'Spec. on $H_1$');
plot([0.008 0.04], [8e-3, 1], 'k--', 'HandleVisibility', 'off');
plot([0.04 0.1], [3, 3], 'k--', 'HandleVisibility', 'off');
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
ylabel('Magnitude');
set(gca, 'XTickLabel',[]);
ylim([5e-3, 5]);
set(gca, 'XTickLabel',[]); ylabel('Magnitude');
hold off;
ylim([5e-3, 10]);
leg = legend('location', 'southeast', 'FontSize', 8, 'NumColumns', 2);
leg.ItemTokenSize(1) = 16;
% Phase
ax2 = nexttile;
hold on;
plot(w, 180/pi*unwrap(angle(H)), 'k-');
plot(w, 180/pi*unwrap(angle(1-H)), 'k--');
plot(w, 180/pi*unwrap(angle(1-H)), '-');
plot(w, 180/pi*unwrap(angle(H)), '-');
hold off;
xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
set(gca, 'XScale', 'log');
yticks([-450:90:180]); ylim([-450, 200]);
yticks([-360:180:180]); ylim([-380, 200]);
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', 'tall');
exportFig('figs/fir_filter_ligo.pdf', 'width', 'wide', 'height', 600);
#+end_src
#+name: fig:fir_filter_ligo
@@ -812,22 +821,22 @@ wL = minreal(ss(wL));
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(freqs, abs(squeeze(freqresp(inv(wL), freqs, 'Hz'))), '-', 'DisplayName', '$|W_L|^{-1}$');
plot([0.1, 10], [0.045, 0.045], 'k:', 'DisplayName', 'Spec. on $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');
set(gca,'ColorOrderIndex',2);
plot(freqs, abs(squeeze(freqresp(inv(wH), freqs, 'Hz'))), '-', 'DisplayName', '$|W_H|^{-1}$');
plot([0.0001, 0.008], [8e-3, 8e-3], 'k--', 'DisplayName', 'Spec. on $H_1$');
plot([0.008 0.04], [8e-3, 1], 'k--', 'HandleVisibility', 'off');
plot([0.04 0.1], [3, 3], '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]);
leg = legend('location', 'southeast', 'FontSize', 8);
leg.ItemTokenSize(1) = 18;
ylim([5e-3, 10]);
leg = legend('location', 'southeast', 'FontSize', 8, 'NumColumns', 2);
leg.ItemTokenSize(1) = 16;
#+end_src
#+begin_src matlab :tangle no :exports results :results file replace
@@ -905,27 +914,22 @@ The bode plot of the obtained filters as shown on figure [[fig:hinf_synthesis_li
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');
plot(freqs, abs(squeeze(freqresp(Hl, freqs, 'Hz'))), '-', 'DisplayName', '$L_1^\prime$');
plot([0.1, 10], [0.045, 0.045], 'k:', 'DisplayName', 'Spec. on $L_1$');
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$');
plot(freqs, abs(squeeze(freqresp(Hh, freqs, 'Hz'))), '-', 'DisplayName', '$H_1^\prime$');
plot([0.0001, 0.008], [8e-3, 8e-3], 'k--', 'DisplayName', 'Spec. on $H_1$');
plot([0.008 0.04], [8e-3, 1], 'k--', 'HandleVisibility', 'off');
plot([0.04 0.1], [3, 3], '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]);
leg = legend('location', 'southeast', 'FontSize', 8);
leg.ItemTokenSize(1) = 18;
ylim([5e-3, 10]);
leg = legend('location', 'southeast', 'FontSize', 8, 'NumColumns', 2);
leg.ItemTokenSize(1) = 16;
#+end_src
#+begin_src matlab :tangle no :exports results :results file replace
@@ -948,18 +952,18 @@ tiledlayout(3, 1, 'TileSpacing', 'None', 'Padding', 'None');
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$');
'DisplayName', '$L_1(s)$ - $\mathcal{H}_\infty$');
set(gca,'ColorOrderIndex',2);
plot(freqs, abs(squeeze(freqresp(Hh, freqs, 'Hz'))), '-', ...
'DisplayName', '$H_1(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');
'DisplayName', '$L_1(s)$ - FIR');
set(gca,'ColorOrderIndex',2);
plot(w, abs(H), '--', ...
'DisplayName', '$H_1(s)$ - FIR');
hold off;
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
ylabel('Magnitude');
@@ -972,25 +976,25 @@ leg.ItemTokenSize(1) = 16;
ax2 = nexttile;
hold on;
set(gca,'ColorOrderIndex',1);
plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(Hh, freqs, 'Hz')))), '-');
set(gca,'ColorOrderIndex',2);
plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(Hl, freqs, 'Hz')))), '-');
set(gca,'ColorOrderIndex',2);
plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(Hh, freqs, 'Hz')))), '-');
set(gca,'ColorOrderIndex',1);
plot(w, 180/pi*unwrap(angle(H)), '--');
set(gca,'ColorOrderIndex',2);
plot(w, 180/pi*unwrap(angle(1-H)), '--');
set(gca,'ColorOrderIndex',2);
plot(w, 180/pi*unwrap(angle(H)), '--');
set(gca, 'XScale', 'log');
xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
hold off;
yticks([-450:90:180]); ylim([-450, 200]);
yticks([-360:180:180]); ylim([-380, 200]);
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', 600);
#+end_src
#+name: fig:comp_fir_ligo_hinf

61
matlab/matlab/comp_filters_ligo.m
View File

@@ -6,6 +6,8 @@ s = zpk('s');
freqs = logspace(-3, 0, 1000);
addpath('./src');
% Specifications
% The specifications for the filters are:
% 1. From $0$ to $0.008\text{ Hz}$,the magnitude of the filter’s transfer function should be less than or equal to $8 \times 10^{-3}$
@@ -31,7 +33,8 @@ xlabel('Frequency [Hz]'); ylabel('Magnitude');
hold off;
xlim([freqs(1), freqs(end)]);
ylim([1e-4, 10]);
legend('location', 'northeast');
leg = legend('location', 'southeast', 'FontSize', 8);
leg.ItemTokenSize(1) = 18;
% FIR Filter
% We here try to implement the FIR complementary filter synthesis as explained in cite:hua05_low_ligo.
@@ -184,8 +187,10 @@ w = [w1 w2 w3 w4];
H = [exp(-j*kron(w'.*2*pi,[0:n-1]))]*h;
figure;
tiledlayout(3, 1, 'TileSpacing', 'None', 'Padding', 'None');
ax1 = subplot(2,1,1);
% Magnitude
ax1 = nexttile([2, 1]);
hold on;
plot(w, abs(H), 'k-');
plot(w, abs(1-H), 'k--');
@@ -195,18 +200,18 @@ ylabel('Magnitude');
set(gca, 'XTickLabel',[]);
ylim([5e-3, 5]);
ax2 = subplot(2,1,2);
% Phase
ax2 = nexttile;
hold on;
plot(w, 180/pi*angle(H), 'k-');
plot(w, 180/pi*angle(1-H), 'k--');
plot(w, 180/pi*unwrap(angle(H)), 'k-');
plot(w, 180/pi*unwrap(angle(1-H)), 'k--');
hold off;
xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
set(gca, 'XScale', 'log');
yticks([-540:90:360]);
yticks([-450:90:180]); ylim([-450, 200]);
linkaxes([ax1,ax2],'x');
xlim([1e-3, 1]);
xticks([0.01, 0.1, 1, 10, 100, 1000]);
% Weights
% We design weights that will be used for the $\mathcal{H}_\infty$ synthesis of the complementary filters.
@@ -260,7 +265,8 @@ xlabel('Frequency [Hz]'); ylabel('Magnitude');
hold off;
xlim([freqs(1), freqs(end)]);
ylim([1e-3, 10]);
legend('location', 'southeast');
leg = legend('location', 'southeast', 'FontSize', 8);
leg.ItemTokenSize(1) = 18;
% H-Infinity Synthesis
% We define the generalized plant as shown on figure [[fig:h_infinity_robst_fusion]].
@@ -349,47 +355,56 @@ xlabel('Frequency [Hz]'); ylabel('Magnitude');
hold off;
xlim([freqs(1), freqs(end)]);
ylim([1e-3, 10]);
legend('location', 'southeast');
leg = legend('location', 'southeast', 'FontSize', 8);
leg.ItemTokenSize(1) = 18;
% Compare FIR and H-Infinity Filters
% 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]].
figure;
ax1 = subplot(2,1,1);
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'))), '-');
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'))), '-');
plot(freqs, abs(squeeze(freqresp(Hl, freqs, 'Hz'))), '-', ...
'DisplayName', '$H_L(s)$ - $\mathcal{H}_\infty$');
set(gca,'ColorOrderIndex',1);
plot(w, abs(H), '--');
plot(w, abs(H), '--', ...
'DisplayName', '$H_H(s)$ - FIR');
set(gca,'ColorOrderIndex',2);
plot(w, abs(1-H), '--');
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([1e-3, 10]);
ylim([5e-3, 10]);
leg = legend('location', 'southeast', 'FontSize', 8, 'NumColumns', 2);
leg.ItemTokenSize(1) = 16;
ax2 = subplot(2,1,2);
% Phase
ax2 = nexttile;
hold on;
set(gca,'ColorOrderIndex',1);
plot(freqs, 180/pi*angle(squeeze(freqresp(Hh, freqs, 'Hz'))), '-', 'DisplayName', '$\mathcal{H}_\infty$ filters');
plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(Hh, freqs, 'Hz')))), '-');
set(gca,'ColorOrderIndex',2);
plot(freqs, 180/pi*angle(squeeze(freqresp(Hl, freqs, 'Hz'))), '-', 'HandleVisibility', 'off');
plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(Hl, freqs, 'Hz')))), '-');
set(gca,'ColorOrderIndex',1);
plot(w, 180/pi*angle(H), '--', 'DisplayName', 'FIR filters');
plot(w, 180/pi*unwrap(angle(H)), '--');
set(gca,'ColorOrderIndex',2);
plot(w, 180/pi*angle(1-H), '--', 'HandleVisibility', 'off');
plot(w, 180/pi*unwrap(angle(1-H)), '--');
set(gca, 'XScale', 'log');
xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
hold off;
yticks([-540:90:360]);
legend('location', 'northeast');
yticks([-450:90:180]); ylim([-450, 200]);
linkaxes([ax1,ax2],'x');
xlim([freqs(1), freqs(end)]);
xticks([0.001, 0.01, 0.1, 1]);

109
matlab/matlab/h_inf_synthesis_complementary_filters.m
View File

@@ -6,6 +6,8 @@ s = zpk('s');
freqs = logspace(-1, 3, 1000);
addpath('./src');
% Design of Weighting Function
% A formula is proposed to help the design of the weighting functions:
% \begin{equation}
@@ -24,15 +26,48 @@ freqs = logspace(-1, 3, 1000);
% The general shape of a weighting function generated using the formula is shown in figure [[fig:weight_formula]].
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;
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([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/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/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]);
% #+name: fig:weight_formula
% #+caption: Amplitude of the proposed formula for the weighting functions
% [[file:figs-tikz/weight_formula.png]]
% #+caption: Gain of the Weighting Function formula
% #+RESULTS:
% [[file:figs/weight_formula.png]]
n = 2; w0 = 2*pi*11; G0 = 1/10; G1 = 1000; Gc = 1/2;
n = 2; w0 = 2*pi*10; G0 = 1/10; G1 = 1000; Gc = 0.45;
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;
n = 3; w0 = 2*pi*10; G0 = 1000; G1 = 0.1; Gc = 0.45;
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;
figure;
@@ -45,9 +80,10 @@ set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
xlabel('Frequency [Hz]'); ylabel('Magnitude');
hold off;
xlim([freqs(1), freqs(end)]);
ylim([5e-4, 20]);
ylim([1e-4, 20]);
xticks([0.1, 1, 10, 100, 1000]);
legend('location', 'northeast');
leg = legend('location', 'southeast', 'FontSize', 8);
leg.ItemTokenSize(1) = 18;
% H-Infinity Synthesis
% We define the generalized plant $P$ on matlab.
@@ -67,35 +103,28 @@ P = [W1 -W1;
% #+RESULTS:
% #+begin_example
% [H2, ~, gamma, ~] = hinfsyn(P, 1, 1,'TOLGAM', 0.001, 'METHOD', 'ric', 'DISPLAY', 'on');
% Resetting value of Gamma min based on D_11, D_12, D_21 terms
% Test bounds: 0.1000 < gamma <= 1050.0000
% Test bounds: 0.3223 <= gamma <= 1000
% gamma hamx_eig xinf_eig hamy_eig yinf_eig nrho_xy p/f
% 1.050e+03 2.8e+01 2.4e-07 4.1e+00 0.0e+00 0.0000 p
% 525.050 2.8e+01 2.4e-07 4.1e+00 0.0e+00 0.0000 p
% 262.575 2.8e+01 2.4e-07 4.1e+00 0.0e+00 0.0000 p
% 131.337 2.8e+01 2.4e-07 4.1e+00 -1.0e-13 0.0000 p
% 65.719 2.8e+01 2.4e-07 4.1e+00 -9.5e-14 0.0000 p
% 32.909 2.8e+01 2.4e-07 4.1e+00 0.0e+00 0.0000 p
% 16.505 2.8e+01 2.4e-07 4.1e+00 -1.0e-13 0.0000 p
% 8.302 2.8e+01 2.4e-07 4.1e+00 -7.2e-14 0.0000 p
% 4.201 2.8e+01 2.4e-07 4.1e+00 -2.5e-25 0.0000 p
% 2.151 2.7e+01 2.4e-07 4.1e+00 -3.8e-14 0.0000 p
% 1.125 2.6e+01 2.4e-07 4.1e+00 -5.4e-24 0.0000 p
% 0.613 2.3e+01 -3.7e+01# 4.1e+00 0.0e+00 0.0000 f
% 0.869 2.6e+01 -3.7e+02# 4.1e+00 0.0e+00 0.0000 f
% 0.997 2.6e+01 -1.1e+04# 4.1e+00 0.0e+00 0.0000 f
% 1.061 2.6e+01 2.4e-07 4.1e+00 0.0e+00 0.0000 p
% 1.029 2.6e+01 2.4e-07 4.1e+00 0.0e+00 0.0000 p
% 1.013 2.6e+01 2.4e-07 4.1e+00 0.0e+00 0.0000 p
% 1.005 2.6e+01 2.4e-07 4.1e+00 0.0e+00 0.0000 p
% 1.001 2.6e+01 -3.1e+04# 4.1e+00 -3.8e-14 0.0000 f
% 1.003 2.6e+01 -2.8e+05# 4.1e+00 0.0e+00 0.0000 f
% 1.004 2.6e+01 2.4e-07 4.1e+00 -5.8e-24 0.0000 p
% 1.004 2.6e+01 2.4e-07 4.1e+00 0.0e+00 0.0000 p
% gamma X>=0 Y>=0 rho(XY)<1 p/f
% 1.795e+01 1.4e-07 0.0e+00 1.481e-16 p
% 2.406e+00 1.4e-07 0.0e+00 3.604e-15 p
% 8.806e-01 -3.1e+02 # -1.4e-16 7.370e-19 f
% 1.456e+00 1.4e-07 0.0e+00 1.499e-18 p
% 1.132e+00 1.4e-07 0.0e+00 8.587e-15 p
% 9.985e-01 1.4e-07 0.0e+00 2.331e-13 p
% 9.377e-01 -7.7e+02 # -6.6e-17 3.744e-14 f
% 9.676e-01 -2.0e+03 # -5.7e-17 1.046e-13 f
% 9.829e-01 -6.6e+03 # -1.1e-16 2.949e-13 f
% 9.907e-01 1.4e-07 0.0e+00 2.374e-19 p
% 9.868e-01 -1.6e+04 # -6.4e-17 5.331e-14 f
% 9.887e-01 -5.1e+04 # -1.5e-17 2.703e-19 f
% 9.897e-01 1.4e-07 0.0e+00 1.583e-11 p
% Limiting gains...
% 9.897e-01 1.5e-07 0.0e+00 1.183e-12 p
% 9.897e-01 6.9e-07 0.0e+00 1.365e-12 p
% Gamma value achieved: 1.0036
% Best performance (actual): 0.9897
% #+end_example
% 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]].
@@ -108,8 +137,10 @@ H1 = 1 - H2;
figure;
tiledlayout(3, 1, 'TileSpacing', 'None', 'Padding', 'None');
ax1 = subplot(2,1,1);
% Magnitude
ax1 = nexttile([2, 1]);
hold on;
set(gca,'ColorOrderIndex',1)
plot(freqs, 1./abs(squeeze(freqresp(W1, freqs, 'Hz'))), '--', 'DisplayName', '$w_1$');
@@ -125,10 +156,13 @@ hold off;
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
ylabel('Magnitude');
set(gca, 'XTickLabel',[]);
ylim([5e-4, 20]);
legend('location', 'northeast');
ylim([1e-4, 20]);
yticks([1e-4, 1e-3, 1e-2, 1e-1, 1, 1e1]);
leg = legend('location', 'southeast', 'FontSize', 8, 'NumColumns', 2);
leg.ItemTokenSize(1) = 18;
ax2 = subplot(2,1,2);
% Phase
ax2 = nexttile;
hold on;
set(gca,'ColorOrderIndex',1)
plot(freqs, 180/pi*phase(squeeze(freqresp(H1, freqs, 'Hz'))), '-');
@@ -137,8 +171,7 @@ plot(freqs, 180/pi*phase(squeeze(freqresp(H2, freqs, 'Hz'))), '-');
hold off;
xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
set(gca, 'XScale', 'log');
yticks([-360:90:360]);
yticks([-180:90:180]);
linkaxes([ax1,ax2],'x');
xlim([freqs(1), freqs(end)]);
xticks([0.1, 1, 10, 100, 1000]);

23
matlab/matlab/three_comp_filters.m
View File

@@ -6,6 +6,8 @@ s = zpk('s');
freqs = logspace(-2, 4, 1000);
addpath('./src');
% Weights
% First we define the weights.
@@ -28,9 +30,9 @@ plot(freqs, 1./abs(squeeze(freqresp(W3, freqs, 'Hz'))), '--', 'DisplayName', '$|
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
xlabel('Frequency [Hz]'); ylabel('Magnitude');
hold off;
xlim([freqs(1), freqs(end)]);
xticks([0.01, 0.1, 1, 10, 100, 1000]);
legend('location', 'northeast');
xlim([freqs(1), freqs(end)]); ylim([2e-4, 1.3e1])
leg = legend('location', 'northeast', 'FontSize', 8);
leg.ItemTokenSize(1) = 18;
% H-Infinity Synthesis
% Then we create the generalized plant =P=.
@@ -54,8 +56,10 @@ H3 = tf(H(2));
H1 = 1 - H2 - H3;
figure;
tiledlayout(3, 1, 'TileSpacing', 'None', 'Padding', 'None');
ax1 = subplot(2,1,1);
% Magnitude
ax1 = nexttile([2, 1]);
hold on;
set(gca,'ColorOrderIndex',1)
plot(freqs, 1./abs(squeeze(freqresp(W1, freqs, 'Hz'))), '--', 'DisplayName', '$|W_1|^{-1}$');
@@ -74,10 +78,12 @@ hold off;
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
ylabel('Magnitude');
set(gca, 'XTickLabel',[]);
ylim([5e-4, 20]);
legend('location', 'northeast');
ylim([1e-4, 20]);
leg = legend('location', 'northeast', 'FontSize', 8, 'NumColumns', 2);
leg.ItemTokenSize(1) = 18;
ax2 = subplot(2,1,2);
% Phase
ax2 = nexttile;
hold on;
set(gca,'ColorOrderIndex',1)
plot(freqs, 180/pi*phase(squeeze(freqresp(H1, freqs, 'Hz'))));
@@ -88,8 +94,7 @@ 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]);
yticks([-360:90:360]); ylim([-270, 270]);
linkaxes([ax1,ax2],'x');
xlim([freqs(1), freqs(end)]);
xticks([0.1, 1, 10, 100, 1000]);