www.gusucode.com > signal 工具箱matlab源码程序 > signal/goertzel.m

    function Y = goertzel(X,INDVEC,DIM)
%GOERTZEL Second-order Goertzel algorithm.  
%   GOERTZEL(X,INDVEC) computes the discrete Fourier transform (DFT)
%   of X at indices contained in the vector INDVEC, 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, INDVEC is assumed to be 1:N.
%
%   GOERTZEL(X,[],DIM) or GOERTZEL(X,INDVEC,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.
%
%   EXAMPLE:
%      % Resolve the 1.24 kHz and 1.26 kHz components in the following
%      % noisy cosine which also has a 10 kHz component.
%      Fs = 32e3;   t = 0:1/Fs:2.96;
%      x  = cos(2*pi*t*10e3)+cos(2*pi*t*1.24e3)+cos(2*pi*t*1.26e3)...
%           + randn(size(t));
%
%      N = (length(x)+1)/2;
%      f = (Fs/2)/N*(0:N-1);              % Generate frequency vector
%      indxs = find(f>1.2e3 & f<1.3e3);   % Find frequencies of interest
%      X = goertzel(x,indxs);
%      
%      plot(f(indxs)/1e3,20*log10(abs(X)/length(X)));
%      title('Mean Squared Spectrum');
%      xlabel('Frequency (kHz)');
%      ylabel('Power (dB)');
%      grid on;
%      set(gca,'XLim',[f(indxs(1)) f(indxs(end))]/1e3);
%
%   See also FFT, FFT2.

%   Copyright 1988-2013 The MathWorks, Inc.

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

if ~(exist('goertzelmex', 'file') == 3), 
	error(message('signal:goertzel:NotSupported'));
end

narginchk(1,3);
if nargin < 2, INDVEC = []; end
if nargin < 3, DIM = []; end

% Cache input data type to cast output at end of computations
isDataSingle = false;
if isa(X,'single')
  isDataSingle = true;
end
X = signal.internal.sigcasttofloat(X,'double','goertzel','X');
INDVEC = signal.internal.sigcasttofloat(INDVEC,'double','goertzel',...
  'INDVEC','allownumeric');
DIM = signal.internal.sigcasttofloat(DIM,'double','goertzel','DIM',...
  'allownumeric');

if ~isempty(DIM) && DIM > ndims(X)
	error(message('signal:goertzel:InvalidDimensions'))
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 INDVEC are valid.
siz = size(X);
if isempty(INDVEC),
	INDVEC = 1:siz(1); % siz(1) is the number of rows of X
else
	INDVEC = INDVEC(:);
	if max(INDVEC) > siz(1),
		error(message('signal:goertzel:IdxGtBound'));
	elseif min(INDVEC) < 1
		error(message('signal:goertzel:IdxLtBound'));
	elseif all(INDVEC-fix(INDVEC))
		error(message('signal:goertzel:MustBeInteger'));
	end
end

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

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

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

% If X is single, then relevant output should be single
if isDataSingle
  Y = single(Y);
end

% [EOF] goertzel.m