zaza

%% first (calcs)
clear

A = 7.2058; % We already found it

% Make mechana (the same for all)
mechane_a = [1 -0.5];
mechane_b = [1 -4 .* exp(1i .* pi ./ 8)];
mechane_c = [1 -4 .* exp(- 1i .* pi ./ 8)];

mechane = poly([0.5 4 .* exp(1i .* pi ./ 8) 4 .* exp(- 1i .* pi ./ 8)]);

% Make mone - first, second, third, ...
mone_a = [1 -0.3];
mone_b = [1 -0.95 .* exp(1i .* pi ./ 3)];
mone_c = [1 -0.95 .* exp(- 1i .* pi ./ 3)];

H1mone = A .* poly([0.3 0.95 .* exp(1i .* pi ./ 3) 0.95 .* exp(- 1i .* pi ./ 3)]);

mone_a = [-0.3 1];
mone_b = [1 -0.95 .* exp(1i .* pi ./ 3)];
mone_c = [1 -0.95 .* exp(- 1i .* pi ./ 3)];

H2mone = A .* 2 .* poly([inv(0.3) 0.95 .* exp(1i .* pi ./ 3) 0.95 .* exp(- 1i .* pi ./ 3)]);

mone_a = [1 -0.3];
mone_b = [-0.95 .* exp(1i .* pi ./ 3) 1];
mone_c = [-0.95 .* exp(- 1i .* pi ./ 3) 1];

H3mone = A .* (1.5 .^ 2) .* poly([0.3 inv(0.95 .* exp(1i .* pi ./ 3)) inv(0.95 .* exp(- 1i .* pi ./ 3))]);

mone_a = [-0.3 1];
mone_b = [-0.95 .* exp(1i .* pi ./ 3) 1];
mone_c = [-0.95 .* exp(- 1i .* pi ./ 3) 1];

H4mone = A .* 2 .* (1.5 .^ 2) .* poly([inv(0.3) inv(0.95 .* exp(1i .* pi ./ 3)) inv(0.95 .* exp(- 1i .* pi ./ 3))]);

H1 = freqz(mechane,H1mone);
H2 = freqz(mechane,H2mone);
H3 = freqz(mechane,H3mone);
H4 = freqz(mechane,H4mone);

%% first plot
w = linspace(0, pi , length(H1));

subplot(4,2,1)
plot(w,abs(H1))
xlabel('\omega (rad/s)')
ylabel('|H_1(e^{j\omega})|')
title('|H1(e^{j\omega})|');

subplot(4,2,2)
zplane(H1mone,mechane);
%---------------------------------
subplot(4,2,3)
plot(w,abs(H2))
xlabel('\omega (rad/s)')
ylabel('|H_2(e^{j\omega})|')
title('|H2(e^{j\omega})|');

subplot(4,2,4)
zplane(H2mone,mechane);
%----------------------------------
subplot(4,2,5)
plot(w,abs(H3))
xlabel('\omega (rad/s)')
ylabel('|H_3(e^{j\omega})|')
title('|H3(e^{j\omega})|');

subplot(4,2,6)
zplane(H3mone,mechane);

%---------------------------------
subplot(4,2,7)
plot(w,abs(H4))
xlabel('\omega (rad/s)')
ylabel('|H_4(e^{j\omega})|')
title('|H4(e^{j\omega})|');

subplot(4,2,8)
zplane(H4mone,mechane);

%---------------------------------

%% second plot
subplot(2,2,1)
hold on
plot(w, angle(H1),'r');
plot(w, angle(H2),'g');
plot(w, angle(H3),'b');
plot(w, angle(H4),'c');
hold off

xlabel('\omega (rad/s)')
ylabel('Angle (rad)')
title('Phase Response');
legend('H1','H2','H3','H4')

%-------------------------------
subplot(2,2,2)

hold on
plot(w, grpdelay(H1mone,mechane),'r');
plot(w, grpdelay(H2mone,mechane),'g');
plot(w, grpdelay(H3mone,mechane),'b');
plot(w, grpdelay(H4mone,mechane),'c');
hold off

title('Group Delay');
xlabel('\omega (rad/s)')
ylabel('t grd')
legend('H1','H2','H3','H4')

%-------------------------------
subplot(2,2,3)

[H1d,W1] = impz(H1mone,mechane);
[H2d,W2] = impz(H2mone,mechane);
[H3d,W3] = impz(H3mone,mechane);
[H4d,W4] = impz(H4mone,mechane);
hold on
plot(W1, H1d,'r');
plot(W2, H2d,'g');
plot(W3, H3d,'b');
plot(W4, H4d,'c');
hold off
xlabel('\omega (rad/s)')
title('Response to delta');
legend('H1','H2','H3','H4')

%-------------------------------
subplot(2,2,4)

Z1 = zeros(1,50);
Z2 = zeros(1,50);
Z3 = zeros(1,50);
Z4 = zeros(1,50);

for m = 1:51
 Z1(m)=sum((H1d(W1<=m)).^2);
 Z2(m)=sum((H2d(W2<=m)).^2);
 Z3(m)=sum((H3d(W3<=m)).^2);
 Z4(m)=sum((H4d(W4<=m)).^2);
end

hold on
plot(0:1:50, Z1, 'r');
plot(0:1:50, Z2, 'g');
plot(0:1:50, Z3, 'b');
plot(0:1:50, Z4, 'c');
hold off

title('Energy');
legend('H1','H2','H3','H4')