Mercurial > hg > ltpda

comparison m-toolbox/classes/@pzmodel/pzm2ab.m @ 0:f0afece42f48

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Import.
author Daniele Nicolodi <nicolodi@science.unitn.it>
date Wed, 23 Nov 2011 19:22:13 +0100
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:f0afece42f48
1 % PZM2AB convert pzmodel to IIR filter coefficients using bilinear transform.
2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3 %
4 % DESCRIPTION: PZM2AB convert pzmodel to IIR filter coefficients using bilinear
5 % transform.
6 %
7 % CALL: [a,b] = pzm2ab(pzm, fs)
8 %
9 % VERSION: $Id: pzm2ab.m,v 1.10 2009/08/31 17:11:41 ingo Exp $
10 %
11 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
12
13 function varargout = pzm2ab(varargin)
14
15 import utils.const.*
16
17 %%% Collect input variable names
18 in_names = cell(size(varargin));
19 for ii = 1:nargin,in_names{ii} = inputname(ii);end
20
21 %%% Collect all AOs and plists
22 [pzm, pzm_invars, fs] = utils.helper.collect_objects(varargin(:), 'pzmodel', in_names);
23
24 %%% Check inputs
25 if numel(pzm) ~= 1
26 error('### Please use only one PZ-Model.');
27 end
28 if numel(fs) ~= 1 && ~isnumeric(fs)
29 error('### Please define ''fs''.');
30 else
31 fs = fs{1};
32 end
33
34 gain = pzm.gain;
35 poles = pzm.poles;
36 zeros = pzm.zeros;
37 np = numel(poles);
38 nz = numel(zeros);
39
40 ao = [];
41 bo = [];
42
43 utils.helper.msg(utils.const.msg.OPROC1, 'converting %s', pzm.name)
44
45 % First we should do complex pole/zero pairs
46 cpoles = [];
47 for j=1:np
48 if poles(j).q > 0.5
49 cpoles = [cpoles poles(j)];
50 end
51 end
52 czeros = [];
53 for j=1:nz
54 if zeros(j).q > 0.5
55 czeros = [czeros zeros(j)];
56 end
57 end
58
59 czi = 1;
60 for j=1:length(cpoles)
61 if czi <= length(czeros)
62 % we have a pair
63 p = cpoles(j);
64 z = czeros(czi);
65
66 [ai,bi] = cpolezero(p.f, p.q, z.f, z.q, fs);
67 if ~isempty(ao)
68 [ao,bo] = pzmodel.abcascade(ao,bo,ai,bi);
69 else
70 ao = ai;
71 bo = bi;
72 end
73
74 % increment zero counter
75 czi = czi + 1;
76 end
77 end
78
79 if length(cpoles) > length(czeros)
80 % do remaining cpoles
81 for j=length(czeros)+1:length(cpoles)
82 utils.helper.msg(msg.OPROC2, 'computing complex pole');
83 [ai,bi] = cp2iir(cpoles(j), fs);
84 if ~isempty(ao)
85 [ao,bo] = pzmodel.abcascade(ao,bo,ai,bi);
86 else
87 ao = ai;
88 bo = bi;
89 end
90 end
91 else
92 % do remaining czeros
93 for j=length(cpoles)+1:length(czeros)
94 utils.helper.msg(msg.OPROC2, 'computing complex zero');
95 [ai,bi] = cz2iir(czeros(j), fs);
96 if ~isempty(ao)
97 [ao,bo] = pzmodel.abcascade(ao,bo,ai,bi);
98 else
99 ao = ai;
100 bo = bi;
101 end
102 end
103 end
104
105 % Now do the real poles and zeros
106 for j=1:np
107 pole = poles(j);
108 if isnan(pole.q) || pole.q < 0.5
109 utils.helper.msg(msg.OPROC2, 'computing real pole');
110 [ai,bi] = rp2iir(pole, fs);
111 if ~isempty(ao)
112 [ao,bo] = pzmodel.abcascade(ao,bo,ai,bi);
113 else
114 ao = ai;
115 bo = bi;
116 end
117 end
118 end
119
120 for j=1:nz
121 zero = zeros(j);
122 if isnan(zero.q) || zero.q < 0.5
123 utils.helper.msg(msg.OPROC2, 'computing real zero');
124 [ai,bi] = rz2iir(zero, fs);
125 if ~isempty(ao)
126 [ao,bo] = pzmodel.abcascade(ao,bo,ai,bi);
127 else
128 ao = ai;
129 bo = bi;
130 end
131 end
132 end
133
134 ao = ao.*gain;
135
136 varargout{1} = ao;
137 varargout{2} = bo;
138 end
139
140 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
141 % Local Functions %
142 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
143
144 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
145 %
146 % FUNCTION: cpolezero
147 %
148 % DESCRIPTION: Return IIR filter coefficients for a complex pole and
149 % complex zero designed using the bilinear transform.
150 %
151 % CALL: [a,b] = cpolezero(pf, pq, zf, zq, fs)
152 %
153 % HISTORY: 18-02-2007 M Hewitson
154 % Creation.
155 %
156 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
157
158 function [a,b] = cpolezero(pf, pq, zf, zq, fs)
159
160 import utils.const.*
161 utils.helper.msg(msg.OPROC1, 'computing complex pole/zero pair');
162
163 wp = pf*2*pi;
164 wp2 = wp^2;
165 wz = zf*2*pi;
166 wz2 = wz^2;
167
168 k = 4*fs*fs + 2*wp*fs/pq + wp2;
169
170 a(1) = (4*fs*fs + 2*wz*fs/zq + wz2)/k;
171 a(2) = (2*wz2 - 8*fs*fs)/k;
172 a(3) = (4*fs*fs - 2*wz*fs/zq + wz2)/k;
173 b(1) = 1;
174 b(2) = (2*wp2 - 8*fs*fs)/k;
175 b(3) = (wp2 + 4*fs*fs - 2*wp*fs/pq)/k;
176
177 % normalise dc gain to 1
178 g = iirdcgain(a,b);
179 a = a / g;
180 end
181
182
183 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
184 %
185 % FUNCTION: iirdcgain
186 %
187 % DESCRIPTION: Work out the DC gain of an IIR filter given the coefficients.
188 %
189 % CALL: g = iirdcgain(a,b)
190 %
191 % INPUTS: a - numerator coefficients
192 % b - denominator coefficients
193 %
194 % OUTPUTS g - gain
195 %
196 % HISTORY: 03-07-2002 M Hewitson
197 % Creation.
198 %
199 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
200
201 function g = iirdcgain(a,b)
202 suma = sum(a);
203 if(length(b)>1)
204 sumb = sum(b);
205 g = suma / sumb;
206 else
207 g = suma;
208 end
209 end