--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/m-toolbox/classes/+utils/@math/pfallpsyms2.m	Wed Nov 23 19:22:13 2011 +0100
@@ -0,0 +1,106 @@
+% PFALLPSYMS2 all pass filtering to stabilize TF poles and zeros.
+%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% DESCRIPTION:
+% 
+%     All pass filtering in order to stabilize transfer function poles and
+%     zeros. It inputs a partial fraction expanded discrete model and
+%     outputs a pole-zero minimum phase system
+% 
+% CALL:
+% 
+%     resp= pfallpsyms2(ip,mresp,f,fs)
+% 
+% INPUTS:
+% 
+%     ip: is a struct with fields named poles
+%     mresp: is a vector with functions response
+%     f: is the frequancies vector in (Hz)
+%     fs: is the sampling frequency in (Hz)
+%     
+% OUTPUTS:
+% 
+%     resp: is the stable functions frequency response
+%
+% NOTE:
+%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% VERSION: $Id: pfallpz.m,v 1.6 2009/06/10 15:47:00 luigi Exp $
+% 
+% HISTORY:     12-09-2008 L Ferraioli
+%                 Creation
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+function varargout = pfallpsyms2(ip,mresp,f)
+
+  [a,b] = size(ip);
+  if a<b
+    ip = ip.'; % reshape as a column vector
+  end
+
+  [a,b] = size(f);
+  if a<b
+    f = f.'; % reshape as a column vector
+  end
+  
+
+  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+  Nb = numel(ip);
+  for nn = 1:Nb
+
+    p = ip(nn).poles;
+
+    % stabilizing poles
+    sp = sym(p);
+    unst = real(p) > 0;
+    sp(unst) = conj(sp(unst));
+
+    pp = sym(p(unst));
+    psp = sp(unst);
+    syms s
+    allpsym = 1;
+    for jj=1:numel(psp)
+      allpsym = allpsym.*((s-pp(jj))./(s+psp(jj)));
+    end
+    funcell{nn} = allpsym;
+
+  end
+  
+  fullallprsp = 1;
+  
+  for nn = 1:Nb
+    symexpr = funcell{nn};
+    nterm = subs(symexpr,s,(1i*2*pi).*f);
+    % willing to work with columns
+    if size(nterm,2)>1
+      nterm = nterm.';
+    end
+    
+    fullallprsp = fullallprsp.*nterm;
+    
+  end
+  
+%   rallp = real(fullallprsp);
+%   iallp = imag(fullallprsp);
+%   ang = atan(iallp./rallp);
+  
+  for kk=1:Nb   
+    sresp(:,kk) = mresp(:,kk).*fullallprsp;
+%     sresp(:,kk) = mresp(:,kk).*(cos(ang)+1i.*sin(ang));
+  end
+  
+  for kk=1:Nb   
+    resp(:,kk) = double(sresp(:,kk));
+  end
+  
+  
+  % output
+  if nargout == 1
+    varargout{1} = resp;
+  else
+    error('Too many output arguments!')
+  end
+end
+
+ 
+
+