Home > classes > @pzmodel > private > pzm2ab.m

pzm2ab

PURPOSE ^

PZM2AB convert pzmodel to IIR filter coefficients using bilinear

SYNOPSIS ^

function [ao,bo] = pzm2ab(pzm, fs)

DESCRIPTION ^

 PZM2AB convert pzmodel to IIR filter coefficients using bilinear
 transform.
 
 usage: [a,b] = pzm2ab(pzm, fs)
 
 M Hewitson 03-04-07
 
 $Id: pzm2ab.html,v 1.2 2007/07/10 05:37:12 hewitson Exp $

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 function [ao,bo] = pzm2ab(pzm, fs)
0002 % PZM2AB convert pzmodel to IIR filter coefficients using bilinear
0003 % transform.
0004 %
0005 % usage: [a,b] = pzm2ab(pzm, fs)
0006 %
0007 % M Hewitson 03-04-07
0008 %
0009 % $Id: pzm2ab.html,v 1.2 2007/07/10 05:37:12 hewitson Exp $
0010 %
0011 
0012 disp(sprintf('$$$ converting %s', get(pzm, 'name')))
0013 
0014 gain  = get(pzm, 'gain');
0015 poles = get(pzm, 'poles');
0016 zeros = get(pzm, 'zeros');
0017 np = length(poles);
0018 nz = length(zeros);
0019 
0020 ao = [];
0021 bo = [];
0022 
0023 % First we should do complex pole/zero pairs
0024 cpoles = [];
0025 for j=1:np
0026   if get(poles(j), 'q') > 0.5
0027     cpoles = [cpoles poles(j)];
0028   end
0029 end
0030 czeros = [];
0031 for j=1:nz
0032   if get(zeros(j), 'q') > 0.5
0033     czeros = [czeros zeros(j)];
0034   end
0035 end
0036 
0037 czi = 1;
0038 for j=1:length(cpoles)
0039   if czi <= length(czeros)
0040     % we have a pair
0041     p = cpoles(j);
0042     z = czeros(czi);
0043     
0044     [ai,bi] = cpolezero(get(p, 'f'), get(p, 'q'), get(z, 'f'), get(z, 'q'), fs);
0045     if ~isempty(ao)>0
0046       [ao,bo] = abcascade(ao,bo,ai,bi);
0047     else
0048       ao = ai;
0049       bo = bi;
0050     end
0051     
0052     % increment zero counter
0053     czi = czi + 1;
0054   end
0055 end
0056 
0057 if length(cpoles) > length(czeros)
0058   % do remaining cpoles
0059   for j=length(czeros)+1:length(cpoles)
0060     disp('$$$ computing complex pole');
0061     [ai,bi] = cp2iir(cpoles(j), fs);
0062     if ~isempty(ao)>0
0063       [ao,bo] = abcascade(ao,bo,ai,bi);
0064     else
0065       ao = ai;
0066       bo = bi;
0067     end
0068   end
0069 else
0070   % do remaining czeros
0071   for j=length(cpoles)+1:length(czeros)
0072     disp('$$$ computing complex zero');
0073     [ai,bi] = cz2iir(czeros(j), fs);
0074     if ~isempty(ao)>0
0075       [ao,bo] = abcascade(ao,bo,ai,bi);
0076     else
0077       ao = ai;
0078       bo = bi;
0079     end
0080   end
0081 end
0082   
0083 % Now do the real poles and zeros
0084 
0085 
0086 for j=1:np
0087   pole = poles(j);
0088   if get(pole, 'q') <= 0.5
0089     disp('$$$ computing real pole');
0090     [ai,bi] = rp2iir(pole, fs);
0091     if ~isempty(ao)>0
0092       [ao,bo] = abcascade(ao,bo,ai,bi);
0093     else
0094       ao = ai;
0095       bo = bi;
0096     end
0097   end
0098 end
0099 
0100 for j=1:nz
0101   zero = zeros(j);
0102   if get(zeros, 'q') <= 0.5
0103     disp('$$$ computing real zero');
0104     [ai,bi] = rz2iir(zero, fs);
0105     if ~isempty(ao)>0
0106       [ao,bo] = abcascade(ao,bo,ai,bi);
0107     else
0108       ao = ai;
0109       bo = bi;
0110     end
0111   end
0112 end
0113 
0114 % Old version that did everything individually - not so good.
0115 % for j=1:np
0116 %   pole = poles(j);
0117 %   if get(pole, 'q') > 0.5
0118 %     [ai,bi] = cp2iir(pole, fs);
0119 %   else
0120 %     [ai,bi] = rp2iir(pole, fs);
0121 %   end
0122 %   if ~isempty(ao)>0
0123 %     [ao,bo] = abcascade(ao,bo,ai,bi);
0124 %   else
0125 %     ao = ai;
0126 %     bo = bi;
0127 %   end
0128 % end
0129 %
0130 % for j=1:nz
0131 %   zero = zeros(j);
0132 %   if get(zeros, 'q') > 0.5
0133 %     [ai,bi] = cz2iir(zero, fs);
0134 %   else
0135 %     [ai,bi] = rz2iir(zero, fs);
0136 %   end
0137 %   if ~isempty(ao)>0
0138 %     [ao,bo] = abcascade(ao,bo,ai,bi);
0139 %   else
0140 %     ao = ai;
0141 %     bo = bi;
0142 %   end
0143 % end
0144 
0145 ao = ao.*gain;
0146 
0147 function [a,b] = cpolezero(pf, pq, zf, zq, fs)
0148 %
0149 % Return IIR filter coefficients for a complex pole
0150 % and complex zero designed using the bilinear transform.
0151 %
0152 % usage: [a,b] = cpolezero(pf, pq, zf, zq, fs)
0153 %
0154 %
0155 %  M Hewitson 2003-02-18
0156 %
0157 
0158 disp('$$$ computing complex pole/zero pair');
0159 
0160 wp = pf*2*pi;
0161 wp2 = wp^2;
0162 wz = zf*2*pi;
0163 wz2 = wz^2;
0164 
0165 k = 4*fs*fs + 2*wp*fs/pq + wp2;
0166 
0167 a(1) = (4*fs*fs + 2*wz*fs/zq + wz2)/k;
0168 a(2) = (2*wz2 - 8*fs*fs)/k;
0169 a(3) = (4*fs*fs - 2*wz*fs/zq + wz2)/k;
0170 b(1) = 1;
0171 b(2) = (2*wp2 - 8*fs*fs)/k;
0172 b(3) = (wp2 + 4*fs*fs - 2*wp*fs/pq)/k;
0173 
0174 % normalise dc gain to 1
0175 g = iirdcgain(a,b);
0176 a = a / g;
0177 
0178 
0179 function g = iirdcgain(a,b)
0180 
0181 % Work out the DC gain of an IIR
0182 % filter given the coefficients.
0183 %
0184 % usage:
0185 %   g = iirdcgain(a,b)
0186 %
0187 % inputs:
0188 %   a - numerator coefficients
0189 %   b - denominator coefficients
0190 %
0191 % outputs:
0192 %   g - gain
0193 %
0194 % M Hewitson 03-07-02
0195 %
0196 %
0197 % $Id: pzm2ab.html,v 1.2 2007/07/10 05:37:12 hewitson Exp $
0198 %
0199 
0200 g = 0;
0201 
0202 suma = sum(a);
0203 if(length(b)>1)
0204   sumb = sum(b);
0205   g = suma / sumb;
0206 else
0207   g = suma;
0208 end
0209 
0210 % END

Generated on Wed 04-Jul-2007 19:03:10 by m2html © 2003