Home > m > timetools > statespacefunctions > ltpda_ss_pz2ssVjordan.m

ltpda_ss_pz2ssVjordan

PURPOSE ^

% conversion from a pole zero model

SYNOPSIS ^

function subsys = ltpda_ss_pz2ss(pzm)

DESCRIPTION ^

% conversion from a pole zero model

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 %% conversion from a pole zero model
0002 
0003 function subsys = ltpda_ss_pz2ss(pzm)
0004 % ltpda_pz2ss converts a pole-zero model into a subsystem.
0005 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0006 %
0007 % DESCRIPTION: ltpda_pz2ss converts a pole-zero model into a subsystem.
0008 %
0009 % CALL: subsys = ltpda_pz2ss(pzm)
0010 %
0011 % INPUTS: 3 possibilites
0012 % pzm - a pole zero model object
0013 %
0014 % OUTPUTS: subsys - subsystem dscribed by a plist
0015 %
0016 % subsystem plist format (More updated infos might be available in the
0017 % ltpda_ss_check function):
0018 % plist('TYPE', TYPE ,'NAME', NAME ,'TIMESTEP', TIMESTEP ,...
0019 %     'PARAMNAMES', PARAMNAMES ,'PARAMVALUES', PARAMVALUES ,'PARAMSIGMAS', PARAMSIGMAS ,...
0020 %     'NBINPUTS', NBINPUTS ,'INPUTNAMES', INPUTNAMES ,'INPUTSIZES', INPUTSIZES ,'INPUTISUSED', INPUTISUSED ,...
0021 %     'AMAT', AMAT ,'BMATS', BMATS ,'CMAT', CMAT ,'DMATS', DMATS );% PARAMETERS:
0022 % 'TYPE' should be 'SUBSYSTEM'
0023 % 'NAME' is a string
0024 % 'TIMESTEP' is a real positive integer, set to zero in continuous case
0025 % 'PARAMNAMES' cell array of strings describing parameter variables
0026 % 'PARAMVALUES' array of doubles (means expected of parameters)
0027 % 'PARAMSIGMAS' array of doubles (variance expected of parameters)
0028 % 'NBINPUTS' integer number of inputs
0029 % 'INPUTNAMES' cell array of strings (name of the inputs)
0030 % 'INPUTSIZES'array of integers (nb of dimensions of the inputs),
0031 % 'INPUTISUSED' double array (binary to indicate which inputs are used)
0032 % 'AMAT' cell array (contains A matrix)
0033 % 'BMATS' cell array (contains B matrices)
0034 % 'CMAT' cell array (contains C matrix)
0035 % 'DMATS' cell array (contains D matrices)
0036 % ***** THERE ARE NO DEFAULT PARAMETERS *****
0037 %
0038 % VERSION: $Id: ltpda_ss_pz2ssVjordan.html,v 1.1 2008/03/01 12:29:32 hewitson Exp $
0039 %
0040 % TO DO : lots of testing including numerical checks.
0041 % Discuss the possiblity of jordan block form and PFD form for the output
0042 % pz model
0043 %
0044 % HISTORY: 01-02-2008 A Grynagier
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 % for next part check must be made there is no pole zero cancelation
0053 %% building numerator
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); %multplication of polynomials to get denominator
0067 end
0068 %% building denominator
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); %multplication of polynomials to get numerator
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 % is residue ok??  MUST BE CHANGED !!!! PART WRITTEN TOECONOMISE CODE
0090 % WIRTING: CURRENTLY DENOMINATOR IS EXPANDED AND THEN FACTORIZED AGAIN!!!
0091 [r,p,k] = residue(a,b); %!! here b are for the poles, a the zeros
0092 D = k;
0093 %% detecting blocks of same pole
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); %for last pole which is not indexed by the loop above
0110 %% splitting
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 %% for each block build matrix
0119 for i=blockNumber
0120     if isRealLoc(i) %Real pole
0121         A{i} = diag(pol{i}) + [zeros(nssLoc(i)-1,1) eye(nssLoc(i)-1); zeros(1,nssLoc(i)) ];
0122         B{i} = res(i); %% to check !!!!!! (take opposite ?)
0123         C{i} = [1 zeros(1,nssLoc(i)-1)];
0124     else %complex pole
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         % Building block in A
0132         % pol_loc_2 = pol_loc(1:nss_loc_2,1:nss_loc_2);
0133         % A_loc_2 = pol_loc_2*conj(pol_loc_2);
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         % Building Block in B
0143         B_loc = [ 2*Real(pol_loc_2) pol_loc_2.*conj(pol_loc_2)];
0144         % Building block in C
0145         C_loc = [1 zeros(1,nssLoc(i)-1)];
0146 
0147         A{i} = A_loc;
0148         B{i} = B_loc; %% to check !!!!!!
0149         C{i} = C_loc;
0150     end
0151 end
0152 
0153 %% put all matrices altoghether
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 %% parsing subsystem
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

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