# zaza

May 4th, 2021
469
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
1. %% first (calcs)
2. clear
3.
4. A = 7.2058; % We already found it
5.
6. % Make mechana (the same for all)
7. mechane_a = [1 -0.5];
8. mechane_b = [1 -4 .* exp(1i .* pi ./ 8)];
9. mechane_c = [1 -4 .* exp(- 1i .* pi ./ 8)];
10.
11. mechane = poly([0.5 4 .* exp(1i .* pi ./ 8) 4 .* exp(- 1i .* pi ./ 8)]);
12.
13. % Make mone - first, second, third, ...
14. mone_a = [1 -0.3];
15. mone_b = [1 -0.95 .* exp(1i .* pi ./ 3)];
16. mone_c = [1 -0.95 .* exp(- 1i .* pi ./ 3)];
17.
18. H1mone = A .* poly([0.3 0.95 .* exp(1i .* pi ./ 3) 0.95 .* exp(- 1i .* pi ./ 3)]);
19.
20. mone_a = [-0.3 1];
21. mone_b = [1 -0.95 .* exp(1i .* pi ./ 3)];
22. mone_c = [1 -0.95 .* exp(- 1i .* pi ./ 3)];
23.
24. H2mone = A .* 2 .* poly([inv(0.3) 0.95 .* exp(1i .* pi ./ 3) 0.95 .* exp(- 1i .* pi ./ 3)]);
25.
26. mone_a = [1 -0.3];
27. mone_b = [-0.95 .* exp(1i .* pi ./ 3) 1];
28. mone_c = [-0.95 .* exp(- 1i .* pi ./ 3) 1];
29.
30. H3mone = A .* (1.5 .^ 2) .* poly([0.3 inv(0.95 .* exp(1i .* pi ./ 3)) inv(0.95 .* exp(- 1i .* pi ./ 3))]);
31.
32. mone_a = [-0.3 1];
33. mone_b = [-0.95 .* exp(1i .* pi ./ 3) 1];
34. mone_c = [-0.95 .* exp(- 1i .* pi ./ 3) 1];
35.
36. H4mone = A .* 2 .* (1.5 .^ 2) .* poly([inv(0.3) inv(0.95 .* exp(1i .* pi ./ 3)) inv(0.95 .* exp(- 1i .* pi ./ 3))]);
37.
38. H1 = freqz(mechane,H1mone);
39. H2 = freqz(mechane,H2mone);
40. H3 = freqz(mechane,H3mone);
41. H4 = freqz(mechane,H4mone);
42.
43. %% first plot
44. w = linspace(0, pi , length(H1));
45.
46. subplot(4,2,1)
47. plot(w,abs(H1))
49. ylabel('|H_1(e^{j\omega})|')
50. title('|H1(e^{j\omega})|');
51.
52. subplot(4,2,2)
53. zplane(H1mone,mechane);
54. %---------------------------------
55. subplot(4,2,3)
56. plot(w,abs(H2))
58. ylabel('|H_2(e^{j\omega})|')
59. title('|H2(e^{j\omega})|');
60.
61. subplot(4,2,4)
62. zplane(H2mone,mechane);
63. %----------------------------------
64. subplot(4,2,5)
65. plot(w,abs(H3))
67. ylabel('|H_3(e^{j\omega})|')
68. title('|H3(e^{j\omega})|');
69.
70. subplot(4,2,6)
71. zplane(H3mone,mechane);
72.
73. %---------------------------------
74. subplot(4,2,7)
75. plot(w,abs(H4))
77. ylabel('|H_4(e^{j\omega})|')
78. title('|H4(e^{j\omega})|');
79.
80. subplot(4,2,8)
81. zplane(H4mone,mechane);
82.
83. %---------------------------------
84.
85. %% second plot
86. subplot(2,2,1)
87. hold on
88. plot(w, angle(H1),'r');
89. plot(w, angle(H2),'g');
90. plot(w, angle(H3),'b');
91. plot(w, angle(H4),'c');
92. hold off
93.
96. title('Phase Response');
97. legend('H1','H2','H3','H4')
98.
99. %-------------------------------
100. subplot(2,2,2)
101.
102. hold on
103. plot(w, grpdelay(H1mone,mechane),'r');
104. plot(w, grpdelay(H2mone,mechane),'g');
105. plot(w, grpdelay(H3mone,mechane),'b');
106. plot(w, grpdelay(H4mone,mechane),'c');
107. hold off
108.
109. title('Group Delay');
111. ylabel('t grd')
112. legend('H1','H2','H3','H4')
113.
114. %-------------------------------
115. subplot(2,2,3)
116.
117. [H1d,W1] = impz(H1mone,mechane);
118. [H2d,W2] = impz(H2mone,mechane);
119. [H3d,W3] = impz(H3mone,mechane);
120. [H4d,W4] = impz(H4mone,mechane);
121. hold on
122. plot(W1, H1d,'r');
123. plot(W2, H2d,'g');
124. plot(W3, H3d,'b');
125. plot(W4, H4d,'c');
126. hold off
128. title('Response to delta');
129. legend('H1','H2','H3','H4')
130.
131. %-------------------------------
132. subplot(2,2,4)
133.
134. Z1 = zeros(1,50);
135. Z2 = zeros(1,50);
136. Z3 = zeros(1,50);
137. Z4 = zeros(1,50);
138.
139. for m = 1:51
140.  Z1(m)=sum((H1d(W1<=m)).^2);
141.  Z2(m)=sum((H2d(W2<=m)).^2);
142.  Z3(m)=sum((H3d(W3<=m)).^2);
143.  Z4(m)=sum((H4d(W4<=m)).^2);
144. end
145.
146. hold on
147. plot(0:1:50, Z1, 'r');
148. plot(0:1:50, Z2, 'g');
149. plot(0:1:50, Z3, 'b');
150. plot(0:1:50, Z4, 'c');
151. hold off
152.
153. title('Energy');
154. legend('H1','H2','H3','H4')
RAW Paste Data