Sign in

Not a member? | Forgot your Password?

Search code

Search tips

Free PDF Downloads

A Quadrature Signals Tutorial: Complex, But Not Complicated

Understanding the 'Phasing Method' of Single Sideband Demodulation

Complex Digital Signal Processing in Telecommunications

Introduction to Sound Processing

C++ Tutorial

Introduction of C Programming for DSP Applications

Fixed-Point Arithmetic: An Introduction

Cascaded Integrator-Comb (CIC) Filter Introduction

FIR Filter Design Software

See Also

Embedded SystemsFPGA

DSP Code Sharing > FFT for arbitrary time and frequency grids (Chirp-z Transform)


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

Language: Matlab

Processor: Not Relevant

Submitted by Jesus Selva on Nov 17 2011

Licensed under a Creative Commons Attribution 3.0 Unported License

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

 
 
Rate this code snippet:
4.5
Rating: 4.5 | Votes: 2
 
   
 
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).


Comments


No comments yet for this code


Add a Comment
You need to login before you can post a comment (best way to prevent spam). ( Not a member? )