**Language:** Matlab

**Processor:** Not Relevant

**Submitted by Jesus Selva on Nov 17 2011**

Licensed under a Creative Commons Attribution 3.0 Unported License

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

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

Jesus Selva's background is in Communications Engineering and Mathematics. He is currently a research fellow at University of Alicante (Spain).

Comments

No comments yet for this code