0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022 function varargout = pzm2ab(varargin)
0023
0024 import utils.const.*
0025
0026
0027 if utils.helper.isinfocall(varargin{:})
0028 varargout{1} = getInfo(varargin{3});
0029 return
0030 end
0031
0032
0033 in_names = cell(size(varargin));
0034 for ii = 1:nargin,in_names{ii} = inputname(ii);end
0035
0036
0037 [pzm, pzm_invars, fs] = utils.helper.collect_objects(varargin(:), 'pzmodel', in_names);
0038
0039
0040 if numel(pzm) ~= 1
0041 error('### Please use only one PZ-Model.');
0042 end
0043 if numel(fs) ~= 1 && ~isnumeric(fs)
0044 error('### Please define ''fs''.');
0045 else
0046 fs = fs{1};
0047 end
0048
0049 gain = pzm.gain;
0050 poles = pzm.poles;
0051 zeros = pzm.zeros;
0052 np = length(poles);
0053 nz = length(zeros);
0054
0055 ao = [];
0056 bo = [];
0057
0058 utils.helper.msg(utils.const.msg.OPROC1, 'converting %s', pzm.name)
0059
0060
0061 cpoles = [];
0062 for j=1:np
0063 if poles(j).q > 0.5
0064 cpoles = [cpoles poles(j)];
0065 end
0066 end
0067 czeros = [];
0068 for j=1:nz
0069 if zeros(j).q > 0.5
0070 czeros = [czeros zeros(j)];
0071 end
0072 end
0073
0074 czi = 1;
0075 for j=1:length(cpoles)
0076 if czi <= length(czeros)
0077
0078 p = cpoles(j);
0079 z = czeros(czi);
0080
0081 [ai,bi] = cpolezero(p.f, p.q, z.f, z.q, fs);
0082 if ~isempty(ao)>0
0083 [ao,bo] = pzmodel.abcascade(ao,bo,ai,bi);
0084 else
0085 ao = ai;
0086 bo = bi;
0087 end
0088
0089
0090 czi = czi + 1;
0091 end
0092 end
0093
0094 if length(cpoles) > length(czeros)
0095
0096 for j=length(czeros)+1:length(cpoles)
0097 utils.helper.msg(msg.OPROC2, 'computing complex pole');
0098 [ai,bi] = cp2iir(cpoles(j), fs);
0099 if ~isempty(ao)>0
0100 [ao,bo] = pzmodel.abcascade(ao,bo,ai,bi);
0101 else
0102 ao = ai;
0103 bo = bi;
0104 end
0105 end
0106 else
0107
0108 for j=length(cpoles)+1:length(czeros)
0109 utils.helper.msg(msg.OPROC2, 'computing complex zero');
0110 [ai,bi] = cz2iir(czeros(j), fs);
0111 if ~isempty(ao)>0
0112 [ao,bo] = pzmodel.abcascade(ao,bo,ai,bi);
0113 else
0114 ao = ai;
0115 bo = bi;
0116 end
0117 end
0118 end
0119
0120
0121 for j=1:np
0122 pole = poles(j);
0123 if isnan(pole.q) || pole.q < 0.5
0124 utils.helper.msg(msg.OPROC2, 'computing real pole');
0125 [ai,bi] = rp2iir(pole, fs);
0126 if ~isempty(ao)>0
0127 [ao,bo] = pzmodel.abcascade(ao,bo,ai,bi);
0128 else
0129 ao = ai;
0130 bo = bi;
0131 end
0132 end
0133 end
0134
0135 for j=1:nz
0136 zero = zeros(j);
0137 if isnan(zero.q) || zero.q < 0.5
0138 utils.helper.msg(msg.OPROC2, 'computing real zero');
0139 [ai,bi] = rz2iir(zero, fs);
0140 if ~isempty(ao)>0
0141 [ao,bo] = pzmodel.abcascade(ao,bo,ai,bi);
0142 else
0143 ao = ai;
0144 bo = bi;
0145 end
0146 end
0147 end
0148
0149 ao = ao.*gain;
0150
0151 varargout{1} = ao;
0152 varargout{2} = bo;
0153 end
0154
0155
0156
0157
0158
0159
0160
0161
0162
0163
0164
0165
0166
0167
0168
0169
0170
0171
0172
0173 function [a,b] = cpolezero(pf, pq, zf, zq, fs)
0174
0175 utils.helper.msg(msg.OPROC1, 'computing complex pole/zero pair');
0176
0177 wp = pf*2*pi;
0178 wp2 = wp^2;
0179 wz = zf*2*pi;
0180 wz2 = wz^2;
0181
0182 k = 4*fs*fs + 2*wp*fs/pq + wp2;
0183
0184 a(1) = (4*fs*fs + 2*wz*fs/zq + wz2)/k;
0185 a(2) = (2*wz2 - 8*fs*fs)/k;
0186 a(3) = (4*fs*fs - 2*wz*fs/zq + wz2)/k;
0187 b(1) = 1;
0188 b(2) = (2*wp2 - 8*fs*fs)/k;
0189 b(3) = (wp2 + 4*fs*fs - 2*wp*fs/pq)/k;
0190
0191
0192 g = iirdcgain(a,b);
0193 a = a / g;
0194 end
0195
0196
0197
0198
0199
0200
0201
0202
0203
0204
0205
0206
0207
0208
0209
0210
0211
0212
0213
0214
0215 function g = iirdcgain(a,b)
0216 suma = sum(a);
0217 if(length(b)>1)
0218 sumb = sum(b);
0219 g = suma / sumb;
0220 else
0221 g = suma;
0222 end
0223 end
0224
0225
0226
0227
0228
0229
0230
0231
0232
0233
0234
0235
0236 function ii = getInfo(varargin)
0237 if nargin == 1 && strcmpi(varargin{1}, 'None')
0238 sets = {};
0239 pl = [];
0240 else
0241 sets = {'Default'};
0242 pl = getDefaultPlist;
0243 end
0244
0245 ii = minfo(mfilename, 'pzmodel', '', utils.const.categories.internal, '$Id: pzm2ab.m,v 1.4 2008/09/04 15:29:31 ingo Exp $', sets, pl);
0246 end
0247
0248
0249
0250
0251
0252
0253
0254
0255
0256
0257
0258
0259 function plo = getDefaultPlist()
0260 plo = plist();
0261 end
0262