% FFTFILT_CORE Simple core method which computes the fft filter.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% DESCRIPTION: Simple core method which computes the fft filter.
%
% CALL:        ao = fftfilt_core(ao, filt, Npad)
%
% INPUTS:      ao:   Single input analysis object
%              Npad: Number of bins for zero padding
%              filt: The filter to apply to the data
%                      smodel - a model to filter with.
%                      mfir   - an FIR filter
%                      miir   - an IIR filter
%                      tf     - an ltpda_tf object. Including:
%                                - pzmodel
%                                - rational
%                                - parfrac
%
% VERSION:     $Id: fftfilt_core.m,v 1.17 2011/05/30 10:42:32 mauro Exp $
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function bs = fftfilt_core(varargin)
  
  bs    =   varargin{1};
  filt  =   varargin{2};
  Npad  =   varargin{3};
  
  if nargin == 4
    inCondsMdl = varargin{4};
  else
    inCondsMdl = [];
  end
  
  [m, n] = size(bs.data.y);
  % FFT time-series data
  fs   = bs.data.fs;
  % zero padding data before fft
  if isempty(Npad)
    Npad = length(bs.data.y) - 1;
  end
  if n == 1
    tdat = ao(tsdata([bs.data.y;zeros(Npad,1)], fs));
  else
    tdat = ao(tsdata([bs.data.y zeros(1,Npad)], fs));
  end
  % get onesided fft
  ft = fft_core(tdat, 'one');
  
  switch class(filt)
    case 'smodel'
      % Evaluate model at the given frequencies
      
      amdly = filt.setXvals(ft.data.x).double;
      amdly = reshape(amdly, size(ft.data.y));
      
      amdl = ao(fsdata(ft.data.x, amdly, fs));
      
      % set units
      bs.setYunits(simplify(bs.data.yunits .* filt.yunits));
      
    case {'miir', 'mfir', 'pzmodel', 'parfrac', 'rational'}
      
      % Check if the frequency of the filter is the same as the frequency
      % of the AO.
      if isa (filt, 'ltpda_filter') && fs ~= filt.fs
        error('### Please use a filter with the same frequency as the AO [%dHz]', fs);
      end
      
      % get filter response on given frequencies
      amdl = resp(filt, plist('f', ft.data.x));
      
      % set units
      bs.setYunits(simplify(bs.data.yunits .* filt.ounits ./ filt.iunits));
      
    case 'ao'
      
      % check if filter and data have the same shape
      if size(ft.data.y)~=size(filt.data.y)
        % reshape
        amdl = copy(filt,1);
        amdl.setX(ft.data.x);
        amdl.setY(reshape(filt.data.y,size(ft.data.x)));
        amdl.setName(filt.name);
      else
        amdl = copy(filt,1);
        amdl.setName(filt.name);
      end
      
      % set units
      bs.setYunits(simplify(bs.data.yunits .* amdl.data.yunits));
      
    case 'collection'
      
      % run over collection elements
      amdl = ao(fsdata(ft.data.x, ones(size(ft.data.x)), fs));
      for ii=1:numel(filt.objs)
        switch class(filt.objs{ii})
          case 'smodel'
            % Evaluate model at the given frequencies
            amdly = filt.objs{ii}.setXvals(ft.data.x).double;
            amdly = reshape(amdly, size(ft.data.y));
            amdl_temp = ao(fsdata(ft.data.x, amdly, fs));
            % set units
            bs.setYunits(simplify(bs.data.yunits .* filt.objs{ii}.yunits));

          case {'miir'}
            % get filter response on given frequencies
            amdly = utils.math.mtxiirresp(filt.objs{ii},ft.data.x,fs,[]);
            amdl_temp = ao(fsdata(ft.data.x, amdly, fs));
            % set units
            bs.setYunits(simplify(bs.data.yunits .* filt.objs{ii}.ounits ./ filt.objs{ii}.iunits));
            
          case 'ao'
            % check if filter and data have the same shape
            if size(ft.data.y)~=size(filt.objs{ii}.data.y)
              % reshape
              amdl_temp = copy(filt.objs{ii},1);
              amdl_temp.setX(ft.data.x);
              amdl_temp.setY(reshape(filt.objs{ii}.data.y,size(ft.data.x)));
              amdl_temp.setName(filt.objs{ii}.name);
            else
              amdl_temp = copy(filt.objs{ii},1);
              amdl_temp.setName(filt.objs{ii}.name);
            end
            % set units
            bs.setYunits(simplify(bs.data.yunits .* amdl_temp.data.yunits));
            
          case 'filterbank'
            % get filter response on given frequencies
            amdly = utils.math.mtxiirresp(filt.objs{ii}.filters,ft.data.x,fs,filt.objs{ii}.type);
            amdl_temp = ao(fsdata(ft.data.x, amdly, fs));
            % handle units
            switch lower(filt.objs{ii}.type)
              case 'parallel'
                % set units of the output object
                bs.setYunits(simplify(bs.data.yunits .* filt.objs{ii}.filters(1).ounits ./ filt.objs{ii}.filters(1).iunits));
              case 'series'
                % get units from the series
                sunits = filt.objs{ii}.filters(1).ounits ./ filt.objs{ii}.filters(1).iunits;
                for jj = 2:numel(filt.objs{ii}.filters)
                  sunits = sunits.*filt.objs{ii}.filters(jj).ounits ./ filt.objs{ii}.filters(jj).iunits;
                end
                % set units of the output object
                bs.setYunits(simplify(bs.data.yunits .* sunits));
            end
        end
        % update response
        amdl = amdl .* amdl_temp;
      end
      
    otherwise
      
      error('### Unknown filter mode.');
      
  end
  
  % Add initial conditions
  if ~isempty(inCondsMdl) && ~isempty(inCondsMdl.expr.s)
    inCondsMdl.setXvals(amdl.x);
    inCondsMdl.setXunits(amdl.xunits);
    inCondsMdl.setYunits(amdl.yunits*ft.yunits);
    inCondsEval = inCondsMdl.eval;
    % Multiply by model and take inverse FFT
    y = ifft_core(ft.*amdl+inCondsEval, 'symmetric');
  else
    % Multiply by model and take inverse FFT
    y = ifft_core(ft.*amdl, 'symmetric');
  end
  
  % split and reshape the data
  if m == 1
    bs.setY(y.data.getY(1:n));
  else
    bs.setY(y.data.getY(1:m));
  end
  
  % clear errors
  bs.clearErrors;
  
end