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.
%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
; %Start value.
Df = rand(1)
; %Frequency spacing.
M = round(rand(1)
; %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
; %Start value; (instant of first sample).
Dt = rand(1)
; %Time spacing.
N = round(rand(1)
; %Number of samples.
%Random vector of samples.
x = randn(1
%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.
, XD = DirectGridDFT(
; tmD = toc
, XF = chirpzDFT(
; tmF = toc
disp(['Timing direct method: '
, ' sec.'])
disp(['Timing chirp z: '
disp(['Relative difference: '
input('(Press a key for another trial or Ctrl-C) ')
X = chirpzDFT(
%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
% X(k) = sum x(n)*exp(-j*2*pi*(fo+Df*(k-1))*(to+Dt*(n-1))), 1 <= k <= M.
% 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(
P = 2
%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
b = exp(-1i
phi = fft(exp(1i
; %FFT of chirp pulse.
%Weigh x using a and perform FFT convolution with phi.
X = ifft( fft(
.* phi )
%Truncate the convolution tails and weigh using b.
X = X(
X = DirectGridDFT(
%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 = zeros(1
k = 1
n = 1
Rate this code snippet:
posted by Jesus Selva
Jesus Selva's background is in Communications Engineering and Mathematics. He is currently a research fellow at University of Alicante (Spain).