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.m,v 1.4 2007/09/03 09:27:04 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.m,v 1.4 2007/09/03 09:27:04 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 for j=1:np
0085   pole = poles(j);
0086   if isnan(get(pole, 'q')) || get(pole, 'q') < 0.5
0087     disp('$$$ computing real pole');
0088     [ai,bi] = rp2iir(pole, fs);
0089     if ~isempty(ao)>0
0090       [ao,bo] = abcascade(ao,bo,ai,bi);
0091     else
0092       ao = ai;
0093       bo = bi;
0094     end
0095   end
0096 end
0097 
0098 for j=1:nz
0099   zero = zeros(j);
0100   if isnan(get(zero, 'q')) || get(zero, 'q') < 0.5
0101     disp('$$$ computing real zero');
0102     [ai,bi] = rz2iir(zero, fs);
0103     if ~isempty(ao)>0
0104       [ao,bo] = abcascade(ao,bo,ai,bi);
0105     else
0106       ao = ai;
0107       bo = bi;
0108     end
0109   end
0110 end
0111 
0112 % Old version that did everything individually - not so good.
0113 % for j=1:np
0114 %   pole = poles(j);
0115 %   if get(pole, 'q') > 0.5
0116 %     [ai,bi] = cp2iir(pole, fs);
0117 %   else
0118 %     [ai,bi] = rp2iir(pole, fs);
0119 %   end
0120 %   if ~isempty(ao)>0
0121 %     [ao,bo] = abcascade(ao,bo,ai,bi);
0122 %   else
0123 %     ao = ai;
0124 %     bo = bi;
0125 %   end
0126 % end
0127 %
0128 % for j=1:nz
0129 %   zero = zeros(j);
0130 %   if get(zeros, 'q') > 0.5
0131 %     [ai,bi] = cz2iir(zero, fs);
0132 %   else
0133 %     [ai,bi] = rz2iir(zero, fs);
0134 %   end
0135 %   if ~isempty(ao)>0
0136 %     [ao,bo] = abcascade(ao,bo,ai,bi);
0137 %   else
0138 %     ao = ai;
0139 %     bo = bi;
0140 %   end
0141 % end
0142 
0143 ao = ao.*gain;
0144 
0145 function [a,b] = cpolezero(pf, pq, zf, zq, fs)
0146 %
0147 % Return IIR filter coefficients for a complex pole
0148 % and complex zero designed using the bilinear transform.
0149 %
0150 % usage: [a,b] = cpolezero(pf, pq, zf, zq, fs)
0151 %
0152 %
0153 %  M Hewitson 2003-02-18
0154 %
0155 
0156 disp('$$$ computing complex pole/zero pair');
0157 
0158 wp = pf*2*pi;
0159 wp2 = wp^2;
0160 wz = zf*2*pi;
0161 wz2 = wz^2;
0162 
0163 k = 4*fs*fs + 2*wp*fs/pq + wp2;
0164 
0165 a(1) = (4*fs*fs + 2*wz*fs/zq + wz2)/k;
0166 a(2) = (2*wz2 - 8*fs*fs)/k;
0167 a(3) = (4*fs*fs - 2*wz*fs/zq + wz2)/k;
0168 b(1) = 1;
0169 b(2) = (2*wp2 - 8*fs*fs)/k;
0170 b(3) = (wp2 + 4*fs*fs - 2*wp*fs/pq)/k;
0171 
0172 % normalise dc gain to 1
0173 g = iirdcgain(a,b);
0174 a = a / g;
0175 
0176 
0177 function g = iirdcgain(a,b)
0178 
0179 % Work out the DC gain of an IIR
0180 % filter given the coefficients.
0181 %
0182 % usage:
0183 %   g = iirdcgain(a,b)
0184 %
0185 % inputs:
0186 %   a - numerator coefficients
0187 %   b - denominator coefficients
0188 %
0189 % outputs:
0190 %   g - gain
0191 %
0192 % M Hewitson 03-07-02
0193 %
0194 %
0195 % $Id: pzm2ab.m,v 1.4 2007/09/03 09:27:04 hewitson Exp $
0196 %
0197 
0198 g = 0;
0199 
0200 suma = sum(a);
0201 if(length(b)>1)
0202   sumb = sum(b);
0203   g = suma / sumb;
0204 else
0205   g = suma;
0206 end
0207 
0208 % END

Generated on Mon 03-Sep-2007 12:12:34 by m2html © 2003