--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/m-toolbox/classes/+utils/@math/vdfit.m	Wed Nov 23 19:22:13 2011 +0100
@@ -0,0 +1,674 @@
+% VDFIT: Fit discrete models to frequency responses
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% DESCRIPTION:
+% 
+%     Fits a discrete model to a frequency response using relaxed z-domain
+%     vector fitting algorithm [1 - 3]. The function is able to fit more
+%     than one frequency response per time. In case that more than one
+%     frequency response is passed as input, they are fitted with a set of
+%     common poles. Model functions are expanded in partial fractions:
+% 
+%                r1                  rN
+%     f(z) = ----------- + ... + ----------- + d
+%            1-p1*z^{-1}         1-pN*z^{-1}
+% 
+% CALL:
+% 
+%     [res,poles,dterm,mresp,rdl,mse] = vdfit(y,f,poles,weight,fitin)
+% 
+% INPUTS:
+% 
+%     - y: Is a vector with the frequency response data.
+%     - f: Is the frequency vector in Hz.
+%     - poles: are a set of starting poles.
+%     - weight: are a set of weights used in the fitting procedure.
+%     - fitin: is a struct containing fitting options and parameters. fitin
+%     fields are:
+%       - fitin.stable = 0; fit without forcing poles to be stable.
+%       - fitin.stable = 1; force poles to be stable flipping unstable
+%       poles in the unit circle. z -> 1/z*.
+%       - fitin.dterm = 0; fit with d = 0.
+%       - fitin.dterm = 1; fit with d different from 0.
+%       - fitin.fs = fs; input the sampling frequency in Hz (default value
+%       is 1 Hz).
+%       - fitin.polt = 0; fit without plotting results.
+%       - fitin.plot = 1; plot fit results in loglog scale.
+%       - fitin.plot = 2; plot fit results in semilogx scale.
+%       - fitin.plot = 3; plot fit results in semilogy scale.
+%       - fitin.plot = 4; plot fit results in linear xy scale.
+%       - fitin.ploth = #; a plot handle to define the figure target for
+%       plotting. Default: [1]
+% 
+% OUTPUT:
+% 
+%     - res: vector or residues.
+%     - poles: vector of poles.
+%     - dterm: direct term d.
+%     - mresp: frequency response of the fitted model
+%     - rdl: residuals y - mresp
+%     - mse: normalized men squared error
+% 
+% EXAMPLES:
+% 
+%     - Fit on a single transfer function:
+% 
+%       INPUT
+%       y is a (Nx1) or (1xN) vector
+%       f is a (Nx1) or (1xN) vector
+%       poles is a (Npx1) or (1xNp) vector
+%       weight is a (Nx1) or (1xN) vector
+% 
+%       OUTPUT
+%       res is a (Npx1) vector
+%       poles is a (Npx1) vector
+%       dterm is a constant
+%       mresp is a (Nx1) vector
+%       rdl is a (Nx1) vector
+%       mse is a constant
+% 
+%       - Fit Nt transfer function at the same time:
+% 
+%       INPUT
+%       y is a (NxNt) or (NtxN) vector
+%       f is a (Nx1) or (1xN) vector
+%       poles is a (Npx1) or (1xNp) vector
+%       weight is a (NxNt) or (NtxN) vector
+% 
+%       OUTPUT
+%       res is a (NpxNt) vector
+%       poles is a (Npx1) vector
+%       dterm is a (1xNt) vector
+%       mresp is a (NxNt) vector
+%       rdl is a (NxNt) vector
+%       mse is a (1xNt) vector
+% 
+% REFERENCES:
+% 
+%     [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] Y. S. Mekonnen and J. E. Schutt-Aine, "Fast broadband
+%         macromodeling technique of sampled time/frequency data using
+%         z-domain vector-fitting method", Electronic Components and
+%         Technology Conference, 2008. ECTC 2008. 58th 27-30 May 2008 pp.
+%         1231 - 1235.
+% 
+% 
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% HISTORY:     12-09-2008 L Ferraioli
+%                 Creation
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% VERSION: '$Id: vdfit.m,v 1.12 2009/08/06 09:57:37 miquel Exp $';
+% 
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+function [res,poles,dterm,mresp,rdl,mse] = vdfit(y,f,poles,weight,fitin)
+
+  warning off all
+  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+  % Collecting inputs
+
+  % Default input struct
+  defaultparams = struct('stable',0, 'dterm',0, 'fs',1, 'plot',0, 'ploth',1,...
+    'regout',-1, 'idsamp',1e3, 'idamp', 1e-3);
+
+  names = {'stable','dterm','fs','plot','ploth','regout','idsamp','idamp'};
+
+  % collecting input and default params
+  if ~isempty(fitin)
+    for jj=1:length(names)
+      if isfield(fitin, names(jj)) && ~isempty(fitin.(names{1,jj}))
+       defaultparams.(names{1,jj}) = fitin.(names{1,jj});
+      end
+    end
+  end
+
+  stab = defaultparams.stable; % Enforce pole stability is is 1
+  dt = defaultparams.dterm; % 1 to fit with direct term
+  fs = defaultparams.fs; % sampling frequency
+  plotting = defaultparams.plot; % set to 1 if plotting is required
+  plth = defaultparams.ploth; % set the figure target
+  regout = defaultparams.regout; % set the strategy for complex plane fitting
+  idsamp = defaultparams.idsamp; % number of samples to define impulse response
+  idamp = defaultparams.idamp; % maximum aplitude of impulse response adimtted at idsamp
+  
+  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+  % Inputs in row vectors
+
+  [a,b] = size(y);
+  if a > b % shifting to row
+    y = y.';
+  end
+
+  [a,b] = size(f);
+  if a > b % shifting to row
+    f = f.';
+  end
+
+  [a,b] = size(poles);
+  if a > b % shifting to row
+    poles = poles.';
+  end
+
+  clear w
+  w = weight;
+  [a,b] = size(w);
+  if a > b % shifting to row
+    w = w.';
+  end
+
+  N = length(poles); % Model order
+
+  if dt
+    dl = 1; % Fit with direct term
+  else
+    dl = 0; % Fit without direct term
+  end
+
+  % definition of z
+  z = cos(2.*pi.*f./fs)+1i.*sin(2.*pi.*f./fs);
+
+  Nz = length(z);
+
+  [Nc,Ny] = size(y);
+  if Ny ~= Nz
+    error('### The number of data points is different from the number of frequency points.')
+  end
+
+  %Tolerances used by relaxed version of vector fitting
+  TOLlow=1e-8;
+  TOLhigh=1e8;
+  
+  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+  % Normalizing y
+
+    for nn = 1:Nc
+      y(nn,:) = y(nn,:)./z;
+    end
+
+  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+  % Marking complex and real poles
+
+  % cindex = 1; pole is complex, next conjugate pole is marked with cindex
+  % = 2. cindex = 0; pole is real
+  cindex=zeros(1,N);
+  for m=1:N 
+    if imag(poles(m))~=0  
+      if m==1 
+        cindex(m)=1;
+      else
+        if cindex(m-1)==0 || cindex(m-1)==2
+          cindex(m)=1; cindex(m+1)=2; 
+        else
+          cindex(m)=2;
+        end
+      end 
+    end
+  end
+
+  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+  % Augmented problem
+  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+  % Matrix initialinzation
+  BA = zeros(Nc*Nz+1,1);
+  Ak=zeros(Nz,N+1);
+  AA=zeros(Nz*Nc+1,(N+dl)*Nc+N+1);
+  nf = zeros(1,Nc*(N+dl)+N+1); % Normalization factor
+
+  % Defining Ak
+  for jj = 1:N
+    if cindex(jj) == 1 % conjugate complex couple of poles
+      Ak(:,jj) = 1./(z-poles(jj)) + 1./(z-conj(poles(jj)));
+      Ak(:,jj+1) = 1i./(z-poles(jj)) - 1i./(z-conj(poles(jj)));
+    elseif cindex(jj) == 0 % real pole
+      Ak(:,jj) = 1./(z-poles(jj));
+    end
+  end
+
+  Ak(1:Nz,N+1) = 1;
+
+  % Scaling factor
+  sc = 0;
+  for mm = 1:Nc
+    sc = sc + (norm(w(mm,:).*y(mm,:)))^2;
+  end
+  sc=sqrt(sc)/Nz; 
+
+  for nn = 1:Nc
+
+    wg = w(nn,:).'; % Weights
+
+    ida=(nn-1)*Nz+1;
+    idb=nn*Nz;
+    idc=(nn-1)*(N+dl)+1;
+
+    for mm =1:N+dl % Diagonal blocks
+      AA(ida:idb,idc-1+mm) = wg.*Ak(1:Nz,mm);
+    end
+    for mm =1:N+1 % Last right blocks
+      AA(ida:idb,Nc*(N+dl)+mm) = -wg.*(Ak(1:Nz,mm).*y(nn,1:Nz).');
+    end
+
+  end
+
+  % setting the last row of AA and BA for the relaxion condition
+  for qq = 1:N+1
+    AA(Nc*Nz+1,Nc*(N+dl)+qq) = real(sc*sum(Ak(:,qq)));
+  end
+
+  AA = [real(AA);imag(AA)];
+
+  % Last element of the solution vector
+  BA(Nc*Nz+1) = Nz*sc;
+
+  xBA = real(BA);
+  xxBA = imag(BA);
+  
+  Nrow = Nz*Nc+1;
+  
+  BA = zeros(2*Nrow,1);
+  
+  BA(1:Nrow,1) = xBA;
+  BA(Nrow+1:2*Nrow,1) = xxBA;
+
+  % Normalization factor
+  % nf = zeros(2*N+dl+1,1);
+  for pp = 1:length(AA(1,:))
+    nf(pp) = norm(AA(:,pp),2); % Euclidean norm
+    AA(:,pp) = AA(:,pp)./nf(pp); % Normalization
+  end
+
+  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+  % Solving augmented problem
+
+%   XA = pinv(AA)*BA;
+  % XA = inv((AA.')*AA)*(AA.')*BA;
+
+  % XA = AA.'*AA\AA.'*BA;
+
+  XA = AA\BA;
+
+  XA = XA./nf.'; % renormalization
+  
+  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+  % Checking the tolerance
+  
+  if abs(XA(end))<TOLlow || abs(XA(end))>TOLhigh
+    
+    if XA(end)==0
+      Dnew=1;
+    elseif abs(XA(end))<TOLlow
+      Dnew=sign(XA(end))*TOLlow;
+    elseif abs(XA(end))>TOLhigh
+      Dnew=sign(XA(end))*TOLhigh;
+    end
+
+    for pp = 1:length(AA(1,:))
+      AA(:,pp) = AA(:,pp).*nf(pp); %removing previous scaling
+    end
+
+    ind=length(AA(:,1))/2; %index to additional row related to relaxation
+    
+    AA(ind,:)=[]; % removing relaxation term
+    
+    BA=-Dnew*AA(:,end);  %new right side
+    
+    AA(:,end)=[];
+    
+    nf(end)=[];
+    
+    for pp = 1:length(AA(1,:))
+      nf(pp) = norm(AA(:,pp),2); % Euclidean norm
+      AA(:,pp) = AA(:,pp)./nf(pp); % Normalization
+    end
+          
+    % XA=(AA.'*AA)\(AA.'*BA); % using normal equation
+      
+    XA=AA\BA;
+    
+    XA = XA./nf.'; % renormalization
+    
+    XA=[XA;Dnew];
+    
+  end
+
+  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+  % Finding zeros of sigma
+
+  lsr = XA((N+dl)*Nc+1:(N+dl)*Nc+N,1); % collect the least square results
+
+  D = XA(end); % direct term of sigma
+  
+  CPOLES = diag(poles);
+  B = ones(N,1);
+  C = lsr.';
+  
+  for kk = 1:N
+    if cindex(kk) == 1
+      CPOLES(kk,kk)=real(poles(kk));
+      CPOLES(kk,kk+1)=imag(poles(kk));
+      CPOLES(kk+1,kk)=-1*imag(poles(kk));
+      CPOLES(kk+1,kk+1)=real(poles(kk));
+      B(kk,1) = 2;
+      B(kk+1,1) = 0;
+    end
+  end
+  
+  H = CPOLES-B*C/D;
+  
+  szeros=eig(H);
+  
+  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+  % Exclude a region of the complex plane
+  switch regout
+    case -1
+      % do nothing
+    case 0
+      % do nothing
+    case 1
+      % set the maximum admitted value for stable poles
+      target_pole = (idamp)^(1/idsamp);
+      % get stable poles outside the fixed limit
+      uptgr = ((abs(szeros) > target_pole) & (abs(szeros) <= 1) & (imag(szeros)==0));
+      uptgi = ((abs(szeros) > target_pole) & (abs(szeros) <= 1) & (imag(szeros)~=0));
+      % get unstable polse smaller than minimum value
+      lwtgr = ((abs(szeros) < 1/target_pole) & (abs(szeros) > 1) & (imag(szeros)==0));
+      lwtgi = ((abs(szeros) < 1/target_pole) & (abs(szeros) > 1) & (imag(szeros)~=0));
+      % get the maximum shift needed
+      ushiftr = max(abs(abs(szeros(uptgr))-target_pole));
+      ushifti = max(abs(abs(szeros(uptgi))-target_pole));
+      lshiftr = max(abs(abs(szeros(lwtgr))-1/target_pole));
+      lshifti = max(abs(abs(szeros(lwtgi))-1/target_pole));
+      % shifting inside
+      szeros(uptgr) = (abs(szeros(uptgr))-ushiftr).*sign(szeros(uptgr));
+      szeros(uptgi) = (abs(szeros(uptgi))-ushifti).*exp(1i.*angle(szeros(uptgi)));
+      % shifting outside
+      szeros(lwtgr) = (abs(szeros(lwtgr))+lshiftr).*sign(szeros(lwtgr));
+      szeros(lwtgi) = (abs(szeros(lwtgi))+lshifti).*exp(1i.*angle(szeros(lwtgi)));
+  end
+      
+
+  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+  % Ruling out unstable poles
+
+  % This option force the poles of f to stay inside the unit circle
+  if stab
+    unst = abs(szeros) > 1;
+    szeros(unst) = 1./conj(szeros(unst));
+  end
+  szeros = sort(szeros);
+  N = length(szeros);
+
+  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+  % Separating complex poles from real poles and ordering
+
+  rnpoles = [];
+  inpoles = [];
+  for tt = 1:N
+    if imag(szeros(tt)) == 0
+      % collecting real poles
+      rnpoles = [rnpoles; szeros(tt)];
+    else
+      % collecting complex poles
+      inpoles = [inpoles; szeros(tt)];
+    end
+  end
+
+  % Sorting complex poles in order to have them in the expected order a+jb
+  % and a-jb a>0 b>0
+  inpoles = sort(inpoles);
+  npoles = [rnpoles;inpoles];
+  npoles = npoles - 2.*1i.*imag(npoles);
+
+  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+  % Marking complex and real poles
+
+  cindex=zeros(N,1);
+  for m=1:N 
+    if imag(npoles(m))~=0  
+      if m==1 
+        cindex(m)=1;
+      else
+        if cindex(m-1)==0 || cindex(m-1)==2
+          cindex(m)=1; cindex(m+1)=2; 
+        else
+          cindex(m)=2;
+        end
+      end 
+    end
+  end
+
+  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+  % Direct problem
+  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+  % Matrix initialinzation
+  nB(1:Nz,1:Nc) = real(w.*y).';
+  nB(Nz+1:2*Nz,1:Nc) = imag(w.*y).';
+
+  B = zeros(2*Nz,1);
+  nAD = zeros(Nz,N+dl);
+  AD = zeros(2*Nz,N+dl);
+  Ak = zeros(Nz,N+dl);
+
+  for jj = 1:N
+    if cindex(jj) == 1 % conjugate complex couple of poles
+      Ak(:,jj) = 1./(z-npoles(jj)) + 1./(z-npoles(jj+1));
+      Ak(:,jj+1) = 1i./(z-npoles(jj)) - 1i./(z-npoles(jj+1));
+    elseif cindex(jj) == 0 % real pole
+      Ak(:,jj) = 1./(z-npoles(jj));
+    end
+  end
+
+  if dt
+  %   Ak(1:Nz,N+dl) = ones(Nz,1); % considering the direct term
+    Ak(1:Nz,N+dl) = 1./z;
+  end
+
+  XX = zeros(Nc,N+dl);
+  for nn = 1:Nc
+
+    % Defining AD
+    for m=1:N+dl
+      nAD(1:Nz,m) = w(nn,:).'.*Ak(1:Nz,m); 
+    end 
+
+    B(1:2*Nz,1) = nB(1:2*Nz,nn);
+
+    AD(1:Nz,:) = real(nAD);
+    AD(Nz+1:2*Nz,:) = imag(nAD);
+
+    % Normalization factor
+    nf = zeros(N+dl,1);
+    for pp = 1:N+dl
+      nf(pp,1) = norm(AD(:,pp),2); % Euclidean norm
+      AD(:,pp) = AD(:,pp)./nf(pp,1); % Normalization
+    end
+
+    % Solving direct problem
+
+    % XD = inv((AD.')*AD)*(AD.')*B;
+    % XD = AD.'*AD\AD.'*B;
+%     XD = pinv(AD)*B;
+    XD = AD\B;
+
+    XD = XD./nf; % Renormalization
+    XX(nn,1:N) = XD(1:N).';
+
+  end
+
+  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+  % Final residues and poles of f
+
+  lsr = XX(:,1:N);
+
+  res = zeros(N,Nc);
+  % Real poles have real residues, complex poles have comples residues
+  for nn = 1:Nc
+    for tt = 1:N
+      if cindex(tt) == 1 % conjugate complex couple of poles
+        res(tt,nn) = lsr(nn,tt)+1i*lsr(nn,tt+1);
+        res(tt+1,nn) = lsr(nn,tt)-1i*lsr(nn,tt+1);
+      elseif cindex(tt) == 0 % real pole
+        res(tt,nn) = lsr(nn,tt);
+      end
+    end
+  end
+
+  poles = npoles;
+
+  if dt
+    dterm = XX(:,N+dl).';
+  else
+    dterm = zeros(1,Nc);
+  end
+
+  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+  % Calculating responses and residuals
+
+  mresp = zeros(Nz,Nc);
+  rdl = zeros(Nz,Nc);
+  yr = zeros(Nz,Nc);
+  mse = zeros(1,Nc);
+
+  for nn = 1:Nc
+    % freq resp of the fit model
+    r = res(:,nn);
+    p = poles;
+    d = dterm(:,nn);
+
+    Nf = length(f);
+    N = length(p);
+
+    % Defining normalized frequencies
+    fn = f./fs;
+
+    rsp = zeros(Nf,1);
+    indx = 0:length(d)-1;
+    for ii = 1:Nf
+      for jj = 1:N
+        rsptemp = exp(1i*2*pi*fn(ii))*r(jj)/(exp(1i*2*pi*fn(ii))-p(jj));
+        rsp(ii) = rsp(ii) + rsptemp;
+      end
+      % Direct terms response
+      rsp(ii) = rsp(ii) + sum(((exp((1i*2*pi*f(ii))*ones(length(d),1))).^(-1.*indx)).*d);
+    end
+
+    % Model response
+    mresp(:,nn) = rsp;
+
+    % Residual
+    yr(:,nn) = (y(nn,:).*z).';
+    rdl(:,nn) = yr(:,nn) - rsp;
+    
+    % RMS error
+%     rmse(:,nn) = sqrt(sum((abs(rdl(:,nn)./yr(:,nn)).^2))/(Nf-N));
+    
+    % Chi Square or mean squared error
+    % Note that this error is normalized to the input data in order to
+    % comparable between different sets of data
+    mse(:,nn) = sum((rdl(:,nn)./yr(:,nn)).*conj((rdl(:,nn)./yr(:,nn))))/(Nf-N);
+    
+  end
+
+  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+  % Plotting response
+  nf = f./fs;
+
+  switch plotting
+
+    case 0
+      % No plot
+
+    case 1
+      % LogLog plot for absolute value
+      figure(plth)
+      subplot(2,1,1);
+      p1 = loglog(nf,abs(yr),'k');
+      hold on
+      p2 = loglog(nf,abs(mresp),'r');
+      p3 = loglog(nf,abs(rdl),'b');
+      xlabel('Normalized Frequency [f/fs]')
+      ylabel('Amplitude')
+      legend([p1(1) p2(1) p3(1)],'Original', 'VDFIT','Residual')
+      hold off
+
+      subplot(2,1,2);
+      p4 = semilogx(nf,(180/pi).*unwrap(angle(yr)),'k');
+      hold on
+      p5 = semilogx(nf,(180/pi).*unwrap(angle(mresp)),'r');
+      xlabel('Normalized Frequency [f/fs]')
+      ylabel('Phase [Deg]')
+      legend([p4(1) p5(1)],'Original', 'VDFIT')
+      hold off
+
+    case 2
+      % Semilogx plot for absolute value
+      figure(plth)
+      subplot(2,1,1);
+      p1 = semilogx(nf,abs(yr),'k');
+      hold on
+      p2 = semilogx(nf,abs(mresp),'r');
+      p3 = semilogx(nf,abs(rdl),'b');
+      xlabel('Normalized Frequency [f/fs]')
+      ylabel('Amplitude')
+      legend([p1(1) p2(1) p3(1)],'Original', 'VDFIT','Residual')
+      hold off
+
+      subplot(2,1,2);
+      p4 = semilogx(nf,(180/pi).*unwrap(angle(yr)),'k');
+      hold on
+      p5 = semilogx(nf,(180/pi).*unwrap(angle(mresp)),'r');
+      xlabel('Normalized Frequency [f/fs]')
+      ylabel('Phase [Deg]')
+      legend([p4(1) p5(1)],'Original', 'VDFIT')
+      hold off
+
+    case 3
+      % Semilogy plot for absolute value
+      figure(plth)
+      subplot(2,1,1);
+      p1 = semilogy(nf,abs(yr),'k');
+      hold on
+      p2 = semilogy(nf,abs(mresp),'r');
+      p3 = semilogy(nf,abs(rdl),'b');
+      xlabel('Normalized Frequency [f/fs]')
+      ylabel('Amplitude')
+      legend([p1(1) p2(1) p3(1)],'Original', 'VDFIT','Residual')
+      hold off
+
+      subplot(2,1,2);
+      p4 = semilogy(nf,(180/pi).*unwrap(angle(yr)),'k');
+      hold on
+      p5 = semilogy(nf,(180/pi).*unwrap(angle(mresp)),'r');
+      xlabel('Normalized Frequency [f/fs]')
+      ylabel('Phase [Deg]')
+      legend([p4(1) p5(1)],'Original', 'VDFIT')
+      hold off
+
+    case 4
+      % Linear plot for absolute value
+      figure(plth)
+      subplot(2,1,1);
+      p1 = plot(nf,abs(yr),'k');
+      hold on
+      p2 = plot(nf,abs(mresp),'r');
+      p3 = plot(nf,abs(rdl),'b');
+      xlabel('Normalized Frequency [f/fs]')
+      ylabel('Amplitude')
+      legend([p1(1) p2(1) p3(1)],'Original', 'VDFIT','Residual')
+      hold off
+
+      subplot(2,1,2);
+      p4 = plot(nf,(180/pi).*unwrap(angle(yr)),'k');
+      hold on
+      p5 = plot(nf,(180/pi).*unwrap(angle(mresp)),'r');
+      xlabel('Normalized Frequency [f/fs]')
+      ylabel('Phase [Deg]')
+      legend([p4(1) p5(1)],'Original', 'VDFIT')
+      hold off
+
+  end
\ No newline at end of file