--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/m-toolbox/classes/+utils/@math/pfresp.m	Wed Nov 23 19:22:13 2011 +0100
@@ -0,0 +1,153 @@
+% PFRESP returns frequency response of a partial fraction TF.
+% 
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%
+% DESCRIPTION:
+%
+% Returns frequency response of a partial fraction expanded function
+% (continuous or discrete).
+%
+% Continuous case
+% The expected model is:
+%
+%          r1             rN
+% f(s) = ------ + ... + ------ + d1 + d2*s + ... + dK*s^{K-1}
+%        s - p1         s - p1
+%
+% Discrete case
+% The expected model is:
+%
+%         z*r1           z*rN
+% f(z) = ------ + ... + ------ + d1 + d2*z^{-1} + ... + dK*z^{-(K-1)}
+%        z - p1         z - p1
+%
+% NOTE: The function cannot handle poles multiplicity higher than 1 in
+% z domain. Multiple poles in s-domain are accepted.
+%
+% CALL:   pfr = pfresp(pfparams)
+%
+% INPUT:
+%
+% pfparams is a struct containing input parameters
+%   pfparams.type = 'cont' Assumes a continuous model
+%   pfparams.type = 'disc' Assumes a discrete model
+%   pfparams.freq set the frequencies vector in Hz
+%   pfparams.res - set the vector of residues
+%   pfparams.pol - set the vector of poles
+%   pfparams.pmul - set the vectr flag with poles multiplicity (this option
+%   is used only for continuous models)
+%   pfparams.dterm - set the vector of direct terms
+%   pfparams.fs - set the sampling frequency (Necessary for the discrete
+%   case)
+%
+% OUTPUT:
+%
+% pfr is a struct containing output data and parameters
+%   pfr.type = 'cont' if the model is continuous
+%   pfr.type = 'disc' if the model is discrete
+%   pfr.freq - frequencies vector
+%   pfr.nfreq - normalized frequencies vector (Discrete case)
+%   pfr.angfreq - angular frequencies vector
+%   pfr.resp - frequency response data
+%
+%
+% VERSION: $Id: pfresp.m,v 1.5 2011/02/03 15:09:01 luigi Exp $
+%
+% HISTORY:     12-09-2008 L Ferraioli
+%                 Creation
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+function pfr = pfresp(pfparams)
+  
+  %%% switching between continuous and discrete
+  
+  switch pfparams.type
+    case 'cont'
+      % collecting input parameters
+      f = pfparams.freq;
+      % willing to work with row
+      if size(f,1)>size(f,2)
+        f = f.';
+      end
+      r = pfparams.res;
+      p = pfparams.pol;
+      pmul = pfparams.pmul;
+      d = pfparams.dterm;
+      % willing to work with row
+      if ~isempty(d) && size(d,1)>size(d,2)
+        d = d.';
+      end
+      
+      N = length(p);
+      
+% substituted by faster code 03-Feb-2011      
+%       Nf = length(f);
+%       rsp = zeros(Nf,1);
+%       indx = (0:length(d)-1).';
+%       for ii = 1:Nf
+%         for jj = 1:N
+%           m = pmul(jj);
+%           rsptemp = r(jj)/(1i*2*pi*f(ii)-p(jj))^m;
+%           rsp(ii) = rsp(ii) + rsptemp;
+%         end
+%         % Direct terms response
+%         rsp(ii) = rsp(ii) + sum((((1i*2*pi*f(ii))*ones(length(d),1)).^indx).*d);
+%       end
+
+% new code for a faster response calculation 03-Feb-2011
+      rsp = zeros(size(f));
+      for jj = 1:N
+        m = pmul(jj);
+        rsptemp = r(jj)./(1i*2*pi.*f-p(jj)).^m;
+        rsp = rsp + rsptemp;
+      end
+      % get direct term response
+      if ~isempty(d)
+        Z = ones(numel(d),numel(f));
+        ss = 2.*pi.*1i.*f;
+        for jj=2:size(Z,1)
+          Z(jj,:) = ss.^(jj-1);
+        end
+        rdtemp = d*Z;
+        rsp = rsp + sum(rdtemp,1);
+      end
+      
+      % Output
+      pfr.type = 'cont';
+      pfr.freq = f;
+      pfr.angfreq = 2*pi*f;
+      pfr.resp = rsp;
+      
+    case 'disc'
+      % collecting input parameters
+      f = pfparams.freq;
+      fs = pfparams.fs;
+      r = pfparams.res;
+      p = pfparams.pol;
+      d = pfparams.dterm;
+      
+      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
+      
+      % Output
+      pfr.type = 'disc';
+      pfr.freq = f;
+      pfr.nfreq = fn;
+      pfr.angfreq = 2*pi*f;
+      pfr.resp = rsp;
+  end
+end
\ No newline at end of file