function sys = fit_zpk(f, h, order, weight);
%function sys = fit_zpk(f, h, order, weight);
%
%Convenience wrapper around vector-fitting algorithm.
%Fits the transfer function h(f) with a zpk-system
%of given order. Weighting vector is optional.
%
%B. Swinkels 2008

fprintf('Fitting %ith order transfer function ', order)

if nargin < 4
  weight = ones(size(f));
  %weight=1 ./ f;
end

%convert to row vects
f = f(:).';
h = h(:).';
weight = weight(:).';

s = 2*pi*i*f;

%complex conjugate pairs, logarithmically spaced :
beta = logspace(log10(2*pi*f(1)),log10(2*pi*f(end)),floor(order/2));
poles = [];
for n = 1:length(beta)
    alpha = -beta(n)*1e-2;
    poles = [poles (alpha-i*beta(n)) (alpha+i*beta(n)) ];
end

if mod(order,2) ~= 0 
  poles = [poles, mean(f)] % add one real pole
end

n_iter = 4;

options.relax = 1;      %Use vector fitting with relaxed non-triviality constraint
options.kill = 2;       %Enforce stable poles

options.asymp = 2;      %Include both D, E in fitting    
options.skip_pole = 0;  %Do NOT skip pole identification
options.use_normal = 1; %Use Normal Equations instead of QR decomp.
options.use_sparse = 1; %Use sparse computations (pole identification)
options.cmplx_ss = 0;   %Create complex state space model

options.spy1 = 0;       %No plotting for first stage of vector fitting
options.spy2 = 0;       %Create magnitude plot for fitting of f(s) 
options.logx = 1;       %Use logarithmic abscissa axis
options.logy = 1;       %Use logarithmic ordinate axis 
options.errplot = 1;    %Include deviation in magnitude plot
options.phaseplot = 1;  %Also produce plot of phase angle (in addition to magnitiude)
options.legend = 1;     %Do include legends in plots

%only compute residuals in last iteration
options.skip_res = 1;

for iter = 1:n_iter
  if iter == n_iter
      options.skip_res = 0;
  end
  [SER, poles, rmserr, fit] = vectfit2(h, s, poles, weight, options);
    
  fprintf('.')
end
fprintf('\n')

%convert to zpk
[z, p, k] = ss2zp(full(SER.A), SER.B, SER.C, SER.D);
sys = zpk(z, p, k);