function Y = goertzel(X,I,DIM)
%GOERTZEL Second-order Goertzel algorithm.  
%   GOERTZEL(X,I) computes the discrete Fourier transform (DFT) of 
%   X at indices contained in the vector I, using the second-order 
%   Goertzel algorithm.  The indices must be integer values from 1 
%   to N where N is the length of the first non-singleton dimension.  
%   If empty or omitted, I is assumed to be 1:N.
%
%   GOERTZEL(X,[],DIM) or GOERTZEL(X,I,DIM) computes the DFT along 
%   the dimension DIM.
%
%   In general, GOERTZEL is slower than FFT when computing all the possible 
%   DFT indices, but is most useful when X is a long vector and the DFT 
%   computation is required for only a subset of indices less than log2(length(X)). 
%   Indices 1:length(X) correspond to the frequency span [0, 2*pi) radians.
%
%   See also FFT, FFT2.

%   Author: P. Costa
%   Copyright 1988-2002 The MathWorks, Inc.
%   $Revision: 1.5 $  $Date: 2002/03/28 17:28:15 $

%   Reference:
%     C.S.Burrus and T.W.Parks, DFT/FFT and Convolution Algorithms, 
%     John Wiley & Sons, 1985

%if ~(exist('goertzelmex') == 3), 
%	error('Unable to find goertzelmex executable');
%end

%error(nargchk(1,3,nargin));
if nargin < 2, I = []; end
if nargin < 3, DIM = []; end

% Inputs should all be of class 'double'
%error(checkdatatype('goertzel','double',X,I,DIM));

if ~isempty(DIM) & DIM > ndims(X)
	error('Dimension specified exceeds the dimensions of X.')
end

% Reshape X into the right dimension.
if isempty(DIM)
	% Work along the first non-singleton dimension
	[X, nshifts] = shiftdim(X);
else
	% Put DIM in the first dimension (this matches the order 
	% that the built-in filter function uses)
	perm = [DIM,1:DIM-1,DIM+1:ndims(X)];
	X = permute(X,perm);
end

% Verify that the indices in I are valid.
siz = size(X);
if isempty(I),
	I = 1:siz(1); % siz(1) is the number of rows of X
else
	I = I(:);
	if max(I) > siz(1),
		error('Maximum Index exceeds the dimensions of X.');
	elseif min(I) < 1
		error('Index must greater than one.');
	elseif all(I-fix(I))
		error('Indices must be integer values.');
	end
end

% Initialize Y with the correct dimension
Y = zeros([length(I),siz(2:end)]); 

% Call goertzelmex 
for k = 1:prod(siz(2:end)),
	Y(:,k) = goertzelmex(X(:,k),I);
end

% Convert Y to the original shape of X
if isempty(DIM)
	Y = shiftdim(Y, -nshifts);
else
	Y = ipermute(Y,perm);
end


%-------------------------------------------------------------------
%                       Utility Function
%-------------------------------------------------------------------
function msg = checkdatatype(fname,cname,varargin)
%CHECKDATATYPE Check data type class
%   Checks for exact match on class, subclasses don't count.  i.e.,
%   sparse will not be a double, even though isa(sparse(2),'double')
%   evaluates to true.
%
% Inputs:
%   fname    - Function name
%   cname    - Class name
%   varargin - Inputs to check

msg = '';
for k=1:length(varargin),
  if ~strcmpi(class(varargin{k}),cname)
    msg = ['Function ''',fname,...
          ''' is not defined for variables of class ''',...
          class(varargin{k}),''''];
    break;
  end
end

% [EOF] goertzel.m
