FFT for arbitrary time and frequency grids (Chirp-z Transform)

The FFT is usable only for specific values of the time and frequency spacings;  (their
product must be the inverse of a highly composite integer).  In practice, however, we find
arbitrary spacings, and thus  the FFT is not directly applicable. The chirp-z  transform
allows one to compute the DFT  in these cases with the complexity  of the FFT. Following
there is an efficient implementation of this transform for generic time and  frequency
grids, and a short demo for validation.  Just place all the code in a file named
Demo.m and execute it.

function Demo

while 1
  
  %Frequency grid with arbitrary start value, spacing, and number of elements. ...
  %Mathematically, the grid is
  %
  %   fo+Df*m, m=0, 1,...,M-1.
  %
  
  fo = 10*rand(1)-0.5; %Start value.
  Df = rand(1)+0.01;   %Frequency spacing. 
  M = round(rand(1)*1000); %Number of frequencies. 

  %Time grid with arbitrary start value, spacing, and number of elements. Mathematically, ...
  %the grid is
  %
  %   to+Dt*n, n=0, 1,...,N-1.
  %

  to = 10*rand(1)-0.5; %Start value; (instant of first sample).
  Dt = rand(1)+0.01;   %Time spacing. 
  N = round(rand(1)*1000); %Number of samples. 
  
  %Random vector of samples.
  x = randn(1,N)+1i*randn(1,N);
  
  %x(n) is the sample at instant to+Dt*(n-1), n=1,2,...,N. 

  %We proceed to compute the DFT of x at frequencies fo+m*Df, m=0, 1,..., M-1.

  %Compute DFT using the direct and chirp-z methods.
  
  tic, XD = DirectGridDFT(x,to,Dt,fo,Df,M); tmD = toc;
  tic, XF = chirpzDFT(x,to,Dt,fo,Df,M); tmF = toc;
  
  disp(['Timing direct method: ', num2str(tmD), ' sec.'])
  disp(['Timing chirp z: ',num2str(tmF),' sec.'])
  disp(['Relative difference: ', num2str(max(abs(XD-XF))/max(abs(XD)))])
  disp(' ')
  input('(Press a key for another trial or Ctrl-C) ')
  disp(' ')
end

function X = chirpzDFT(x,to,Dt,fo,Df,M)

%Author: J. Selva.
%Date: November 2011.
%
%This function computes the DFT for two regular time and frequency grids with arbitrary
%starting points and spacings, using the Chirp-z Transform. Specifically, it computes
%
%               N
%       X(k) = sum  x(n)*exp(-j*2*pi*(fo+Df*(k-1))*(to+Dt*(n-1))), 1 <= k <= M.
%              n=1
%
%Input parameters:
%
%  x: Signal samples. 
% to: Instant of first sample in x.
% Dt: Time spacing between consecutive samples of x.
% fo: First frequency in which the Fourier sum is computed.
% Df: Spacing of the desired frequency grid.
%  M: Desired number of output samples. 
%
%The algorithm is explained in Sec. 9.6.2, page 656 of
%
% A.V. Oppenheim, R. W. Schafer, J. R. Buck, "Discrete-time signal processing," second
% edition, 1998.
%

x = x(:).';

N = numel(x);
P = 2^ceil(log2(M+N-1));

%The next three lines do not depend on the vector x, and so they can be pre-computed if 
%the time and frequency grids do not change, (i.e, for computing the transform of 
%additional vectors). Thus, this algorithm just involves two FFTs for fixed grids and 
%three if the grids change. 

a = exp(-1i*2*pi*((1:N)*Dt*(fo-Df)+Dt*Df*(1:N).^2/2));
b = exp(-1i*2*pi*((to-Dt)*(fo+(0:M-1)*Df)+Dt*Df*(1:M).^2/2));
phi = fft(exp(1i*2*pi*Dt*Df*(1-N:M-1).^2/2),P); %FFT of chirp pulse. 

%Weigh x using a and perform FFT convolution with phi. 
X = ifft( fft(x.*a,P) .* phi );

%Truncate the convolution tails and weigh using b. 
X = X(N:M+N-1) .* b;

function X = DirectGridDFT(x,to,Dt,fo,Df,M)

%Direct (and slow) version of the Fourier sum with arbitrary and regular time and
%frequency grids. Used for validation of chirpzDFT.

x = x(:).';

N = length(x);

X = zeros(1,M);

for k = 1:M
  for n = 1:N
    X(k) = X(k) + x(n)*exp(-1i*2*pi*(to+Dt*(n-1))*(fo+Df*(k-1)));
  end
end