tp/scripts/partial_fractions_expansion.m (1913B)
1 clc; clear all; close all; 2 3 format shorteng 4 5 function s = eng(x, prec) 6 7 if nargin < 2 8 prec = 3; 9 end 10 11 if x == 0 12 s = sprintf(sprintf("%%.%df", prec), 0); 13 return; 14 end 15 16 e = 3*floor(log10(abs(x))/3); 17 m = x/10^e; 18 19 fmt = sprintf("%%.%dfe%%+03d", prec); 20 s = sprintf(fmt, m, e); 21 22 end 23 24 function print_cosine_form(pole, residue_value) 25 26 sigma = real(pole); 27 omega = imag(pole); 28 29 amplitude = 2*abs(residue_value); 30 phase = angle(residue_value); 31 32 fprintf("\nComplex pair\n"); 33 fprintf("------------\n"); 34 fprintf("p = %s + %sj\n", ... 35 eng(sigma), eng(omega)); 36 37 fprintf("k = %s + %sj\n", ... 38 eng(real(residue_value)), ... 39 eng(imag(residue_value))); 40 41 fprintf("\nTime domain\n"); 42 fprintf("-----------\n"); 43 fprintf("%s * exp(%st) * cos(%st %+0.3f)\n\n", ... 44 eng(amplitude), ... 45 eng(sigma), ... 46 eng(omega), ... 47 phase); 48 end 49 50 % Respuesta al escalón 51 52 fprintf("------------------\n"); 53 fprintf("Heaviside response\n"); 54 fprintf("------------------\n"); 55 56 num = [1 250]; 57 den1 = [1 500]; 58 den2 = [1 4000^2/9000 4000^2]; 59 den3 = [1 0]; 60 den = conv(conv(den1, den2), den3); 61 62 [r, p] = residue(num, den) 63 64 k1 = r(3); p1 = p(3); 65 k2 = r(1); p2 = p(1); 66 k3 = r(2); p3 = p(2); 67 k4 = r(4); p4 = p(4); 68 69 print_cosine_form(p2, k2) 70 71 % Respuesta a seno 72 73 w0 = 2*pi*250; % Aproximadamente 250Hz 74 fprintf("\n---------------------------------------\n"); 75 fprintf("Sine function response with f=%.1f Hz\n", w0/(2*pi)); 76 fprintf("---------------------------------------\n"); 77 78 num = w0*[1 250]; 79 den1 = [1 500]; 80 den2 = [1 4000^2/9000 4000^2]; 81 den3 = [1 0 w0^2]; 82 den = conv(conv(den1, den2), den3); 83 84 [r, p] = residue(num, den) 85 86 k1 = r(5); p1 = p(5); 87 k2 = r(3); p2 = p(3); 88 k3 = r(4); p3 = p(4); 89 k4 = r(1); p4 = p(1); 90 k5 = r(2); p5 = p(2); 91 92 print_cosine_form(p2, k2) 93 print_cosine_form(p4, k4) 94