0001
0002
0003 function subsys = ltpda_ss_pz2ss(pzm)
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028
0029
0030
0031
0032
0033
0034
0035
0036
0037
0038
0039
0040
0041
0042
0043
0044
0045
0046
0047 ALGONAME = mfilename;
0048 VERSION = '$Id: ltpda_ss_pz2ssVjordan.html,v 1.1 2008/03/01 12:29:32 hewitson Exp $';
0049 display(['starting ' mfilename]);
0050
0051
0052
0053
0054 zeros = pzm.zeros;
0055 b = 1;
0056 for i = 1:length(poles)
0057 w0 = 2*pi*get(poles(i), 'f');
0058 if isequal(p.name,'real pole')
0059 b_p = [1 -1/(w0)];
0060 elseif isequal(p.name,'complex pole')
0061 q = get(poles(i), 'q');
0062 b_p = [1 -1/(w0*q) 1/(w0^2)];
0063 else
0064 error('type of pole not expected')
0065 end
0066 b = conv(b,b_p);
0067 end
0068
0069 poles = pzm.poles;
0070 a = 1;
0071 for i = 1:length(zeros)
0072 w0 = 2*pi*get(zeros(i), 'f');
0073 if isequal(p.name,'real pole')
0074 a_z = [1 -1/(w0)];
0075 elseif isequal(p.name,'complex pole')
0076 q = get(zeros(i), 'q');
0077 a_z = [1 -1/(w0*q) 1/(w0^2)];
0078 else
0079 error('type of zero not expected')
0080 end
0081 a = conv(a,a_z);
0082 end
0083
0084 nss = size(b,2)-1;
0085 A = zeros(nss, nss);
0086 B = zeros(nss, 1);
0087 C = zeros(1, nss);
0088
0089
0090
0091 [r,p,k] = residue(a,b);
0092 D = k;
0093
0094 blockNumber = 1;
0095 numbering = zeros(1,nss);
0096 isReal = zeros(1,nss);
0097 for i=1:(nss-1)
0098 if p(i)==p(i+1)
0099 numbering(i) = blockNumber;
0100 numbering(i+1) = blockNumber;
0101 elseif p(i)==conj(p(i+1))
0102 numbering(i) = blockNumber;
0103 numbering(i+1) = blockNumber;
0104 else
0105 numbering = numbering + 1;
0106 end
0107 isReal(i) = (imag(p(i))==0);
0108 end
0109 isReal(nss) = (imag(p(nss))==0);
0110
0111 for i=blockNumber
0112 indices = find(numbering==blockNumber);
0113 res{i} = r(indices);
0114 pol{i} = p(indices);
0115 isRealLoc(i) = isReal(indices(1));
0116 nssLoc(i) = elngth(indices);
0117 end
0118
0119 for i=blockNumber
0120 if isRealLoc(i)
0121 A{i} = diag(pol{i}) + [zeros(nssLoc(i)-1,1) eye(nssLoc(i)-1); zeros(1,nssLoc(i)) ];
0122 B{i} = res(i);
0123 C{i} = [1 zeros(1,nssLoc(i)-1)];
0124 else
0125 nss_loc_2 = nss_loc(i)/2;
0126 if not(nnsloc2==floor(nss_loc2))
0127 error('complex pole has odd number of roots!')
0128 end
0129 pol_loc = pol{i};
0130 pol_loc_2 = pol_loc(1:nss_loc_2);
0131
0132
0133
0134 A_loc_2 = diag(pol_loc_2.*conj(pol_loc_2));
0135 if norm(imag(A_loc_2))>10e-15
0136 error('error because calculus in conjugate was wrong')
0137 end
0138 A_loc_2 = real(A_loc_2);
0139 A_loc_1 = eye(nss_loc_2) + [zeros(nss_loc_2,1) eye(nss_loc_2-1); zeros(1,nss_loc_2) ];
0140 A_loc_3 = diag(2*Real(pol_loc_2));
0141 A_loc = [zeros(nss_loc_2,nss_loc_2) A_loc_1 ; A_loc_2 A_loc_3];
0142
0143 B_loc = [ 2*Real(pol_loc_2) pol_loc_2.*conj(pol_loc_2)];
0144
0145 C_loc = [1 zeros(1,nssLoc(i)-1)];
0146
0147 A{i} = A_loc;
0148 B{i} = B_loc;
0149 C{i} = C_loc;
0150 end
0151 end
0152
0153
0154 A2=[];
0155 B2=[];
0156 C2=[];
0157 for i=1:blockNumber
0158 A2 = blkdiag(A2,A{i});
0159 B2 = vertcat(B2,B{i});
0160 C2 = horzcat(C2,C{i});
0161 end
0162
0163 subsys = plist('TYPE', 'SUMBSYSTEM' ,'NAME', iir.name ,'TIMESTEP', 1/(iir.fs), ...
0164 'XISOUTPUT', 0,'YISOUTPUT',1,'XINI', zeros(nss,1) , ...
0165 'PARAMNAMES', {},'PARAMVALUES', [],'PARAMSIGMAS', [],...
0166 'NBINPUTS', f.ntaps ,'INPUTNAMES', {'U_1'} ,'INPUTSIZES', [1] , 'INPUTISUSED', [1] ,...
0167 'AMAT', {A2} ,'BMATS', {B2} ,'CMAT', {C2} ,'DMATS', {D} );
0168 end
0169
0170
0171 subsys = plist('TYPE', 'SUMBSYSTEM' ,'NAME', pzm.name ,'TIMESTEP', 0 , ...
0172 'XISOUTPUT', 0,'YISOUTPUT',1,'XINI', zeros(nss,1) , ...
0173 'PARAMNAMES', {},'PARAMVALUES', [],'PARAMSIGMAS', [],...
0174 'NBINPUTS', 1 ,'INPUTNAMES', {'U_1'} ,'INPUTSIZES', [1] , 'INPUTISUSED', [1] ,...
0175 'AMAT', {A} ,'BMATS', {B} ,'CMAT', {C} ,'DMATS', {D} );
0176 end