classdef FRDModel < DynamicSystem
   % Frequency Response Data Model objects.
   %
   %   Frequency Response Data (FRD) models are frequency-domain models
   %   that describe the behavior of a dynamic system in terms of its
   %   input/output frequency response. FRD models are useful to manipulate
   %   experimental data collected from a frequency analyzer and you can use
   %   them as surrogates for TF or SS models when designing control systems
   %   in the frequency domain.
   %
   %   All FRD model types derive from the @FRDModel superclass. This class
   %   is not user-facing and cannot be instantiated. Use the FRD command to
   %   construct ordinary FRD models. Combining such models with Control Design
   %   blocks produces generalized FRD models (see GENFRD).
   %
   %   See also FRD, GENFRD, CONTROLDESIGNBLOCK.
   
   %   Copyright 1986-2014 The MathWorks, Inc.
   
   %   All FRD model types derive from the @FRDModel superclass. This class
   %   is not user-facing and cannot be instantiated. Use the @frd class to
   %   construct FRD models. Combining FRD models with Control Design blocks
   %   produces GENFRD or UFRD models.
   
   properties (Access = public, Dependent)
      % Vector of frequency points.
      %
      % This property stores the frequency points (vector of positive real numbers).
      % Frequency values should be expressed in the units specified by the
      % "FrequencyUnit" property.
      Frequency
      % Frequency units (default = rad/TimeUnit).
      %
      % This property specifies in what units the frequency vector is expressed (see
      % "Frequency" property). It can be set to any of the strings 'rad/TimeUnit',
      % 'cycles/TimeUnit', 'rad/s', 'Hz', 'kHz', 'MHz', 'GHz', or 'rpm'. The units
      % 'rad/TimeUnit' and 'cycles/TimeUnit' are relative to the time units specified
      % in the "TimeUnit" property (see help for DynamicSystem.TimeUnit property).
      FrequencyUnit;
   end

   properties (Hidden, Dependent)
      % Support common abbreviations of Frequency for backward compatibility
      % (treated as ambiguous since adopting RedefinesParen API)
      F
      Freq  
   end

   properties (Access = protected)
      % Compact storage of frequency units (default stored as [])
      FrequencyUnit_
   end
   
   %%%%%%%%%%%%%%% PROTECTED INTERFACE %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   
   % TO BE IMPLEMENTED BY ALL SUBCLASSES
   methods (Abstract, Access = protected)
      nf = nfreqs_(sys)
      f = getFrequency_(sys)
      sys = setFrequency_(sys,f)
      sys = scaleFrequency_(sys,sf)
      sys = fselect_(sys,index)
      sys = fdel_(sys,ind2remove)
      sys = fcat_(sys,sys2)
   end
   
   % TO BE OPTIONALLY IMPLEMENTED BY SUBCLASSES
   methods (Access = protected)
      sysfn = fnorm_(sys,ntype)
      
      function fsys = frdfun_(fun,sys)
         % Redirect to FRD implementation
         fsys = frdfun_(fun,frd(sys));
      end
   end
   
   %%%%%%%%%%%%%%%%% END PROTECTED INTERFACE %%%%%%%%%%%%%%%%%%%%%%
   
   
   % PUBLIC FRD METHODS
   methods
      
      function Value = get.Frequency(sys)
         % GET method for Frequency property
         if nmodels(sys)==0
            Value = zeros(0,1);
         else
            Value = getFrequency_(sys);
         end
      end
      
      function sys = set.Frequency(sys,Value)
         % SET method for Frequency property.
         try
            sys = setFrequency_(sys,FRDModel.checkFrequencyVector(Value));
         catch ME
            throw(ME)
         end
      end
      
      function Value = get.FrequencyUnit(sys)
         % GET method for FrequencyUnit property
         Value = sys.FrequencyUnit_;
         if isempty(Value)
            Value = 'rad/TimeUnit';  % default
         end
      end
      
      function sys = set.FrequencyUnit(sys,Value)
         % SET method for FrequencyUnit property
         Unit = ltipack.matchKey(Value,ltipack.getValidFrequencyUnits());
         if isempty(Unit)
            error(message('Control:ltiobject:setFRD1'))
         end
         % Update stored value of frequency vector so that value of SYS.Frequency
         % remains the same in the new units
         TU = sys.TimeUnit;
         sys = scaleFrequency_(sys,funitconv(Unit,sys.FrequencyUnit,TU));
         % Store new value
         if strcmp(Unit,'rad/TimeUnit')
            sys.FrequencyUnit_ = [];
         else
            sys.FrequencyUnit_ = Unit;
         end
      end
      
      function Value = get.F(sys)
         Value = sys.Frequency;
      end
      
      function sys = set.F(sys,Value)
         try 
            sys.Frequency = Value;
         catch ME
            throw(ME)
         end
      end

      function Value = get.Freq(sys)
         Value = sys.Frequency;
      end
      
      function sys = set.Freq(sys,Value)
         try 
            sys.Frequency = Value;
         catch ME
            throw(ME)
         end
      end

      function nf = nfreqs(sys)
         %NFREQS  Number of frequency points in FRD model.
         %
         %   NF = NFREQS(SYS) returns the number NF of frequency points in the
         %   FRD model SYS.
         %
         %   See also SIZE, FRD.
         nf = nfreqs_(sys);
      end
      
      function sys = chgFreqUnit(sys,newUnits)
         %CHGFREQUNIT  Changes frequency units in FRD model.
         %
         %   SYS = CHGFREQUNIT(SYS,NEWUNITS) changes the units of the frequency
         %   points in the FRD model SYS. The string NEWUNITS specifies the new
         %   units among the following: 'rad/TimeUnit', 'cycles/TimeUnit', 'rad/s',
         %   'Hz', 'kHz', 'MHz', 'GHz', and 'rpm'. The first two options are
         %   relative to the system's time units defined in SYS.TimeUnit (seconds
         %   by default).
         %
         %   CHGFREQUNIT adjusts the values of both Frequency and FrequencyUnit
         %   to maintain the same frequency response.
         %
         %   See also FRD, CHGTIMEUNIT.
         newUnits = ltipack.matchKey(newUnits,ltipack.getValidFrequencyUnits());
         if isempty(newUnits)
            error(message('Control:ltiobject:setFRD2'))
         end
         % Update storage value for frequency units
         % Note: No data change since frequency vector is always stored in rad/TimeUnit
         if strcmp(newUnits,'rad/TimeUnit')
            sys.FrequencyUnit_ = [];
         else
            sys.FrequencyUnit_ = newUnits;
         end
      end
      
      function sys = fdel(sys,freq2remove)
         %FDEL  Removes specific frequencies from FRD model.
         %
         %   SYS = FDEL(SYS,FREQS) removes the frequency points nearest to the values
         %   specified in the vector FREQS from the FRD model SYS. The output SYS is
         %   the FRD model containing the remaining frequency points. The frequencies
         %   FREQS should be expressed in the FRD model units (see "FrequencyUnit"
         %   property).
         %
         %   To select frequencies in a particular range, use FSELECT rather than FDEL.
         %
         %   See also FSELECT, FCAT, FRD.
         
         % Check number of input & output arguments
         narginchk(2,2)
         sys = fdel_(sys,funitconv(sys.FrequencyUnit,...
            'rad/TimeUnit',sys.TimeUnit) * freq2remove);
      end
      
      function sys = abs(sys)
         %ABS  Frequency response magnitude for FRD models.
         %
         %   ABSFRD = ABS(SYS) computes the magnitude of the frequency
         %   response contained in the FRD model SYS. For MIMO models,
         %   the magnitude is computed for each I/O pair. The output ABSFRD
         %   is an FRD model containing the magnitude data across
         %   frequencies.
         %
         %   See also BODEMAG, SIGMA, FNORM, FRDMODEL/ANGLE, FRDMODEL/REAL,
         %   FRDMODEL/IMAG, FRDFUN.
         sys = frdfun(@abs,sys);
      end
      
      function sys = angle(sys)
         %ANGLE  Frequency response phase for FRD models.
         %
         %   ANGFRD = ANGLE(SYS) computes the phase in radians of the
         %   frequency response contained in the FRD model SYS. For MIMO
         %   models, the phase is computed for each I/O pair. The output
         %   ANGFRD is an FRD model containing the phase data across
         %   frequencies.
         %
         %   See also BODE, FRDMODEL/ABS, FRDMODEL/REAL, FRDMODEL/IMAG, FRDFUN.
         sys = frdfun(@angle,sys);
      end
      
      function sys = real(sys)
         %REAL  Real part of frequency response for FRD models.
         %
         %   REALFRD = REAL(SYS) computes the real part of the frequency
         %   response contained in the FRD model SYS, including the
         %   contribution of input, output, and I/O delays.  The output
         %   REALFRD is an FRD object containing the values of the real
         %   part across frequencies.
         %
         %   See also NYQUIST, FRDMODEL/IMAG, FRDMODEL/ABS, FRDMODEL/ANGLE, FRDFUN.
         sys = frdfun(@real,sys);
      end
      
      function sys = imag(sys)
         %IMAG  Imaginary part of frequency response for FRD models.
         %
         %   IMAGFRD = IMAG(SYS) computes the imaginary part of the
         %   frequency response contained in the FRD model SYS, including
         %   the contribution of input, output, and I/O delays.  The output
         %   IMAGFRD is an FRD object containing the values of the imaginary
         %   part across frequencies.
         %
         %   See also NYQUIST, FRDMODEL/REAL, FRDMODEL/ABS, FRDMODEL/ANGLE, FRDFUN.
         sys = frdfun(@imag,sys);
      end
      
      function sysfn = fnorm(sys, ntype)
         %FNORM  Pointwise gain of FRD model.
         %
         %   FNRM = FNORM(SYS) takes an FRD model SYS and computes the gain of
         %   the frequency response at each frequency point. For MIMO models,
         %   "gain" refers to the largest singular value of the frequency
         %   response matrix. The output FNRM is an FRD object containing the
         %   gain values across frequencies.
         %
         %   FNRM = FNORM(SYS,NTYPE) computes the frequency response gains using
         %   the matrix norm specified by NTYPE. See NORM for valid matrix norms
         %   and corresponding NTYPE values. The 2-norm (largest singular value)
         %   is used by default.
         %
         %   See also DYNAMICSYSTEM/NORM, FRDMODEL/ABS, FRDFUN.
         if nargin<2
            ntype = 2;
         elseif ~(isequal(ntype,1) || isequal(ntype,2) || isequal(ntype,Inf) || ...
               isequal(ntype,'fro') || isequal(ntype,"fro"))
            error(message('Control:analysis:fnorm1'))
         end
         sysfn = frdfun(@(x) norm(x,ntype),sys);
      end
      
      function sys = diag(sys)
         %DIAG  Diagonal frequency response or diagonal of a frequency response.
         %
         %   MSYS = DIAG(VSYS) takes an FRD model VSYS with vector-valued frequency
         %   response V1,...,Vn and constructs the FRD model MSYS with matrix-valued
         %   response DIAG(V1),...,DIAG(Vn).
         %
         %   Conversely, VSYS = DIAG(MSYS) takes a matrix-valued frequency response
         %   and returns the vector of diagonal elements at each frequency.
         %
         %   This operation is equivalent to applying DIAG to the response data
         %   at each frequency point.
         %
         %   See also APPEND, FRD.
         try
            sys = frdfun(@diag,sys);
         catch ME
            throw(ME)
         end
      end
      
      function sys = sqrtm(sys)
         %SQRTM  Takes square root of frequency response.
         %
         %   SYS2 = SQRTM(SYS) takes an FRD model SYS and returns an FRD
         %   model SYS2 whose frequency response is the matrix square root
         %   of the frequency response of SYS (see SQRTM).
         %
         %   See also SQRTM, FRD.
         ios = iosize(sys);
         if ios(1)~=ios(2)
            error(message('Control:transformation:frdfun5'))
         end
         try
            sys = frdfun(@sqrtm,sys);
         catch ME
            throw(ME)
         end
      end
      
   end
   
   %% OBSOLETE
   methods (Hidden)
      
      function sys = chgunits(sys,newUnits)
         %CHGUNITS  Change frequency unit in FRD models.
         %
         %   SYS = CHGUNITS(SYS,UNIT) changes the unit of the frequency
         %   points stored in the FRD model SYS to UNIT, where UNIT
         %   is either 'Hz or 'rad/s'.  A 2*pi scaling factor is applied
         %   to the frequency values and the 'FrequencyUnit' property is
         %   updated.
         %
         %   See also FRD.
         sys = chgFreqUnit(sys,newUnits);
      end
      
   end
   
   
   methods (Static, Hidden)
      
      [Model,w,unit] = parseFRDInputs(FRDfcn,InputList)
      
      function boo = isSameFrequencyGrid(freqs1,freqs2)
         % Checks that two frequency grids are equal within numerical roundoff
         reltol = 1e3*eps;
         boo = (numel(freqs1)==numel(freqs2) && ...
            all(abs(freqs1-freqs2)<=reltol*(1+abs(freqs1))));
      end
      
      function f = checkFrequencyVector(f)
         % Frequency vector validation
         if ~isempty(f)
            if ~(isvector(f) && isnumeric(f) && isreal(f) && allfinite(f))
               error(message('Control:ltiobject:frdProperties1'))
            elseif any(diff(sort(f))==0)
               error(message('Control:ltiobject:frdProperties5'))
            end
         end
         f = double(full(f(:)));
      end
      
      function A = checkResponseData(A,Nf)
         % Response data validation
         if isnumeric(A)
            A = double(full(A));
         else
            error(message('Control:ltiobject:frdProperties2'))
         end
         % Adjust size
         if isequal(A,[])
            % Support frd([],w) as shorthand for frd(zeros(0,0,nf),w)
            A = zeros(0,0,Nf);
         else
            % Fill in missing leading singleton dimensions
            rsize = [size(A) 1];
            if Nf~=rsize(3)
               % Look for missing leading singleton dimensions
               if Nf==rsize(1)
                  rsize = [1 1 rsize];
               elseif rsize(1)==1 && Nf==rsize(2)
                  rsize = [1 rsize];
               end
               A = reshape(A,rsize); % g347287
            end
         end
      end
      
   end
   
   
end
