<< Chapter < Page Chapter >> Page >

Two Matlab programs to calculate an arbitrary length DFT using the chirp z-transform is shown in [link] .

function y = chirpc(x); % function y = chirpc(x)% computes an arbitrary-length DFT with the % chirp z-transform algorithm. csb. 6/12/91% N = length(x); n = 0:N-1; %Sequence lengthW = exp(-j*pi*n.*n/N); %Chirp signal xw = x.*W; %Modulate with chirpWW = [conj(W(N:-1:2)),conj(W)]; %Construct filtery = conv(WW,xw); %Convolve w filter y = y(N:2*N-1).*W; %Demodulate w chirpfunction y = chirp(x); % function y = chirp(x)% computes an arbitrary-length Discrete Fourier Transform (DFT) % with the chirp z transform algorithm. The linear convolution% then required is done with FFTs. % 1988: L. Arevalo; 11.06.91 K. Schwarz, LNT Erlangen; 6/12/91 csb.% N = length(x); %Sequence lengthL = 2^ceil(log((2*N-1))/log(2)); %FFT length n = 0:N-1;W = exp(-j*pi*n.*n/N); %Chirp signal FW = fft([conj(W), zeros(1,L-2*N+1), conj(W(N:-1:2))],L); y = ifft(FW.*fft(x.'.*W,L)); %Convolve using FFTy = y(1:N).*W; %Demodulate

Goertzel's algorithm (or a betterDft algorithm)

Goertzel's algorithm [link] , [link] , [link] is another methods that calculates the DFT by converting it into a digital filtering problem. Themethod looks at the calculation of the DFT as the evaluation of a polynomial on the unit circle in the complex plane. This evaluation isdone by Horner's method which is implemented recursively by an IIR filter.

The first-order goertzel algorithm

The polynomial whose values on the unit circle are the DFT is a slightly modified z-transform of x(n) given by

X ( z ) = n = 0 N - 1 x ( n ) z - n

which for clarity in this development uses a positive exponent .This is illustrated for a length-4 sequence as a third-order polynomial by

X ( z ) = x ( 3 ) z 3 + x ( 2 ) z 2 + x ( 1 ) z + x ( 0 )

The DFT is found by evaluating [link] at z = W k , which can be written as

C ( k ) = X ( z ) | z = W k = D F T { x ( n ) }

where

W = e - j 2 π / N

The most efficient way of evaluating a general polynomial without any pre-processing is by “Horner's rule" [link] which is a nested evaluation. This is illustrated for the polynomial in [link] by

X ( z ) = [ x ( 3 ) z + x ( 2 ) ] z + x ( 1 ) z + x ( 0 )

This nested sequence of operations can be written as a linear difference equation in the form of

y ( m ) = z y ( m - 1 ) + x ( N - m )

with initial condition y ( 0 ) = 0 , and the desired result being the solution at m = N . The value of the polynomial is given by

X ( z ) = y ( N ) .

[link] can be viewed as a first-order IIR filter with the input being the data sequence in reverse order and the value of thepolynomial at z being the filter output sampled at m = N . Applying this to the DFT gives the Goertzel algorithm [link] , [link] which is

y ( m ) = W k y ( m - 1 ) + x ( N - m )

with y ( 0 ) = 0 and

C ( k ) = y ( N )

where

C ( k ) = n = 0 N - 1 x ( n ) W n k .

The flowgraph of the algorithm can be found in [link] , [link] and a simple FORTRAN program is given in the appendix.

When comparing this program with the direct calculation of [link] , it is seen that the number of floating-point multiplications and additionsare the same. In fact, the structures of the two algorithms look similar, but close examination shows that the way the sines and cosines enter thecalculations is different. In [link] , new sine and cosine values are calculated for each frequency and for each data value, while for the Goertzel algorithm in [link] , they are calculated only for each frequency in the outer loop. Because of the recursive or feedback natureof the algorithm, the sine and cosine values are “updated" each loop rather than recalculated. This results in 2 N trigonometric evaluations rather than 2 N 2 . It also results in an increase in accumulated quantization error.

Get Jobilize Job Search Mobile App in your pocket Now!

Get it on Google Play Download on the App Store Now




Source:  OpenStax, Fast fourier transforms. OpenStax CNX. Nov 18, 2012 Download for free at http://cnx.org/content/col10550/1.22
Google Play and the Google Play logo are trademarks of Google Inc.

Notification Switch

Would you like to follow the 'Fast fourier transforms' conversation and receive update notifications?

Ask