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

Add functions to generate weights and CF

Browse Source
This commit is contained in:
Thomas Dehaeze 2021-04-29 15:03:59 +02:00
parent 3fb03c1d29
commit ea35194620
3 changed files with 915 additions and 546 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

304
matlab/index.org
View File

@ -76,6 +76,15 @@ This document is divided into several sections:
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
We here synthesize two complementary filters using the $\mathcal{H}_\infty$ synthesis.
The goal is to specify upper bounds on the norms of the two complementary filters $H_1(s)$ and $H_2(s)$ while ensuring their complementary property ($H_1(s) + H_2(s) = 1$).
@ -321,6 +330,15 @@ The obtained complementary filters are shown on figure [[fig:hinf_filters_result
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
We want:
\begin{align*}
@ -521,6 +539,15 @@ The FIR complementary filters designed in cite:hua05_low_ligo are of order 512.
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
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}$
@ -951,7 +978,7 @@ Let's now compare the FIR filters designed in cite:hua05_low_ligo and the one ob
set(gca,'ColorOrderIndex',1);
plot(w, 180/pi*angle(H), '--');
set(gca,'ColorOrderIndex',2);
set(gca,' H2ColorOrderIndex',2);
plot(w, 180/pi*angle(1-H), '--');
set(gca, 'XScale', 'log');
xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
@ -972,6 +999,7 @@ 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)
<<matlab-dir>>
@ -981,6 +1009,15 @@ Let's now compare the FIR filters designed in cite:hua05_low_ligo and the one ob
<<matlab-init>>
#+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
In order to synthesize the complementary filter using the proposed method, we can use two alternative generalized plant as shown in Figures [[fig:h_infinity_arch_1]] and [[fig:h_infinity_arch_2]].
@ -1273,7 +1310,124 @@ 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)
<<matlab-dir>>
@ -1283,6 +1437,15 @@ And compare $H_1(s)$ with $1 - H_2(s)$ and $H_2(s)$ with $1 - H_1(s)$ in Figure
<<matlab-init>>
#+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
Suppose we want $H_2(s)$ to be an high pass filter with a slope of +2 at low frequency (from 0Hz).
@ -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:
<<sec:generateWF>>
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:
<<sec:generateCF>>
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

28
matlab/matlab/src/generateCF.m Normal file
View File

@ -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;

39
matlab/matlab/src/generateWF.m Normal file
View File

@ -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