0001 function [ao,bo] = pzm2ab(pzm, fs)
0002
0003
0004
0005
0006
0007
0008
0009
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
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
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
0053 czi = czi + 1;
0054 end
0055 end
0056
0057 if length(cpoles) > length(czeros)
0058
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
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
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
0115
0116
0117
0118
0119
0120
0121
0122
0123
0124
0125
0126
0127
0128
0129
0130
0131
0132
0133
0134
0135
0136
0137
0138
0139
0140
0141
0142
0143
0144
0145 ao = ao.*gain;
0146
0147 function [a,b] = cpolezero(pf, pq, zf, zq, fs)
0148
0149
0150
0151
0152
0153
0154
0155
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
0175 g = iirdcgain(a,b);
0176 a = a / g;
0177
0178
0179 function g = iirdcgain(a,b)
0180
0181
0182
0183
0184
0185
0186
0187
0188
0189
0190
0191
0192
0193
0194
0195
0196
0197
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