Using CZT for zoom-FFT and zoom-IFFT in Matlab/Octave

czt.m

%  Copyright (C) 2004 Daniel Gunyan
% 
%  This program is free software; you can redistribute it and/or modify it under
%  the terms of the GNU General Public License as published by the Free Software
%  Foundation; either version 3 of the License, or (at your option) any later
%  version.
% 
%  This program is distributed in the hope that it will be useful, but WITHOUT
%  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
%  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
%  details.
% 
%  You should have received a copy of the GNU General Public License along with
%  this program; if not, see <http://www.gnu.org/licenses/>.

%  -*- texinfo -*-
%  @deftypefn  {Function File} {} czt (@var{x})
%  @deftypefnx {Function File} {} czt (@var{x}, @var{m})
%  @deftypefnx {Function File} {} czt (@var{x}, @var{m}, @var{w})
%  @deftypefnx {Function File} {} czt (@var{x}, @var{m}, @var{w}, @var{a})
%  Chirp z-transform.  Compute the frequency response starting at a and
%  stepping by w for m steps.  a is a point in the complex plane, and
%  w is the ratio between points in each step (i.e., radius increases
%  exponentially, and angle increases linearly).
% 
%  To evaluate the frequency response for the range f1 to f2 in a signal
%  with sampling frequency Fs, use the following:
% 
%  @example
%  @group
%  m = 32;                               %% number of points desired
%  w = exp(-j*2*pi*(f2-f1)/((m-1)*Fs));  %% freq. step of f2-f1/m
%  a = exp(j*2*pi*f1/Fs);                %% starting at frequency f1
%  y = czt(x, m, w, a);
%  @end group
%  @end example
% 
%  If you don't specify them, then the parameters default to a Fourier
%  transform:
%      m=length(x), w=exp(-j*2*pi/m), a=1
% 
%  If x is a matrix, the transform will be performed column-by-column.
%  @end deftypefn

%  Algorithm (based on Oppenheim and Schafer, 'Discrete-Time Signal
%  Processing', pp. 623-628):
%    make chirp of length -N+1 to max(N-1,M-1)
%      chirp => w^([-N+1:max(N-1,M-1)]^2/2)
%    multiply x by chirped a and by N-elements of chirp, and call it g
%    convolve g with inverse chirp, and call it gg
%      pad ffts so that multiplication works
%      ifft(fft(g)*fft(1/chirp))
%    multiply gg by M-elements of chirp and call it done

function y = czt(x, m, w, a)

  if nargin < 1 || nargin > 4, print_usage; end

  [row, col] = size(x);
  if row == 1, x = x(:); col = 1; end

  if nargin < 2 || isempty(m), m = length(x(:,1)); end
  if length(m) > 1, error('czt: m must be a single element\n'); end
  if nargin < 3 || isempty(w), w = exp(-2*j*pi/m); end
  if nargin < 4 || isempty(a), a = 1; end
  if length(w) > 1, error('czt: w must be a single element\n'); end
  if length(a) > 1, error('czt: a must be a single element\n'); end

  % indexing to make the statements a little more compact
  n = length(x(:,1));
  N = [0:n-1]'+n;
  NM = [-(n-1):(m-1)]'+n;
  M = [0:m-1]'+n;

  nfft = 2^nextpow2(n+m-1); % fft pad
  W2 = w.^(([-(n-1):max(m-1,n-1)]'.^2)/2); % chirp

  for idx = 1:col
    fg = fft(x(:,idx).*(a.^-(N-n)).*W2(N), nfft);
    fw = fft(1./W2(NM), nfft);
    gg = ifft(fg.*fw, nfft);

    y(:,idx) = gg(M).*W2(M);
  end

  if row == 1, y = y.'; end

end

%!shared x
%! x = [1,2,4,1,2,3,5,2,3,5,6,7,8,4,3,6,3,2,5,1];
%!assert(fft(x),czt(x),10000*eps);
%!assert(fft(x'),czt(x'),10000*eps);
%!assert(fft([x',x']),czt([x',x']),10000*eps);

czttest.m

clear

N = 10;
idxOfInterest = linspace(4, 5, 13); % Indexes of DFF desired
idxOfInterest = -linspace(-1, 1, 13); % both of these work

forward = diff(idxOfInterest(1:2)) > 0;

x = randn(N, 1) + 1j * randn(N, 1);
F = exp(-1j * 2 * pi / N * idxOfInterest(:) * [0:N-1]);
y = F * x;

% check F:
relerr = @(dirt, gold) abs((dirt(:) - gold(:)) ./ gold(:));
if 1
  if forward
    fx = fft(x);
    [~, fxIdx, yIdx] = intersect(0:N, idxOfInterest);
  else
    fx = ifft(x) * N;
    [~, fxIdx, yIdx] = intersect(0:N, mod(-idxOfInterest, N));
  end
  max(relerr(fx(fxIdx), y(yIdx)))
end

% Generate y via CZT
a = exp(1j * 2 * pi * idxOfInterest(1) / N);
if length(idxOfInterest) >= 2
  w = exp(-1j * 2 * pi / N * diff(idxOfInterest(1:2)));
else
  w = 1; % doesn't matter
end
n = length(idxOfInterest);
yz = czt(x, n, w, a);

% Check yz
max(relerr(yz, y))