classdef (Hidden = true) tfdata < ltipack.ltidata
   % Class definition for @tfdata (transfer function data)

   %   Author(s): P. Gahinet
   %   Copyright 1986-2012 The MathWorks, Inc.

   % IMPORTANT: Code assumes that the NUM,DEN vectors for each I/O pair
   % have equal length and that NUM(1) or DEN(1) is nonzero. Make sure
   % all methods return NUM,DEN arrays that comply with this assumption
   %
   % NOTE: Delay.IO is used to store the matrix of pairwise I/O delays 
   % and Delay.Internal is not used. TFDATA objects do not support 
   % fractional delays and the NUM,DEN arrays are always commensurate 
   % with the I/O size. The rationale for limited delay support is:
   %  * Internal delay support requires multiplying and inverting transfer
   %    matrices
   %  * Setting/getting the nominal NUM,DEN requires complex MINREAL
   %    operations to preserve order:
   %       set H11_nom:   H11(s) = H11_nom(s) - dH11(s)
   %       get H11_nom:   H11(s) + dH11(s) = (H11_nom - dH11) + dH11

   properties
      num  % numerator data
      den  % denominator data
   end

   methods
      function D = tfdata(num,den,Ts)
         % Constructs @tfdata instance
         if nargin==3
            [ny,nu] = size(num);
            D.num = num;
            D.den = den;
            D.Ts = Ts;
            D.Delay = ltipack.utDelayStruct(ny,nu,false);
         end
      end
      
      function D = scaleTime(D,sf)
         % Time vector scaling tnew = sf * told
         D = scaleTime@ltipack.ltidata(D,sf);
         if D.Ts==0
            % Reflect change s_new = s_old/sf (a_n -> a_n sf^n)
            num = D.num;  den = D.den; %#ok<*PROPLC>
            % Note: num{ct} and den{ct} have same length
            lden = cellfun(@length,den); % num and den have same length
            powSF = cumprod([1 sf(:,ones(max(lden(:))-1,1))]);
            for ct=1:numel(num)
               % Note: Leading coefficients preserved
               num{ct} = num{ct} ./ powSF(:,1:lden(ct));
               den{ct} = den{ct} ./ powSF(:,1:lden(ct));
            end
            D.num = num;  D.den = den;
         end
      end
      
      function sys = utCreateLTI(D)
         % Creates LTI object from low-level @xxxdata representation
         sys = tf.make(D);
      end
      
      function D = checkData(D)
         % Checks that num,den data is consistent and NaN free. 
         
         % Determine I/O size from NUM array
         [Ny,Nu] = size(D.num);
         
         % Allow for scalar expansion of denominator
         if (Ny~=1 || Nu~=1) && isscalar(D.den)
            % TF(NUM,[1 0]) (common denominator)
            D.den = repmat(D.den,Ny,Nu);
         end
         
         % Check I/O size consistency
         if ~isequal(size(D.den),[Ny Nu])
            ctrlMsgUtils.error('Control:ltiobject:tfProperties2')
         end
         
         % Checks for NaNs
         Num = D.num;
         Den = D.den;
         for ctio=1:numel(Num)
            if ~(allfinite(Num{ctio}) && allfinite(Den{ctio}))
               D.num(:) = {NaN};
               D.den(:) = {1};
               D.Delay.IO(:) = 0;
               D.Delay.Input(:) = 0;
               D.Delay.Output(:) = 0;
               break
            end
         end
      end
      
      %-----------------
      function [ny,nu] = iosize(D)
         % Returns I/O size.
         %   [NY,NU] = IOSIZE(SYS)
         %   S = IOSIZE(SYS) returns S = [NY NU].
         [ny,nu] = size(D.num);  % NOTE: Data defines I/O size
         if nargout<2
            ny = [ny nu];
         end
      end
      
      function D = transpose(D)
         % Transposition of transfer functions.
         D.Delay = ltipack.util.transposeDelay(D.Delay);
         D.num = D.num.';
         D.den = D.den.';
      end
      
      function D = utGrowIO(D,ny,nu)
         % Grows I/O size of TF model (used in SYS(I,J)=RHS)
         [ny0,nu0] = size(D.num);
         if ny>ny0 || nu>nu0
            num = cell(ny,nu);  num(:) = {0};
            den = cell(ny,nu);  den(:) = {1};
            num(1:ny0,1:nu0) = D.num;
            den(1:ny0,1:nu0) = D.den;
            D.num = num;
            D.den = den;
            % Delays
            D.Delay.IO(ny,nu) = 0;
            if nu>nu0
               D.Delay.Input(nu0+1:nu,:) = 0;
            end
            if ny>ny0
               D.Delay.Output(ny0+1:ny,:) = 0;
            end
         end
      end
      
      function D = createGain(Dref,G)
         % Wraps static gain into @tfdata object
         [ny,nu] = size(G);
         D = ltipack.tfdata(num2cell(G),repmat({1},ny,nu),Dref.Ts);
      end

      function [Dsim,dt,tf,SimInfo] = trespSetUp(D,RespType,dt,tf,SimInfo)
         % Build discrete model for STEP and IMPULSE simulation
         if D.Ts==0 || ~isreal(D)
            % Continuous-time TF or discrete complex-valued TF: Delegate to ssdata
            [Dsim,dt,tf,SimInfo] = trespSetUp(ss(D),RespType,dt,tf,SimInfo);
         else
            % Discrete-time real case: simulate directly in TF form
            Dsim = D;
            % Compute output offset y0
            u0 = SimInfo.uinit;
            if any(u0)
               y0 = dcgain(Dsim) * u0;
               if ~allfinite(y0)
                  % Can't compute IC that ensures steady state with specified
                  % offsets (e.g., due to integrators).
                  error(message('Control:analysis:step2'))
               end
            else
               [ny,~] = iosize(Dsim);
               y0 = zeros(ny,1);
            end
            SimInfo.y0 = y0;
         end
      end

      function [D1,D2,WARN,SepInfo] = stabsep(D,varargin)
         [D1,D2,WARN] = stabsep(ss(D),varargin{:});
         D1 = tf(D1);   D2 = tf(D2);
         SepInfo = [];
      end

      function [D1,D2,WARN,SepInfo] = freqsep(D,varargin)
         [D1,D2,WARN] = freqsep(ss(D),varargin{:});
         D1 = tf(D1);
         D2 = tf(D2);
         SepInfo = [];
      end

      function [D,D0,INFO] = modsep(D,opt)
         % To avoid introducing internal delays in SS realization, cache
         % delays and zero them out before conversion to SS
         Delay = D.Delay;
         [ny,nu] = iosize(D);
         D.Delay = ltipack.utDelayStruct(ny,nu,false);
         % Compute decomposition in state space
         [Dss,D0ss] = modsep(ss(D),opt);
         % Convert back to TF and restore original delays
         D = ltipack.tfdata.array(size(Dss));
         for ct=1:numel(D)
            D(ct) = tf(Dss(ct));   D(ct).Delay = Delay;
         end
         D0 = tf(D0ss); D0.Delay = Delay;
         INFO = [];
      end

      function [D,D0,INFO] = modsepREG(D,varargin)
         % To avoid introducing internal delays in SS realization, cache
         % delays and zero them out before conversion to SS
         Delay = D.Delay;
         [ny,nu] = iosize(D);
         D.Delay = ltipack.utDelayStruct(ny,nu,false);
         % Compute decomposition in state space
         [Dss,D0ss] = modsepREG(ss(D),varargin{:});
         % Convert back to TF and restore original delays
         D = ltipack.tfdata.array(size(Dss));
         for ct=1:numel(D)
            D(ct) = tf(Dss(ct));   D(ct).Delay = Delay;
         end
         D0 = tf(D0ss); D0.Delay = Delay;
         INFO = [];
      end
      
      function DG = modsum(DH)
         % Sums models in DH.
         DG = DH(1);
         Delay = DG.Delay;
         [ny,nu] = iosize(DG);
         ZeroDelay = ltipack.utDelayStruct(ny,nu,false);
         DG.Delay = ZeroDelay;
         DG = ss(DG);
         for ct=2:numel(DH)
            if isequal(DH(ct).Delay,Delay)
               DH(ct).Delay = ZeroDelay;
               DG = plus(DG,ss(DH(ct)));
            else
               error(message('Control:transformation:modalsum1'))
            end
         end
         DG = tf(DG);
         DG.Delay = Delay;
      end

   end
   
   methods(Static)
      function D = array(size)
         % Create a tfdata array of a given size
         D = ltipack.tfdata.newarray(size);
      end
      
      function D = loadobj(S)
         % Load filter for @tfdata
         if isstruct(S)
            % R2013b: Tags removed
            D = ltipack.transferData(ltipack.tfdata,S);
         else
            D = S;
         end
         if isfield(D.Delay,'Internal')
            % Pre-R2009b: Delay structure had Internal field for all model types
            D.Delay = rmfield(D.Delay,'Internal');
         end
      end
      
      [num,den,NeedsPadding] = utFormatNumDen(num,den,RightFlag)
      [num,den] = utTrimNumDen(num,den,RightTrim)
      
   end
         
   % Protected methods (utilities)
   methods(Access=protected)
      
      function [ro,co,D] = getOrder(D)
         % Computes order of TF models
         %
         %   [RO,CO] = GETORDER(D) computes the orders RO and CO of
         %   row-wise and column-wise state-space realizations of D
         %
         %   [RO,CO,D] = GETORDER(D) also returns the normalized
         %   transfer function with all leading denominator coefficients
         %   equal to 1
         num = D.num; %#ok<*PROP>
         den = D.den;
         [ny,nu] = size(num);

         % Normalize denominators and compute order of hij(s)
         % RE: Normalization needed to identify common den in
         % sys = [tf([1 2 5],-2*[1 2 2]) ; tf(-6*[1 0],4*[1 2 2])];
         Normalized = true;
         ioOrders = zeros(ny,nu);
         for ct=1:ny*nu
            nct = num{ct};
            dct = den{ct};
            den1 = dct(find(dct~=0,1));
            if den1~=1
               num{ct} = nct/den1;
               den{ct} = dct/den1;
               Normalized = false;
            end
            % Order of hij(s)
            if all(nct==0)
               % zero entries contribute no dynamics
               ioOrders(ct) = 0;
            elseif dct(1)~=0
               % proper
               ioOrders(ct) = length(dct)-1;
            else
               ioOrders(ct) = length(nct);
            end
         end

         % Determine row-wise order
         ro = 0;
         for i=1:ny
            jdyn = find(ioOrders(i,:));  % non-static entries
            if length(jdyn)>1 && ltipack.isEqualDen(den{i,jdyn})
               % Common denominator
               ro = ro + ioOrders(i,jdyn(1));
            else
               % Sum orders for hij(s), j = 1:nu
               ro = ro + sum(ioOrders(i,jdyn));
            end
         end

         % Determine column-wise order
         co = 0;
         for j=1:nu
            idyn = find(ioOrders(:,j)); % non-static entries
            if length(idyn)>1 && ltipack.isEqualDen(den{idyn,j})
               % Common denominator
               co = co + ioOrders(idyn(1),j);
            else
               co = co + sum(ioOrders(idyn,j));
            end
         end

         if nargout==3 && ~Normalized
            D.num = num;
            D.den = den;
         end
      end
      
      function boo = utIsIOScaling(D)
         % True if D is of the form diag(dj*exp(-s*tauj)) where dj is a gain
         [ny,nu] = size(D.num);
         if ny~=nu || any(cellfun('length',D.num(:))>1)
            boo = false;
         else
            num = cat(1,D.num{:});
            idx = 1:ny*nu;
            idx(1:nu+1:ny*nu) = [];  % off diagonal entries
            boo = ~(any(D.Delay.IO(idx)) || any(num(idx)));
         end
      end

   end

end

