Home > m > timetools > statespacefunctions > ltpda_ss_miir2ssVjordan.m

ltpda_ss_miir2ssVjordan

PURPOSE ^

ltpda_miir2ss converts a iir filter model into a subsystem.

SYNOPSIS ^

function subsys = ltpda_ss_miir2ss(iir)

DESCRIPTION ^

 ltpda_miir2ss converts a iir filter model into a subsystem.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

 DESCRIPTION: ltpda_miir2ss converts a iir filter model into a subsystem.
 Jordan block transformation was dropped due to numerical issues. This
 version could be greatly enhanced using this feature.

 CALL: subsys = ltpda_miir2ss(iir)

 INPUTS: iir - an infinite inpulsr response model object

 OUTPUTS: subsys - subsystem dscribed by a plist

 ***** THERE ARE NO DEFAULT PARAMETERS *****

 VERSION: $Id: ltpda_ss_miir2ssVjordan.html,v 1.1 2008/03/01 12:29:32 hewitson Exp $

 HISTORY: 7-02-2008 A Grynagier
 2-02-2008 A Grynagier
 
 TO DO : Add the jordan block version as an option
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function subsys = ltpda_ss_miir2ss(iir)
0002 % ltpda_miir2ss converts a iir filter model into a subsystem.
0003 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0004 %
0005 % DESCRIPTION: ltpda_miir2ss converts a iir filter model into a subsystem.
0006 % Jordan block transformation was dropped due to numerical issues. This
0007 % version could be greatly enhanced using this feature.
0008 %
0009 % CALL: subsys = ltpda_miir2ss(iir)
0010 %
0011 % INPUTS: iir - an infinite inpulsr response model object
0012 %
0013 % OUTPUTS: subsys - subsystem dscribed by a plist
0014 %
0015 % ***** THERE ARE NO DEFAULT PARAMETERS *****
0016 %
0017 % VERSION: $Id: ltpda_ss_miir2ssVjordan.html,v 1.1 2008/03/01 12:29:32 hewitson Exp $
0018 %
0019 % HISTORY: 7-02-2008 A Grynagier
0020 % 2-02-2008 A Grynagier
0021 %
0022 % TO DO : Add the jordan block version as an option
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 % for next part check must be made there is no pole zero cancelation
0032 nss = size(b,2)-1;
0033 % A = zeros(nss, nss);
0034 % B = zeros(nss, 1);
0035 % C = zeros(1, nss);
0036 
0037 % is residue ok??
0038 [r,p,k] = residue(a,b); %!! here b are for the poles, a the zeros
0039 D = k;
0040 %% detecting blocks of same pole
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); %for last pole which is not indexed by the loop above
0059 %% splitting
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 %% for each block build matrix
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}; %% to check !!!!!! (take opposite ?)
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 %         a = a/b(1);
0093 %         b = b/b(1);
0094 %         A{i} =  [zeros(nss_loc_2*2-1,1) eye(nss_loc_2*2-1); b(2:length(b))];
0095 %         B{i} = [zeros(nss_loc_2*2-1,1) ;1];
0096 %         C{i} = a;
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         % Building block in A
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         % Building Block in B
0112         B_loc = [ 2*real(res_loc_2); res_loc_2.*conj(res_loc_2)];
0113         % Building block in C
0114         C_loc = [1 zeros(1,nssLoc(i)-1)];
0115         A{i} = A_loc;
0116         B{i} = B_loc; %% to check !!!!!!
0117         C{i} = C_loc;
0118         
0119 
0120         
0121 
0122     end
0123 end
0124 
0125 %% put all matrices altoghether
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

Generated on Fri 22-Feb-2008 23:32:26 by m2html © 2003