--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/m-toolbox/classes/@pzmodel/pzm2ab.m	Wed Nov 23 19:22:13 2011 +0100
@@ -0,0 +1,209 @@
+% PZM2AB convert pzmodel to IIR filter coefficients using bilinear transform.
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%
+% DESCRIPTION: PZM2AB convert pzmodel to IIR filter coefficients using bilinear
+%              transform.
+%
+% CALL:        [a,b] = pzm2ab(pzm, fs)
+%
+% VERSION:     $Id: pzm2ab.m,v 1.10 2009/08/31 17:11:41 ingo Exp $
+%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+function varargout = pzm2ab(varargin)
+
+  import utils.const.*
+  
+  %%% Collect input variable names
+  in_names = cell(size(varargin));
+  for ii = 1:nargin,in_names{ii} = inputname(ii);end
+
+  %%% Collect all AOs and plists
+  [pzm, pzm_invars, fs] = utils.helper.collect_objects(varargin(:), 'pzmodel', in_names);
+
+  %%% Check inputs
+  if numel(pzm) ~= 1
+    error('### Please use only one PZ-Model.');
+  end
+  if numel(fs) ~= 1 && ~isnumeric(fs)
+    error('### Please define ''fs''.');
+  else
+    fs = fs{1};
+  end
+
+  gain  = pzm.gain;
+  poles = pzm.poles;
+  zeros = pzm.zeros;
+  np = numel(poles);
+  nz = numel(zeros);
+
+  ao = [];
+  bo = [];
+
+  utils.helper.msg(utils.const.msg.OPROC1, 'converting %s', pzm.name)
+
+  % First we should do complex pole/zero pairs
+  cpoles = [];
+  for j=1:np
+    if poles(j).q > 0.5
+      cpoles = [cpoles poles(j)];
+    end
+  end
+  czeros = [];
+  for j=1:nz
+    if zeros(j).q > 0.5
+      czeros = [czeros zeros(j)];
+    end
+  end
+
+  czi = 1;
+  for j=1:length(cpoles)
+    if czi <= length(czeros)
+      % we have a pair
+      p = cpoles(j);
+      z = czeros(czi);
+
+      [ai,bi] = cpolezero(p.f, p.q, z.f, z.q, fs);
+      if ~isempty(ao)
+        [ao,bo] = pzmodel.abcascade(ao,bo,ai,bi);
+      else
+        ao = ai;
+        bo = bi;
+      end
+
+      % increment zero counter
+      czi = czi + 1;
+    end
+  end
+
+  if length(cpoles) > length(czeros)
+    % do remaining cpoles
+    for j=length(czeros)+1:length(cpoles)
+      utils.helper.msg(msg.OPROC2, 'computing complex pole');
+      [ai,bi] = cp2iir(cpoles(j), fs);
+      if ~isempty(ao)
+        [ao,bo] = pzmodel.abcascade(ao,bo,ai,bi);
+      else
+        ao = ai;
+        bo = bi;
+      end
+    end
+  else
+    % do remaining czeros
+    for j=length(cpoles)+1:length(czeros)
+      utils.helper.msg(msg.OPROC2, 'computing complex zero');
+      [ai,bi] = cz2iir(czeros(j), fs);
+      if ~isempty(ao)
+        [ao,bo] = pzmodel.abcascade(ao,bo,ai,bi);
+      else
+        ao = ai;
+        bo = bi;
+      end
+    end
+  end
+
+  % Now do the real poles and zeros
+  for j=1:np
+    pole = poles(j);
+    if isnan(pole.q) || pole.q < 0.5
+      utils.helper.msg(msg.OPROC2, 'computing real pole');
+      [ai,bi] = rp2iir(pole, fs);
+      if ~isempty(ao)
+        [ao,bo] = pzmodel.abcascade(ao,bo,ai,bi);
+      else
+        ao = ai;
+        bo = bi;
+      end
+    end
+  end
+
+  for j=1:nz
+    zero = zeros(j);
+    if isnan(zero.q) || zero.q < 0.5
+      utils.helper.msg(msg.OPROC2, 'computing real zero');
+      [ai,bi] = rz2iir(zero, fs);
+      if ~isempty(ao)
+        [ao,bo] = pzmodel.abcascade(ao,bo,ai,bi);
+      else
+        ao = ai;
+        bo = bi;
+      end
+    end
+  end
+
+  ao = ao.*gain;
+
+  varargout{1} = ao;
+  varargout{2} = bo;
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%                               Local Functions                               %
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%
+% FUNCTION:    cpolezero
+%
+% DESCRIPTION: Return IIR filter coefficients for a complex pole and
+%              complex zero designed using the bilinear transform.
+%
+% CALL:        [a,b] = cpolezero(pf, pq, zf, zq, fs)
+%
+% HISTORY:     18-02-2007 M Hewitson
+%                Creation.
+%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+function [a,b] = cpolezero(pf, pq, zf, zq, fs)
+
+  import utils.const.*
+  utils.helper.msg(msg.OPROC1, 'computing complex pole/zero pair');
+
+  wp = pf*2*pi;
+  wp2 = wp^2;
+  wz = zf*2*pi;
+  wz2 = wz^2;
+
+  k = 4*fs*fs + 2*wp*fs/pq + wp2;
+
+  a(1) = (4*fs*fs + 2*wz*fs/zq + wz2)/k;
+  a(2) = (2*wz2 - 8*fs*fs)/k;
+  a(3) = (4*fs*fs - 2*wz*fs/zq + wz2)/k;
+  b(1) = 1;
+  b(2) = (2*wp2 - 8*fs*fs)/k;
+  b(3) = (wp2 + 4*fs*fs - 2*wp*fs/pq)/k;
+
+  % normalise dc gain to 1
+  g = iirdcgain(a,b);
+  a = a / g;
+end
+
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%
+% FUNCTION:    iirdcgain
+%
+% DESCRIPTION: Work out the DC gain of an IIR filter given the coefficients.
+%
+% CALL:        g = iirdcgain(a,b)
+%
+% INPUTS:      a - numerator coefficients
+%              b - denominator coefficients
+%
+% OUTPUTS      g - gain
+%
+% HISTORY:     03-07-2002 M Hewitson
+%                Creation.
+%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+function g = iirdcgain(a,b)
+  suma = sum(a);
+  if(length(b)>1)
+    sumb = sum(b);
+    g = suma / sumb;
+  else
+    g = suma;
+  end
+end