0001 function subsys = ltpda_ss_miir2ss(iir)
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025 ALGONAME = mfilename;
0026 VERSION = '$Id: ltpda_ss_miir2ssVjordan.html,v 1.1 2008/03/01 12:29:32 hewitson Exp $';
0027 display(['starting ' mfilename]);
0028 a = iir.a;
0029 b = iir.b;
0030
0031
0032 nss = size(b,2)-1;
0033
0034
0035
0036
0037
0038 [r,p,k] = residue(a,b);
0039 D = k;
0040
0041 blockNumber = 1;
0042 numbering = zeros(1,nss);
0043 isReal = zeros(1,nss);
0044 for i=1:(nss-1)
0045 if norm(p(i)-p(i+1))/norm(p(i))<1e-2
0046 numbering(i) = blockNumber;
0047 numbering(i+1) = blockNumber;
0048 elseif norm(p(i)-conj(p(i+1)))/norm(p(i))<1e-2
0049 numbering(i) = blockNumber;
0050 numbering(i+1) = blockNumber;
0051 else
0052 numbering(i) = blockNumber;
0053 blockNumber = blockNumber + 1;
0054 numbering(i+1) = blockNumber;
0055 end
0056 isReal(i) = (abs(imag(p(i)))<1e-3);
0057 end
0058 isReal(nss) = (imag(p(nss))==0);
0059
0060 for i=1:blockNumber
0061 indices = find(numbering==i);
0062 res{i} = r(indices);
0063 pol{i} = p(indices);
0064 isRealBlock(i) = isReal(indices(1));
0065 nssLoc(i) = length(indices);
0066 end
0067
0068 for i=1:blockNumber
0069 if isRealBlock(i)
0070 A{i} = diag(pol{i}) + [zeros(nssLoc(i)-1,1) eye(nssLoc(i)-1); zeros(1,nssLoc(i)) ];
0071 B{i} = res{i};
0072 C{i} = [1 zeros(1,nssLoc(i)-1)];
0073 else
0074 nss_loc_2 = nssLoc(i)/2;
0075 positions_2 = (1:nss_loc_2)*2;
0076 if not(nss_loc_2==floor(nss_loc_2))
0077 error('complex pole has odd number of roots!')
0078 end
0079
0080 pol_loc = pol{i};
0081 res_loc = res{i};
0082 [a_loc, b_loc] = residue(res_loc, pol_loc,0);
0083 a_loc = real(a_loc);
0084
0085 pol_loc_2 = pol_loc(positions_2);
0086 pol_loc_mean = mean(pol_loc_2);
0087 b = 1;
0088 a = real(a_loc);
0089 for j=1:nss_loc_2
0090 b = conv(b, [1 2*real(pol_loc_mean) norm(pol_loc_mean)^2]);
0091 end
0092
0093
0094
0095
0096
0097
0098 [res_loc, pol_loc ] = residue(a, b);
0099 pol_loc_2 = pol_loc(positions_2);
0100 res_loc_2 = res_loc(positions_2);
0101
0102
0103 A_loc_2 = diag(pol_loc_2.*conj(pol_loc_2));
0104 if norm(imag(A_loc_2))>10e-15
0105 error('error because calculus in conjugate was wrong')
0106 end
0107 A_loc_2 = real(A_loc_2);
0108 A_loc_1 = eye(nss_loc_2) + [zeros(nss_loc_2-1,1) eye(nss_loc_2-1); zeros(1,nss_loc_2) ];
0109 A_loc_3 = diag(2*real(pol_loc_2));
0110 A_loc = [zeros(nss_loc_2,nss_loc_2) A_loc_1 ; A_loc_2 A_loc_3];
0111
0112 B_loc = [ 2*real(res_loc_2); res_loc_2.*conj(res_loc_2)];
0113
0114 C_loc = [1 zeros(1,nssLoc(i)-1)];
0115 A{i} = A_loc;
0116 B{i} = B_loc;
0117 C{i} = C_loc;
0118
0119
0120
0121
0122 end
0123 end
0124
0125
0126 A2=[];
0127 B2=[];
0128 C2=[];
0129 for i=1:blockNumber
0130 A2 = real(blkdiag(A2,A{i}));
0131 B2 = vertcat(B2,B{i});
0132 C2 = horzcat(C2,C{i});
0133 end
0134
0135 subsys = plist('TYPE', 'SUMBSYSTEM' ,'NAME', iir.name ,'TIMESTEP', 1/(iir.fs), ...
0136 'XISOUTPUT', 0,'YISOUTPUT',1,'XINI', zeros(nss,1) , ...
0137 'PARAMNAMES', {},'PARAMVALUES', [],'PARAMSIGMAS', [],...
0138 'NBINPUTS', iir.ntaps ,'INPUTNAMES', {'U_1'} ,'INPUTSIZES', [1] , 'INPUTISUSED', [1] ,...
0139 'AMAT', {A2} ,'BMATS', {B2} ,'CMAT', {C2} ,'DMATS', {D} );
0140
0141 ltpda_ss_check(subsys);
0142 end