From: <par...@us...> - 2009-08-26 12:27:53
|
Revision: 6142 http://octave.svn.sourceforge.net/octave/?rev=6142&view=rev Author: paramaniac Date: 2009-08-26 12:27:47 +0000 (Wed, 26 Aug 2009) Log Message: ----------- margin: add support for discrete systems Modified Paths: -------------- trunk/octave-forge/main/control/inst/margin.m Modified: trunk/octave-forge/main/control/inst/margin.m =================================================================== --- trunk/octave-forge/main/control/inst/margin.m 2009-08-26 04:12:38 UTC (rev 6141) +++ trunk/octave-forge/main/control/inst/margin.m 2009-08-26 12:27:47 UTC (rev 6142) @@ -73,9 +73,9 @@ ## @end example ## @end deftypefn -## Version: 0.2.2 +## Version: 0.3 -function [gamma, phi, w_gamma, w_phi] = margin (sys, tol) +function [gamma, phi, w_gamma, w_phi] = margin1 (sys, tol) ## check whether arguments are OK if (nargin < 1 || nargin > 2) @@ -101,11 +101,14 @@ ## get transfer function [num, den, Ts] = sys2tf (sys); - if (Ts == 0) # continuous system + if (Ts == 0) # CONTINUOUS SYSTEM + ## create polynomials s -> jw l_num = length (num); l_den = length (den); + num_jw = zeros (1, l_num); + den_jw = zeros (1, l_den); for k = 1 : l_num num_jw(k) = num(k) * i^(l_num - k); @@ -204,10 +207,159 @@ w_phi = NaN; endif - else # discrete system - error ("margin: discrete systems are not supported yet"); + else # DISCRETE SYSTEM + ## create polynomials z -> 1/z + l_num = length (num); + l_den = length (den); + num_inv = zeros (1, l_num); + den_inv = zeros (1, l_den); + + for k = 1 : l_num + num_inv(k) = num(l_num + 1 - k); + endfor + + for k = 1 : l_den + den_inv(k) = den(l_den + 1 - k); + endfor + + z_num = zeros (1, l_num); + z_den = zeros (1, l_den); + z_num(1) = 1; + z_den(1) = 1; + + ## GAIN MARGIN + ## create gm polynomial + poly_1 = conv (conv (num, den_inv), z_num); + poly_2 = conv (conv (num_inv, den), z_den); + + ## make polynomials equally long for subtraction + l_p1 = length (poly_1); + l_p2 = length (poly_2); + l_max = max (l_p1, l_p2); + + poly_eq_1 = zeros (1, l_max); + poly_eq_2 = zeros (1, l_max); + + for k = 1 : l_max + if (l_max - k < l_p1) + poly_eq_1(k) = poly_1(k + l_p1 - l_max); + endif + + if (l_max - k < l_p2) + poly_eq_2(k) = poly_2(k + l_p2 - l_max); + endif + endfor + + ## subtract polynomials + gm_poly = poly_eq_1 - poly_eq_2; + + ## find frequencies z + z = roots (gm_poly); + + ## filter results + idx = find (abs (abs (z) - 1) < tol); # find z with magnitude 1 + + if (length (idx) > 0) # if z with magnitude 1 exist + z_gm = z(idx); + w = log (z_gm) / (i*Ts); # get frequencies w from z + + idx_1 = find (abs (imag (w)) < tol); # find real frequencies + idx_2 = find (real (w) > 0); # find positive frequencies + idx = intersect (idx_1, idx_2); # find frequencies in R+ + + if (length (idx) > 0) # if frequencies in R+ exist + w_gm = real (w(idx)); + + for k = 1 : length (w_gm) + f_resp(k) = polyval (num, exp (i*w_gm(k)*Ts)) / polyval (den, exp (i*w_gm(k)*Ts)); + gm(k) = inv (abs (f_resp(k))); + endfor + + ## find crossings between 0 and -1 + idx_3 = find (real (f_resp) < 0); + idx_4 = find (real (f_resp) >= -1); + idx_5 = intersect (idx_3, idx_4); + + if (length (idx_5) > 0) # if crossings between 0 and -1 exist + gm = gm(idx_5); + w_gm = w_gm(idx_5); + + [gamma, idx_6] = min (gm); + w_gamma = w_gm(idx_6); + else # there are no crossings between 0 and -1 + gamma = Inf; + w_gamma = NaN; + endif + else # there are no frequencies in R+ + gamma = Inf; + w_gamma = NaN; + endif + else # there are no z with magnitude 1 + gamma = Inf; + w_gamma = NaN; + endif + + ## PHASE MARGIN + ## create pm polynomials + poly_1 = conv (conv (num, num_inv), z_den); + poly_2 = conv (conv (den, den_inv), z_num); + + ## make polynomials equally long for subtraction + l_p1 = length (poly_1); + l_p2 = length (poly_2); + l_max = max (l_p1, l_p2); + + poly_eq_1 = zeros (1, l_max); + poly_eq_2 = zeros (1, l_max); + + for k = 1 : l_max + if (l_max - k < l_p1) + poly_eq_1(k) = poly_1(k + l_p1 - l_max); + endif + + if (l_max - k < l_p2) + poly_eq_2(k) = poly_2(k + l_p2 - l_max); + endif + endfor + + ## subtract polynomials + pm_poly = poly_eq_1 - poly_eq_2; + + ## find frequencies z + z = roots (pm_poly); + + ## filter results + idx = find (abs (abs (z) - 1) < tol); # find z with magnitude 1 + + if (length (idx) > 0) # if z with magnitude 1 exist + z_gm = z(idx); + w = log (z_gm) / (i*Ts); # get frequencies w from z + + idx_1 = find (abs (imag (w)) < tol); # find real frequencies + idx_2 = find (real (w) > 0); # find positive frequencies + idx = intersect (idx_1, idx_2); # find frequencies in R+ + + if (length (idx) > 0) # if frequencies in R+ exist + w_pm = real (w(idx)); + + for k = 1 : length (w_pm) + f_resp = polyval (num, exp (i*w_pm(k)*Ts)) / polyval (den, exp (i*w_pm(k)*Ts)); + pm(k) = 180 + arg (f_resp) / pi * 180; + endfor + + [phi, idx_3] = min (pm); + w_phi = w_pm(idx_3); + else # there are no frequencies in R+ + phi = 180; + w_phi = NaN; + endif + else # there are no z with magnitude 1 + phi = 180; + w_phi = NaN; + endif + endif endfunction @@ -215,7 +367,7 @@ %!shared gamma, phi, w_gamma, w_phi, gamma_exp, phi_exp, w_gamma_exp, w_phi_exp %! sys = tf ([24], [1, 6, 11, 6]); -%! [gamma, phi, w_gamma, w_phi] = margin (sys); +%! [gamma, phi, w_gamma, w_phi] = margin1 (sys); %! %! gamma_exp = 2.50000000000000; %! phi_exp = 35.4254199887541; This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. |