Ronel101

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))
  48. xlabel('\omega (rad/s)')
  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))
  57. xlabel('\omega (rad/s)')
  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))
  66. xlabel('\omega (rad/s)')
  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))
  76. xlabel('\omega (rad/s)')
  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.  
  94. xlabel('\omega (rad/s)')
  95. ylabel('Angle (rad)')
  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');
  110. xlabel('\omega (rad/s)')
  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
  127. xlabel('\omega (rad/s)')
  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

Adblocker detected! Please consider disabling it...

We've detected AdBlock Plus or some other adblocking software preventing Pastebin.com from fully loading.

We don't have any obnoxious sound, or popup ads, we actively block these annoying types of ads!

Please add Pastebin.com to your ad blocker whitelist or disable your adblocking software.

×