| << 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 [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 polynomial whose values on the unit circle are the DFT is a slightly modified z-transform of x(n) given by
which for clarity in this development uses a positive exponent .This is illustrated for a length-4 sequence as a third-order polynomial by
The DFT is found by evaluating [link] at , which can be written as
where
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
This nested sequence of operations can be written as a linear difference equation in the form of
with initial condition , and the desired result being the solution at . The value of the polynomial is given by
[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 being the filter output sampled at . Applying this to the DFT gives the Goertzel algorithm [link] , [link] which is
with and
where
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 trigonometric evaluations rather than . It also results in an increase in accumulated quantization error.
Notification Switch
Would you like to follow the 'Fast fourier transforms' conversation and receive update notifications?
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|