--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/m-toolbox/classes/+utils/@math/pfallpsymz2.m	Wed Nov 23 19:22:13 2011 +0100
@@ -0,0 +1,116 @@
+% PFALLPSYMZ2 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= pfallpz2(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 = pfallpsymz2(ip,mresp,f,fs)
+
+  [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
+
+  if isempty(fs)
+    fs = 1;
+  end
+  [a,b] = size(fs);
+  if a ~= b
+    disp(' Fs has to be a number. Only first term will be considered! ')
+    fs = fs(1);
+  end
+  
+
+  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+  Nb = numel(ip);
+  for nn = 1:Nb
+
+    p = ip(nn).poles;
+
+    % stabilizing poles
+    sp = sym(p);
+    unst = abs(p) > 1;
+    sp(unst) = conj(sp(unst));
+
+    pp = sym(p(unst));
+    psp = sp(unst);
+    syms z
+    allpsym = 1;
+    for jj=1:numel(psp)
+      allpsym = allpsym.*((z-pp(jj))./(z*psp(jj)-1));
+    end
+    funcell{nn} = allpsym;
+
+  end
+  
+  fullallprsp = 1;
+  
+  for nn = 1:Nb
+    symexpr = funcell{nn};
+    nterm = subs(symexpr,z,cos((2*pi/fs).*f) + 1i.*sin((2*pi/fs).*f));
+    % willing to work with columns
+    if size(nterm,2)>1
+      nterm = nterm.';
+    end
+    
+    fullallprsp = fullallprsp.*nterm;
+    
+  end
+  
+%   dballprsp = double(fullallprsp);
+%   rallp = real(dballprsp);
+%   iallp = imag(dballprsp);
+%   ang = angle(dballprsp);
+  
+  for kk=1:Nb   
+    sresp(:,kk) = mresp(:,kk).*fullallprsp;
+%     resp(:,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
+
+ 
+
+