--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/m-toolbox/classes/@ao/fftfilt_core.m	Wed Nov 23 19:22:13 2011 +0100
@@ -0,0 +1,181 @@
+% 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