
% Script transfer_1.m in which the transfer function of a
% LP filter with three LC sections (6th order) is examined

clear;
% ----- Initializations
Rg = 100;      % Output resistance of the source
L = 30e-6;     % Inductance
C = 20e-12;    % Capacity
R = 10e3;      % Loss Resistance of Capacitance
Rs = 1000;     % Load resistance

w0 = 1 / sqrt (L * C); % Resonance frequency of a sections (rad/s)
f0 = w0 / (2 * pi);    % in Hz

s = tf ('s');  % s as a Laplace variable
% ------- Calculation of the transfer function based on the output
Us = 1;        % Output voltage
I3 = Us * (1 / Rs + 1 / R + C * s);

U2 = I3 * L * s + Us;
I2 = U2 * (1 / R + C * s) + I3;

U1 = I2 * L * s + U2;
I1 = U1 * (1 / R + C * s) + I2;

Ug = I1 * (L * s + Rg) + U1;

Hs = Us / Ug;     % Transfer function (s object)

% step (Hs); impulse (Hs);

[num, den] = tfdata (Hs);   % Cells of the coefficients
b = num{:}; a = den{:};     % Coefficients of the numerator and denominator

f = logspace (round (log10 (f0 / 100)), round (log10 (f0 * 100)), 200);
         % Frequency range with 2 decades left and right
         % of the resonance frequency
[H, w] = freqs (b, a, f*2*pi); % Frequency response

figure (1); clf;
subplot (211),semilogx(w/(2*pi),20*log10(abs(H)),'k-','LineWidth', 1);
title ('Amplitude response'); xlabel ('Hz'); grid on;
ylabel ('dB');

subplot (212),semilogx(w/(2*pi),180*unwrap(angle(H))/pi,'k-','LineWidth',1);
title ('Phase response'); xlabel ('Hz'); grid on;
ylabel ('degree');

% ------- Zero poles
[z, p, k] = zpkdata (Hs);

figure (2); clf;
iopzplot (Hs, 'k-');
grid on;
% ################################################ #########
% ------- 2nd variant
omega = 2*pi*f; % Frequency range (rad / s)

% ------- Calculation of the transfer function based on the output
Us = 1;      % Output voltage
% Us = 1 * exp (j * pi / 3); % Or output voltage
I3 = Us*(1/Rs + 1/R + C*j*omega);

U2 = I3.*(L*j*omega) + Us;
I2 = U2.*(1/R + C*j * omega) + I3;

U1 = I2.*(L*j*omega) + U2;
I1 = U1.*(1/R + C*j*omega) + I2;

Ug = I1.*(L*j*omega + Rg) + U1;

Hs2 = Us./Ug; % Transfer function

figure (3); clf;
subplot(211),semilogx(w/(2*pi),20*log10(abs(Hs2)),'k-','LineWidth',1);
title ('Amplitude response'); xlabel ('Hz'); grid on;
ylabel ('dB');

subplot(212),semilogx(w/(2*pi),180*unwrap(angle(Hs2))/pi,'k-','LineWidth',1);
title ('Phase response'); xlabel ('Hz'); grid on;
ylabel ('degree');

% -------- State model 1
[A1, B1, C1, D1] = tf2ss (b, a);
system1 = ss (A1, B1, C1, D1);

% -------- State model 2
A2 = [-(1/C)*(1/R + 1/Rs),1/C 0 0 0 0;-1/L 0 1/L 0 0 0;
     0 -1/C -1/(R * C) 1/C 0 0;0 0 -1/L 0 1 /L 0;
     0 0 0 -1/C -1/(R * C) 1/C; 0 0 0 0 -1/L -Rg/L];
     
B2 = [0 0 0 0 0 1 / L]';
C2 = [1 0 0 0 0 0];
D2 = 0;
system2 = ss (A2, B2, C2, D2);

% ############
Nt = 1000;
u = ones(Nt, 1);               % Entry jump
delta_t = 1 / 10e8;            % Time step
t = (0:Nt-1)'* delta_t;        % Time vector
x0 = zeros(6,1);               % Initial state vector

[y2, t, x2] = lsim (system2, u, t, x0);

figure (4), clf;
subplot (6,1,1), plot (t, x2 (:, 1), 'k -', 'LineWidth', 1);
title ('State variable x2(:,1) (output voltage us (t))');
subplot (6,1,2), plot (t, x2 (:, 2), 'k -', 'LineWidth', 1);
title ('State variable x2(:,2) (i3 (t))');
subplot (6,1,3), plot (t, x2 (:, 3), 'k -', 'LineWidth', 1);
title ('State variable x2(:,3) (u2 (t))');
subplot (6,1,4), plot (t, x2 (:, 4), 'k -', 'LineWidth', 1);
title ('State variable x2(:,4) (i2 (t))');
subplot (6,1,5), plot (t, x2 (:, 5), 'k -', 'LineWidth', 1);
title ('State variable x2(:,5) (u1 (t))');
subplot (6,1,6), plot (t, x2 (:, 6), 'k -', 'LineWidth', 1);
title ('State variable x2(:,6) (i1 (t))');

[y1, t, x1] = lsim(system1, u, t, x0); % y1 is the output voltage
% equals x1(:, 6)*4.6296e45; % from the multiplication
                     % of the state variables x1 with
                     % C1 = [0 0 0 0 0 4.6296e45]
figure (5), clf;
subplot (6,1,1), plot (t, x1 (:, 1), 'k -', 'LineWidth', 1);
title ('State variable x1(:,1)');
subplot (6,1,2), plot (t, x1 (:, 2), 'k -', 'LineWidth', 1);
title ('State variable x1(:,2)');
subplot (6,1,3), plot (t, x1 (:, 3), 'k -', 'LineWidth', 1);
title ('State variable x1(:,3)');
subplot (6,1,4), plot (t, x1 (:, 4), 'k -', 'LineWidth', 1);
title ('State variable x1(:,4)');
subplot (6,1,5), plot (t, x1 (:, 5), 'k -', 'LineWidth', 1);
title ('State variable x1(:,5)');
subplot (6,1,6), plot (t, y1, 'k -', 'LineWidth', 1);
title ('State variable x1(:,6) * 4.62e45 (output voltage y1 (t))');

% -------- Impulse response with impulse
figure (6), clf;
[y3, timp] = impulse(Hs); % or impulse(system1) or impulse(system2)
plot (timp, y3, 'k -', 'LineWidth', 1);
title ('Impulse response with impulse');

% ------- impulse response with lsim
N = 1000;
u = [zeros(1,50), 1, zeros(1, N-50-1)]';
t = [0:N-1]'* delta_t;
y = lsim(Hs, u, t)/delta_t;

figure (7); clf;
subplot (211), stairs (t, u, 'k -', 'LineWidth', 1);
title ('Input signal (pulse of duration 1/10e8)'); xlabel ('s'); grid on;

subplot (212), plot (t, y, 'k -', 'LineWidth', 1);
title ('Output signal (impulse response)'); xlabel ('s'); grid on;

% -------- Response of the filter to any input signal
% Input signal = white noise
N = 1000;
rng (1238); % Start value for the noise generator
u = randn (1, N);
t = (0:N-1)*delta_t;
x0 = zeros(1,6); % Initial state vector
y = lsim (Hs, u, t, x0);

figure (8); clf;
subplot (211), stairs (t, u, 'k -', 'LineWidth', 1);
title ('Input Signal'); xlabel ('s'); grid on;

subplot (212), plot (t, y, 'k -', 'LineWidth', 1);
title ('Output signal'); xlabel ('s'); grid on;

% Input signal = two sinusoidal signals
fsig1 = 1e6; ampl1 = 1; % Signal in the pass band of the filter
fsig2 = 1e8; ampl2 = 2; % Signal in the stop range of the filter
N = 5000;
t = (0: N-1) * delta_t;

u = ampl1*sin(2*pi*fsig1*t) + ampl1*sin(2*pi*fsig2*t);
y = lsim (Hs, u, t, x0);

figure (9); clf;
subplot (211), stairs (t, u, 'k -', 'LineWidth', 1);
title ('Input Signal'); xlabel ('s'); grid on;

subplot (212), plot (t, y, 'k -', 'LineWidth', 1);
title ('Output signal'); xlabel ('s'); grid on;
