classdef (CaseInsensitiveProperties, TruncatedProperties, SupportExtensionMethods, ...
      InferiorClasses = {? matlab.graphics.axis.Axes, ? matlab.ui.control.UIAxes}) ...
      frd < numlti & FRDModel
   %FRD  Constructs or converts to Frequency Response Data model.
   %
   %   Frequency Response Data (FRD) models store the frequency response of linear 
   %   systems, for example, experimental data collected with a frequency analyzer.
   %
   %  Construction:
   %    SYS = FRD(RESPONSE,FREQS) creates an FRD model SYS with response data
   %    RESPONSE specified at the frequency points in FREQS. The output SYS 
   %    is an object of class @frd.
   %
   %    SYS = FRD(RESPONSE,FREQS,TS) creates a discrete-time FRD model with
   %    sample time TS (a positive value).
   %
   %    SYS = FRD creates an empty FRD model.
   %
   %    You can set additional model properties by using name/value pairs.
   %    For example,
   %       sys = frd(1:10,1:10,'FrequencyUnit','Hz')
   %    further stipulates that the frequency vector is given in Hz. Type 
   %    "properties(frd)" for a complete list of model properties, and type 
   %       help frd.<PropertyName>
   %    for help on a particular property. For example, "help frd.IODelay" 
   %    provides information about the "IODelay" property.
   %
   %  Data format:
   %    For SISO models, FREQS is a vector of real frequencies, and RESPONSE 
   %    is a vector of frequency response values at these frequencies.
   %
   %    For MIMO FRD models with NY outputs, NU inputs, and NF frequency points,
   %    RESPONSE is a double array of size [NY NU NF] where RESPONSE(i,j,k) 
   %    specifies the frequency response from input j to output i at the 
   %    frequency point FREQS(k).
   %
   %    By default, the frequency points FREQS are assumed to be in radians/second.
   %    If FREQS is expressed in other units, set the "FrequencyUnit" and
   %    "TimeUnit" properties accordingly. For example, if FREQS is expressed
   %    in Hz, set FrequencyUnit='Hz'. If FREQS is expressed in radians/minute,
   %    set FrequencyUnit='rad/TimeUnit' and TimeUnit='minutes'.
   %
   %  Arrays of FRD models:
   %    You can create arrays of FRD models by using an ND array for RESPONSE.
   %    For example, if RESPONSE is an array of size [NY NU NF 3 4], then
   %       SYS = FRD(RESPONSE,FREQS)
   %    creates the 3-by-4 array of FRD models, where
   %       SYS(:,:,k,m) = FRD(RESPONSE(:,:,:,k,m),FREQS),  k=1:3,  m=1:4.
   %    Each of these FRD models has NY outputs, NU inputs, and data at
   %    the frequencies FREQS.
   %
   %  Conversion:
   %    SYS = FRD(SYS,FREQS,UNIT) converts any dynamic system SYS to the FRD
   %    representation by computing the system response at each frequency
   %    point in the vector FREQS.  The frequencies FREQS are expressed in
   %    the unit specified by the string UNIT (see "help frd.FrequencyUnit" for
   %    a list of frequency units). When UNIT is omitted, the default is 
   %    'rad/TimeUnit' (relative to the time units specified in SYS.TimeUnit).  
   %    The resulting SYS is of class @frd.
   %
   %  See also FRDATA, CHGFREQUNIT, FRDFUN, TF, ZPK, SS, DYNAMICSYSTEM.
   
%   Author(s): P. Gahinet
%   Copyright 1986-2011 The MathWorks, Inc.
   
   % Add static method to be included for compiler
   %#function ltipack.utValidateTs
   %#function frd.loadobj

   % Public properties with restricted value
   properties (Access = public, Dependent)
      % Frequency response data (numeric array).
      %
      % The "ResponseData" property stores the frequency response data as
      % a 3D array of complex numbers. For SISO systems, set "ResponseData" 
      % to the vector of frequency response values at the specified 
      % frequency points (see "Frequency" property). For MIMO systems with 
      % NU inputs and NY outputs, set "ResponseData" to an array of size 
      % [NY NU NW] where NW is the number of frequency points (the frequency 
      % runs along the third dimension). 
      %
      % For example,
      %    w = logspace(-2,2,50)
      %    h = freqresp(rss(5,2,3),w)  % h is 2x3x50
      %    sys = frd(h,w)   
      % creates an FRD model with 3 inputs, 2 outputs, and 50 frequency 
      % points. 
      ResponseData
      % Transport delays (numeric array, default = all zeros).
      %
      % The "IODelay" property specifies a separate time delay for each 
      % I/O pair. For continuous-time systems, specify I/O delays in
      % the time unit stored in the "TimeUnit" property. For discrete-
      % time systems, specify I/O delays as integer multiples of the
      % sampling period "Ts", for example, IODelay=2 to mean a delay
      % of two sampling periods.
      %
      % For MIMO systems with Ny outputs and Nu inputs, set "IODelay" to
      % a Ny-by-Nu matrix. You can also set "IODelay" to a scalar value to 
      % apply the same delay to all I/O pairs.
      %
      % While delays can be specified as part of the frequency response 
      % data ("ResponseData" property), keeping them separate results in 
      % more accurate phase plots by preventing drifts due to phase 
      % wrapping and the limited frequency grid resolution.
      IODelay
   end

   % OBSOLETE PROPERTIES
   properties (Access = public, Dependent, Hidden)
      % Shortened to IODelay
      ioDelayMatrix
      % Renames to Unit
      Units
   end
   
   % TYPE MANAGEMENT IN BINARY OPERATIONS
   methods (Static, Hidden)
           
      function T = toClosed(~)
         T = 'frd';
      end

      function T = superiorTypes()
         T = {'frd'};
      end

      function A = getAttributes(A)
         % Override default attributes
         A.Varying = false;
         A.Structured = false;
         A.Sparse = false;
      end
      
      function T = toVarying()
         error(message('Control:combination:FRD'))
      end

      function T = toStructured()
         % Use virtual class to ensure op(frd,ufrd)=ufrd
         T = 'neutral_genfrd';
      end
      
   end
   
   methods
      
      function sys = frd(varargin)
          
         % String support, revisit
         varargin = controllib.internal.util.hString2Char(varargin);
            
         ni = nargin;
         
         % Handle conversion FRD(SYS,F,U) where SYS is a FRD object
         if ni>0 && strcmp(class(varargin{1}),'frd') %#ok<STISA>
            if ni==1
               % no-op
               sys = varargin{1};
            else
               try
                  [sys,freq,funit] = FRDModel.parseFRDInputs('frd',varargin);
                  freq = funitconv(funit,'rad/TimeUnit',sys.TimeUnit)*freq;
                  % Note: May error if FREQ differs from frequency grid for some models
                  sys = copyMetaData(sys,frd_(sys,freq));  % frequnit=rad/TimeUnit
                  % Restore frequency units
                  sys = chgFreqUnit(sys,funit);
               catch ME
                  throw(ME)
               end
            end
            return
         end
         
         % Dissect input list
         DataInputs = 0;
         LtiInput = 0;
         PVStart = ni+1;
         for ct=1:ni
            nextarg = varargin{ct};
            if isa(nextarg,'struct') || isa(nextarg,'lti')
               LtiInput = ct;   PVStart = ct+1;   break
            elseif ischar(nextarg)
               PVStart = ct;   break
            else
               DataInputs = DataInputs+1;
            end
         end
         
         % Handle bad calls
         if PVStart==1
            if ni==1
               % Bad conversion
               error(message('Control:ltiobject:construct3','frd'))
            elseif ni>0
               error(message('Control:ltiobject:construct3','frd'))
            end
         elseif DataInputs>3
            error(message('Control:general:InvalidSyntaxForCommand','frd','frd'))
         end
         
         % Process numerical data
         try
            switch DataInputs
               case 0
                  if ni
                     error(message('Control:ltiobject:construct4','frd'))
                  else
                     freq = zeros(0,1);  resp = zeros(0,0,0);
                  end
               case 1
                  error(message('Control:ltiobject:frd2'))
               otherwise
                  % FREQS and RESPONSE specified
                  freq = FRDModel.checkFrequencyVector(varargin{2});
                  resp = FRDModel.checkResponseData(varargin{1},numel(freq));
            end
            
            % Sample time
            if DataInputs==3
               % Discrete SS
               Ts = ltipack.utValidateTs(varargin{3});
            else
               Ts = 0;
            end
         catch ME
            throw(ME)
         end
         
         % Determine I/O and array size
         if ni>0
            rsize = [size(resp) 1 1 1];
            Ny = rsize(1);
            Nu = rsize(2);
            ArraySize = rsize(4:end);
         else
            Ny = 0;  Nu = 0;  ArraySize = [1 1];
         end
         Nsys = prod(ArraySize);
         sys.IOSize_ = [Ny Nu];
         
         % Create @frddata object array
         % RE: Inlined for optimal speed
         if Nsys==1
            Data = ltipack.frddata(resp,freq,Ts);
         else
            Data = ltipack.frddata.array(ArraySize);
            Delay = ltipack.utDelayStruct(Ny,Nu,false);
            for ct=1:Nsys
               Data(ct).Frequency = freq;
               Data(ct).Response = resp(:,:,:,min(ct,end));
               Data(ct).Ts = Ts;
               Data(ct).Delay = Delay;
            end
         end
         sys.Data_ = Data;
         
         % Process additional settings and validate system
         % Note: Skip when just constructing empty instance for efficiency
         if ni>0
            try
               Settings = cell(1,0);
               
               % Properties inherited from other system
               if LtiInput
                  arg = varargin{LtiInput};
                  if isa(arg,'lti')
                     arg = getSettings(arg);
                  end
                  Settings = [Settings , lti.struct2pv(arg)];
               end
               
               % User-defined properties
               Settings = [Settings , varargin(:,PVStart:ni)];
               
               % Apply settings
               if ~isempty(Settings)
                  sys = fastSet(sys,Settings{:});
               end
               
               % Consistency check
               sys = checkConsistency(sys);
            catch ME
               throw(ME)
            end
         end
      end
      
      function Value = get.ResponseData(sys)
         % GET method for ResponseData property
         Data = sys.Data_;
         if isempty(Data)
            Value = zeros([sys.IOSize_,0,size(Data)]);
         else
            Value = zeros([sys.IOSize_,length(Data(1).Frequency),size(Data)]);
         end
         for ct=1:numel(Data)
            Value(:,:,:,ct) = Data(ct).Response;
         end
      end
            
      function sys = set.ResponseData(sys,Value)
         % SET method for ResponseData property
         % Validate type and add missing leading singleton dimensions
         % (needed to push array dimensions to 4 and above)
         Value = FRDModel.checkResponseData(Value,numel(sys.Frequency));
         % Check compatibility with model array sizes
         Data = ltipack.utCheckAssignValueSize(sys.Data_,Value,3);
         sv = size(Value);
         if isequal(sv(1:2),sys.IOSize_)
            % No change in I/O size
            for ct=1:numel(Data)
               Data(ct).Response = Value(:,:,:,min(ct,end));
               if sys.CrossValidation_
                  Data(ct) = checkData(Data(ct));  % Quick validation
               end
            end
            sys.Data_ = Data;
         else
            % I/O size changes
            for ct=1:numel(Data)
               Data(ct).Response = Value(:,:,:,min(ct,end));
            end
            sys.Data_ = Data;
            if sys.CrossValidation_
               % Note: Full validation needed because a single assignment can change I/O size,
               % e.g., sys = frd(1:10,1:10); sys.resp = randn(2,2,10);
               sys = checkConsistency(sys);
            end
         end
      end
            
      function Value = get.Units(sys)
         Value = sys.FrequencyUnit;
      end      
      function sys = set.Units(sys,Value)
         sys.FrequencyUnit = Value;
      end
      
      function Value = get.IODelay(sys)
         % GET method for IODelay property
         Value = getIODelay(sys);
      end
      
      function sys = set.IODelay(sys,Value)
         % SET method for IODelay property
         sys = setIODelay(sys,Value);
      end
      
      function Value = get.ioDelayMatrix(sys)
         Value = getIODelay(sys);
      end     
      function sys = set.ioDelayMatrix(sys,Value)
         sys = setIODelay(sys,Value);
      end
      
   end
   
   %% ABSTRACT SUPERCLASS INTERFACES
   methods (Access=protected)

      % INPUTOUTPUTMODEL
      function displaySize(sys,sizes)
         % Displays SIZE information in SIZE(SYS)
         ny = sizes(1);
         nu = sizes(2);
         nf = nfreqs(sys);
         if length(sizes)==2
            disp(getString(message('Control:ltiobject:SizeFRD1',ny,nu,nf)))
         else
            ArrayDims = sprintf('%dx',sizes(3:end));
            disp(getString(message('Control:ltiobject:SizeFRD2',...
               ArrayDims(1:end-1),ny,nu,nf)))
         end
      end
      
      % FRDMODEL
      function nf = nfreqs_(sys)
         % Returns number of frequency points
         D = sys.Data_;
         if isempty(D)
            nf = 0;
         else
            nf = numel(D(1).Frequency);
         end
      end
      
      function f = getFrequency_(sys)
         % Retrieves frequency vector expressed in current frequency units
         D = sys.Data_;
         f = D(1).Frequency;  % in rad/TimeUnit
         for ct=2:numel(D)
            if ~FRDModel.isSameFrequencyGrid(f,D(ct).Frequency)
               error(message('Control:ltiobject:get7'))
            end
         end
         % Return values in FrequencyUnit
         f = funitconv('rad/TimeUnit',sys.FrequencyUnit,sys.TimeUnit) * f;
      end
      
      function sys = setFrequency_(sys,f)
         % Sets frequency vector (given in the units specified by FrequencyUnit)
         f = funitconv(sys.FrequencyUnit,'rad/TimeUnit',sys.TimeUnit) * f;
         Data = sys.Data_;
         for ct=1:numel(Data)
            Data(ct).Frequency = f;
            if sys.CrossValidation_
               Data(ct) = checkData(Data(ct));
            end
         end
         sys.Data_ = Data;
      end
      
      function sys = scaleFrequency_(sys,sf)
         % Scales frequency vector (typically in reaction to a unit change)
         Data = sys.Data_;
         for ct=1:numel(Data)
            Data(ct).Frequency = sf * Data(ct).Frequency;
         end
         sys.Data_ = Data;
      end
      
      function sys = fcat_(sys,sys2)
         % FCAT(SYS,SYS2) for two FRD systems.
         [sys,sys2] = matchArraySize(sys,sys2);   % must come first
         [sys,sys2] = matchSamplingTime(sys,sys2);
         % Combine data
         sys.Data_ = fcat(sys.Data_,sys2.Data_);
      end
            
      function sys = fselect_(sys,varargin)
         % Select portion of frequency vector in FRD model.
         % FSELECT_(SYS,INDEX) or FSELECT_(SYS,FMIN,FMAX)
         ni = nargin;
         Data = sys.Data_;
         for ct=1:numel(Data)
            f = Data(ct).Frequency;
            if ni<3
               index = varargin{1};
               if any(index>length(f))
                  error(message('Control:transformation:fselect1'))
               end
            else
               % FMIN and FMAX are in rad/TimeUnit 
               index = find(f>=varargin{1} & f<=varargin{2});
            end
            Data(ct).Response = Data(ct).Response(:,:,index);
            Data(ct).Frequency = f(index,:);
         end
         sys.Data_ = Data;
      end
      
      function sys = fdel_(sys,freq2remove)
         % Delete portion of frequency vector in FRD model.
         % FREQ2REMOVE is specified in rad/TimeUnit.
         Data = sys.Data_;
         for ct=1:numel(Data)
            Data(ct) = fdel(Data(ct),freq2remove);
         end
         sys.Data_ = Data;
      end
      
      function sys = frdfun_(fun,sys)
         % Applies function to response data
         D = sys.Data_;
         if isempty(D)
            ios = size(feval(fun,zeros(iosize(sys))));
         else
            for ct=1:numel(D)
               D(ct) = transformResp(D(ct),fun);
            end
            ios = iosize(D(1));
         end
         sys.Data_ = D;
         sys.IOSize_ = ios;
      end
            
      % SINGLERATESYSTEM
      function sys = setTs_(sys,Ts)
         % Implementation of @SingleRateSystem:setTs_
         if Ts==-1
            warning(message('Control:ltiobject:frdAmbiguousRate1'))
         end            
         sys = setTs_@lti(sys,abs(Ts));
      end
      
   end
   
   
   %% DATA ABSTRACTION INTERFACE
   methods (Access=protected)
      
      %% MODEL CHARACTERISTICS
      function sys = checkDataConsistency(sys)
         % Cross validation of system data. Extends @lti implementation
         sys = checkDataConsistency@lti(sys);         
         % Sampling time restriction
         if sys.Ts==-1
            % Ts=-1 is ambiguous for FRD models and may lead to
            % inconsistencies, e.g., if sys1.Ts=-1 and sys2.Ts=.1,
            % frd(sys1,w)+frd(sys2,w) and frd(sys1+sys2,w) differ
            % because the response in frd(sys1,w) is effectively
            % evaluated for Ts=1. Similar problems arise when
            % absorbing delays with Ts=-1. Force Ts=1 and warn.
            warning(message('Control:ltiobject:frdAmbiguousRate1'))
            sys.Ts = 1;
         end
      end
      
      function boo = isempty_(sys)
         % Note: Returns true if there are no frequencies
         boo = any(size(sys)==0);
         if ~boo
            boo = true;
            D = sys.Data_;
            for ct=1:numel(D)
               if ~isempty(D(ct).Frequency)
                  boo = false;  return
               end
            end
         end
      end

      
      %% BINARY OPERATIONS
      function [sys1,SingularFlag] = mldivide_(sys1,sys2)
         % SYS1\SYS2 for two FRD systems.
         [sys1,sys2] = matchArraySize(sys1,sys2);   % must come first
         [sys1,sys2] = matchAttributes(sys1,sys2);  % overloadable
         % Combine data
         Data1 = sys1.Data_;  Data2 = sys2.Data_;  SingularFlag = false;
         for ct=1:numel(Data1)
            [Data1(ct),sf] = mldivide(Data1(ct),Data2(ct));
            SingularFlag = SingularFlag || sf;
         end
         sys1.Data_ = Data1;
         sys1.SamplingGrid_ = ...
            ltipack.SamplingGrid.merge(sys1.SamplingGrid_,sys2.SamplingGrid_);
         sys1.IOSize_ = sys2.IOSize_;
      end

      function [sys1,SingularFlag] = mrdivide_(sys1,sys2)
         % SYS1/SYS2 for two FRD systems.
         [sys1,sys2] = matchArraySize(sys1,sys2);   % must come first
         [sys1,sys2] = matchAttributes(sys1,sys2);  % overloadable
         % Combine data
         Data1 = sys1.Data_;  Data2 = sys2.Data_;  SingularFlag = false;
         for ct=1:numel(Data1)
            [Data1(ct),sf] = mrdivide(Data1(ct),Data2(ct));
            SingularFlag = SingularFlag || sf;
         end
         sys1.Data_ = Data1;
         sys1.SamplingGrid_ = ...
            ltipack.SamplingGrid.merge(sys1.SamplingGrid_,sys2.SamplingGrid_);
      end
      
      function [sys1,sys2] = matchAttributes(sys1,sys2)
         % Enforces matching attributes in binary operations (e.g.,
         % sample time, variable,...). 
         [sys1,sys2] = matchAttributes@lti(sys1,sys2);
         % Check that frequency vectors match
         D1 = sys1.Data_;   D2 = sys2.Data_;
         if ~(isempty(D1) || isempty(D2) || ...
               FRDModel.isSameFrequencyGrid(D1(1).Frequency,D2(1).Frequency))
            error(message('Control:ltiobject:mrgfreq1'))
         end
      end
      
      %% INDEXING
      function sys = indexasgn_(sys,indices,rhs,ioSize,ArrayMask)
         % Data management in SYS(indices) = RHS.
         % ioSize is the new I/O size and ArrayMask tracks which
         % entries in the resulting system array have been reassigned.
         Data = sys.Data_;
         % Construct template initial value for new entries in system array
         if isempty(Data)
            D0 = ltipack.frddata(zeros([ioSize 0]),zeros(0,1),0);
         else
            freqs = Data(1).Frequency;
            D0 = ltipack.frddata(zeros([ioSize length(freqs)]),freqs,Data(1).Ts);
            D0.FreqUnits = Data(1).FreqUnits;
         end
         % Update data
         try
            sys.Data_ = ltipack.reassignData(Data,indices,rhs.Data_,ioSize,ArrayMask,D0);
         catch ME
            if strcmp(ME.identifier,'Control:ltiobject:mrgfreq1')
               % Recast error message to assignment context
               error(message('Control:ltiobject:subsasgn4'))
            else
               throw(ME)
            end
         end
         % Update sampling grid
         if numel(indices)>2 && ~(isempty(sys.SamplingGrid_) && isempty(rhs.SamplingGrid_))
            sys.SamplingGrid_ = reassign(sys.SamplingGrid_,indices(3:end),rhs.SamplingGrid_,size(sys.Data_));
         end
      end
      
   end
      
   %% PROTECTED METHODS
   methods (Access=protected)
      
      function s = getPropStruct(sys)
         % Move "Frequency*" properties to the top
         s = getPropStruct@InputOutputModel(sys);
         n = numel(fieldnames(s));
         s = orderfields(s,[n-1 n 1:n-2]);
      end
      
      function sys = setTimeUnit_(sys,TU)
         % Changing TimeUnit value. 
         % When FrequencyUnit is not relative to TimeUnit, update stored 
         % frequency vector value to preserve user-facing value.
         if isempty(strfind(sys.FrequencyUnit,'TimeUnit'))
            cf = tunitconv(TU,sys.TimeUnit);
            Data = sys.Data_;
            for ct=1:numel(Data)
               Data(ct).Frequency = cf * Data(ct).Frequency;
            end
            sys.Data_ = Data;
         end
         % Update property value
         sys = setTimeUnit_@DynamicSystem(sys,TU);
      end
            
      function S = getSettings(sys)
         % Gets values of public LTI properties. Needed to support frd(R,w,LTI)
         S = getSettings@lti(sys);
         S.FrequencyUnit = sys.FrequencyUnit;
         S.IODelay = sys.IODelay;
      end
      
   end
   
   
   %% STATIC METHODS
   methods(Static, Hidden)
      
      sys = loadobj(s)
      
      function sys = make(D,IOSize)
         % Constructs FRD model from ltipack.frddata array
         sys = frd;
         sys.Data_ = D;
         if nargin>1
            sys.IOSize_ = IOSize;  % support for empty model arrays
         else
            sys.IOSize_ = iosize(D(1));
         end
      end
      
      function S = makeS(S,TU,FU)
         % Applies frd.make to each field of a struct array S.
         % Used by uncertainty analysis functions
         F = fieldnames(S);
         C = struct2cell(S);
         for ct=1:numel(C)
            V = C{ct};
            if isa(V,'ltipack.frddata')
               C{ct} = set(frd.make(V),'TimeUnit',TU,'FrequencyUnit',FU);
            end
         end
         S = cell2struct(C,F,1);
      end
      
      function sys = convert(X,refsys)
         % Safe conversion to FRD.
         %   X = FRD.CONVERT(X,REFSYS) casts the variable X to FRD using the
         %   frequency vector and frequency units of REFSYS. This method is  
         %   used in assignments, conversions to FRD, and binary operations to 
         %   correctly handle numeric arrays, static models, and undefined 
         %   sampling times.
         if nargin==1
            sys = frd(X);
         else
            f = refsys.Frequency;
            Ts = refsys.Ts;
            if isnumeric(X) || isa(X,'StaticModel')
               % Casting numeric array or static model to FRD
               s = size(X);
               if isnumeric(X)
                  X = full(double(X));
               else
                  X = double(X);
               end
               X = repmat(reshape(X,[s(1:2) 1 s(3:end)]),[1 1 length(f)]);
               sys = frd(X,f,Ts);
               sys.TimeUnit = refsys.TimeUnit;
               sys.FrequencyUnit = refsys.FrequencyUnit;
            elseif isa(X,'DynamicSystem')
               % Casting dynamic model to FRD
               if X.Ts==-1 && (Ts>0 || isstatic(X))
                  % Override Ts=-1 before computing frequency response
                  X.Ts = Ts;
               end
               sys = frd(X,f,refsys.FrequencyUnit);
            end
         end
      end
         
   end
   
end
