function [ decData ] = decimateByFactor(data,R,dimension,firSize)
%[ decData ] = decimateByFactor(data,R,dimension,firSize)
%
% Downsamples and applies anti-aliasing filtering, similar to the MATLAB 
% decimate function.  In contrast to the decimate function, however, 
% decimateByFactor first factors R and then performs downsampling/filtering 
% in stages from the largest to the smallest factor (using the MATLAB
% decimate command), thereby minimizing both computational demands and 
% aliasing artifacts.  The anti-aliasing filters use only FIR filters that 
% are selected based on the decimation rate or can be specified explicitly.
%
% Inputs: data        - Input data to decimate (can be vector or matrix)
%         R           - Total decimation factor
%         [dimension] - Dimension to operate on (default is 1st
%                       non-singleton)
%         [firSize]   - FIR filter size (default is based on data size) 
%
% Outputs: decData    - Decimated data
%
% Copyright (c) 2011, Hidden Solutions, LLC
% All rights reserved.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
%
%   * Redistributions of source code must retain the above copyright
%      notice, this list of conditions and the following disclaimer.
%    * Redistributions in binary form must reproduce the above copyright
%      notice, this list of conditions and the following disclaimer in the
%      documentation and/or other materials provided with the distribution.
%    * Neither the name Hidden Solutions, LLC nor the names of any
%      contributors may be used to endorse or promote products derived from 
%      this software without specific prior written permission.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
% ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
% WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
% DISCLAIMED. IN NO EVENT SHALL HIDDEN SOLUTIONS, LLC BE LIABLE FOR ANY
% DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
% (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
% LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
% ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
% (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
% SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
%
% Author: James S. Hall (Hidden Solutions, LLC)
% Date:   10 August 2011
%
% Example (cut and paste into the command line):
%
%   n              = 2e5;                         % Number of samples
%   D              = 60;                          % Decimation factor
%   L              = 10;                          % FIR filter length
%   a              = randn(1,n);                  % Define random set of data
%   R              = length(factor(D));           % Number of factors
%   tic; aDec      = D*decimate(a,D,L^R,'FIR');   % Decimate w/ decimate
%   time_decimate  = toc                          % Mark computation time
%   tic; aDecByF   = D*decimateByFactor(a,D,2,L); % Decimate w/ decimateByFactor
%   time_decimateByFactor  = toc                  % Mark comptation time
%   a_FFT          = fft(a);                      % FFT of original data
%   AD_FFT         = fft(aDec);                   % FFT of decimate results
%   ADBF_FFT       = fft(aDecByF);                % FFT of decimateByFactor results
%   indxLimit      = round(n/D/2);                % Number of freq samples to compare
%   % Measure default decimation error
%   error_decimate = norm(AD_FFT(1:indxLimit)   -  ...
%                         a_FFT(1:indxLimit))   ./ ...
%                    norm(a_FFT(1:indxLimit))           
%   % Measure decimation error with decimateByFactor
%   error_decimateByFactor = norm(ADBF_FFT(1:indxLimit) -  ...
%                                 a_FFT(1:indxLimit))   ./ ...
%                            norm(a_FFT(1:indxLimit))
%

% Force working dimension to 1st dimension
origDims = size(data);
if nargin < 3
  [data dimension] = shiftdim(data);
  dimension = dimension + 1;
else
  data = shiftdim(data,dimension-1);
end

% Set default firSize
if nargin < 4
  firSize = 0;
end

% Reshape into 2D array
data2D = reshape(data,size(data,1),[]);

% Find factors to decimate by
decFactors = sort(factor(R),'descend');

% Create output matrix
newData2D = zeros(ceil(size(data2D,1)/R),size(data2D,2));

% Step through each dataset
for cc = 1:size(data2D,2)
    
  % Assign filter length
  if firSize > 0
    filtLength = firSize;
  else
    filtLength = max(floor(size(data2D,1)/3),30);
  end
  
  % Perform first decimation
  if filtLength > size(data2D,1)/3
    tempData = decimate(data2D(:,cc),decFactors(1),'FIR');
  else
    tempData = decimate(data2D(:,cc),decFactors(1),filtLength,'FIR');
  end
  
  % Decimate by remaining factors
  for ff = 2:length(decFactors)
      
    % Assign filter length
    if firSize > 0
      filtLength = firSize;
    else
      filtLength = max(floor(size(tempData,1)/3),30);
    end
    
    % Perform additional decimation stages
    if filtLength > size(tempData,1)/3
      tempData = decimate(tempData,decFactors(ff),floor((length(tempData)-1)/2)-1,'FIR');
    else
      tempData = decimate(tempData,decFactors(ff),filtLength,'FIR');
    end
    
  end

  % Store results
  newData2D(:,cc) = tempData;
  
end

% Determine output dimensions
newDims            = origDims;
newDims(dimension) = size(newData2D,1);

% Reshape to match input format
newData = reshape(newData2D,circshift(newDims,[0 1-dimension]));
decData = permute(newData,circshift(1:length(size(newData)),[0 dimension-1]));
