vectfit3.zip_matlab vectfit_矢量匹配_矢量拟合

function [SER,poles,rmserr,fit,opts]=vectfit3(f,s,poles,weight,opts); % %function [SER,poles,rmserr,fit,opts]=vectfit3(f,s,poles,weight,opts) %function [SER,poles,rmserr,fit]=vectfit3(f,s,poles,weight,opts) %function [SER,poles,rmserr,fit]=vectfit3(f,s,poles,weight) % % =========================================================== % = Fast Relaxed Vector Fitting = % = Version 1.0 = % = Last revised: 08.08.2008 = % = Written by: Bjorn Gustavsen = % = SINTEF Energy Research, N-7465 Trondheim, NORWAY = % = bjorn.gustavsen@sintef.no = % = http://www.energy.sintef.no/Produkt/VECTFIT/index.asp = % = Note: RESTRICTED to NON-COMMERCIAL use = % =========================================================== % % PURPOSE : Approximate f(s) with a state-space model % % f(s)=C*(s*I-A)^(-1)*B +D +s*E % % where f(s) is a singe element or a vector of elements. % When f(s) is a vector, all elements become fitted with a common % pole set. % % INPUT : % % f(s) : function (vector) to be fitted. % dimension : (Nc,Ns) % Nc : number of elements in vector % Ns : number of frequency samples % % s : vector of frequency points [rad/sec] % dimension : (1,Ns) % % poles : vector of initial poles [rad/sec] % dimension : (1,N) % % weight: the rows in the system matrix are weighted using this array. Can be used % for achieving higher accuracy at desired frequency samples. % If no weighting is desired, use unitary weights: weight=ones(1,Ns). % % Two dimensions are allowed: % dimension : (1,Ns) --> Common weighting for all vector elements. % dimension : (Nc,Ns)--> Individual weighting for vector elements. % % opts.relax==1 --> Use relaxed nontriviality constraint % opts.relax==0 --> Use nontriviality constraint of "standard" vector fitting % % opts.stable=0 --> unstable poles are kept unchanged % opts.stable=1 --> unstable poles are made stable by 'flipping' them % into the left half-plane % % % opts.asymp=1 --> Fitting with D=0, E=0 % opts.asymp=2 --> Fitting with D~=0, E=0 % opts.asymp=3 --> Fitting with D~=0, E~=0 % % opts.spy1=1 --> Plotting, after pole identification (A) % figure(3): magnitude functions % cyan trace : (sigma*f)fit % red trace : (sigma)fit % green trace : f*(sigma)fit - (sigma*f)fit % % opts.spy2=1 --> Plotting, after residue identification (C,D,E) % figure(1): magnitude functions % figure(2): phase angles % % opts.logx=1 --> Plotting using logarithmic absissa axis % % opts.logy=1 --> Plotting using logarithmic ordinate axis % % opts.errplot=1 --> Include deviation in magnitude plot % % opts.phaseplot=1 -->Show plot also for phase angle % % % opts.skip_pole=1 --> The pole identification part is skipped, i.e (C,D,E) % are identified using the initial poles (A) as final poles. % % opts.skip_res =1 --> The residue identification part is skipped, i.e. only the % poles (A) are identified while C,D,E are returned as zero. % % % opts.cmplx_ss =1 -->The returned state-space model has real and complex conjugate % parameters. Output variable A is diagonal (and sparse). % =0 -->The returned state-space model has real parameters only. % Output variable A is square with 2x2 blocks (and sparse). % % OUTPUT : % % fit(s) = C*(s*I-(A)^(-1)*B +D +s.*E % % SER.A(N,N) : A-matrix (sparse). If cmplx_ss==1: Diagonal and complex. % Otherwise, square and real with 2x2 blocks. % % SER.B(N,1) : B-matrix. If cmplx_ss=1: Column of 1's. % If cmplx_ss=0: contains 0's, 1's and 2's) % SER.C(Nc,N) : C-matrix. If cmplx_ss=1: complex % If cmplx_ss=0: real-only % SERD.D(Nc,1) : constant term (real). Is non-zero if asymp=2 or 3. % SERE.E(Nc,1) : proportional term (real). Is non-zero if asymp=3. % % poles(1,N) : new poles % % rmserr(1) : root-mean-square error of approximation for f(s). % (0 is returned if skip_res==1) % fit(Nc,Ns): Rational approximation at samples. (0 is returned if % skip_res==1). % % % APPROACH: % The identification is done using the pole relocating method known as Vector Fitting [1], % with relaxed non-triviality constraint for faster convergence and smaller fitting errors [2], % and utilization of matrix structure for fast solution of the pole identifion step [3]. % %******************************************************************************** % NOTE: The use of this program is limited to NON-COMMERCIAL usage only. % If the program code (or a modified version) is used in a scientific work, % then reference should be made to the following: % % [1] B. Gustavsen and A. Semlyen, "Rational approximation of frequency % domain responses by Vector Fitting", IEEE Trans. Power Delivery, % vol. 14, no. 3, pp. 1052-1061, July 1999. % % [2] B. Gustavsen, "Improving the pole relocating properties of vector % fitting", IEEE Trans. Power Delivery, vol. 21, no. 3, pp. 1587-1592, % July 2006. % % [3] D. Deschrijver, M. Mrozowski, T. Dhaene, and D. De Zutter, % "Macromodeling of Multiport Systems Using a Fast Implementation of % the Vector Fitting Method", IEEE Microwave and Wireless Components % Letters, vol. 18, no. 6, pp. 383-385, June 2008. %******************************************************************************** % This example script is part of the vector fitting package (VFIT3.zip) % Last revised: 08.08.2008. % Created by: Bjorn Gustavsen. % % def.relax=1; %Use vector fitting with relaxed non-triviality constraint def.stable=1; %Enforce stable poles def.asymp=2; %Include only D in fitting (not E) def.skip_pole=0; %Do NOT skip pole identification def.skip_res=0; %Do NOT skip identification of residues (C,D,E) def.cmplx_ss=1; %Create complex state space model def.spy1=0; %No plotting for first stage of vector fitting def.spy2=1; %Create magnitude plot for fitting of f(s) def.logx=1; %Use logarithmic abscissa axis def.logy=1; %Use logarithmic ordinate axis def.errplot=1; %Include deviation in magnitude plot def.phaseplot=0; %exclude plot of phase angle (in addition to magnitiude) def.legend=1; %Do include legends in plots if nargin<5 opts=def; else %Merge default values into opts A=fieldnames(def); for m=1:length(A) if ~isfield(opts,A(m)) dum=char(A(m)); dum2=getfield(def,dum); opts=setfield(opts,dum,dum2); end end end %Tolerances used by relaxed version of vector fitting TOLlow=1e-18; TOLhigh=1e18; [a,b]=size(poles); if s(1)==0 && a==1 if poles(1)==0 && poles(2)~=0 poles(1)=-1; elseif poles(2)==0 && poles(1)~=0 poles(2)=-1; elseif poles(1)==0 && poles(2)==0 poles(1)=-1+i*10; poles(2)=-1-i*10; end end if (opts.relax~=0) && (opts.relax)~=1 disp([' ERROR in vectfit3.m: ==> Illegal value for opts.relax: ' num2str(opts.asymp)]),return end if (opts.asymp~=1) && (opts.asymp)~=2 && (opts.asymp)~=3 disp([' ERROR in vectfit3.m: ==> Illegal value for opts.asymp: ' num2str(opts.asymp)]),return end if (opts.stable~=0) && (opts.stable~=1) disp([' ERROR in vectfit3.m: ==> Ill