Two-dimensional chirp Z-transform

Signal Processing 4 (1982) 17-24 North-Holland Publishing Company

17

T W O - D I M E N S I O N A L CHIRP Z - T R A N S F O R M D. R A G H U

R A M R E D D Y and V.V. R A O

Department of Electrical Engineering, Indian Institute of Technology, Madras-600 036, India

Received 10 October 1980 Revised 16 January 1981

Abstract. It has been established that the CCD implementation of a one-dimensional chirp z-transform (CZT) is fairly simple and inexpensive. The two-dimensional chirp .--transform (2-D CZT) is developed and its CCD implementation is suggested. The computational efficiency of CZT algorithm over the Fast Fourier Transform (FFT) algorithm in evaluating some Discrete Fourier Transforms (DFT) with a specified resolution is established. Zusammenfassung. Die Realisierung der eindimensionalen Chirp-_--Transformation (CZT) mit Hilfe yon Ladungstransferschaltungen (CCD-Technologie) ist verh~iltnism~ii3ig einfach und billig. In vorliegendem Beitrag wird die zweidimensionale CZT hergeleitet und ihre Realisierung mit Hilfe der CCD-Technologie vorgeschlagen. Die Effizienz der CZT beziiglich des Realisierungsaufwandes im Vergleich zur schnellen Fouriertransformation wird nachgewiesen fiir einige F~ille der diskreten Fouriertransformation mit vorgegebener spektraler Aufl6sung. R~sum~. I1 a ~t~ ~tabli que la r~alisation par CCD de la transformation en z utilisant une modulation de fr6quence (chirp z-transform - CZT) est relativement simple et 6conomique. La transformation CZT bidimensionnelle est d6velopp6e et sa r~alisation par CCD est sugg6r~e. L'efficacit~ de calcul de la transformation CZT par rapport h la transformation de Fourier rapide/FFT) est 6tablie dans I'~valuation de certaines transform~es de Fourier discr&es avec une r~solution sp~cifi~e. Keywords. 2-D CZT, CCD implementation, FFT, 2-D FFT, computational efficiency.

1. Introduction There are many signals that are inherently two dimensional viz., photographic data, air reconnaisance photos, medical records, seismic records gravity and magnetic data etc. Although twodimensional signals may be processed by onedimensional systems, it is generally preferable to consider using two-dimensional systems. There are many basic ideas of one-dimensional signal processing that have been extended to twodimensional systems [1], but there are some concepts of one-dimensional systems that have not yet been extended to two-dimensional systems. One such concept is the chirp z-transform which has not been extended to 2-D systems. 0165-1684/82/0000-0000/$02.75

The method proposed by Blustein [2] views the D F T as a convolution with a nonrecursive chirp filter, a potentially cost efficient alternative to the Fast Fourier Transform (FFT) and the more recent Winograd transforms. The formulation was generalized for one-dimensional systems by Rabiner et al. [3] and was given the name CZT. One efficient method for implementing this algorithm is to perform the required convolution through digital FFTs [3]. Direct convolution, however is now also computationally tractable due to the recent technological advances in charge coupled devices (CCDs) [4, 5]. This is because CCDs allow real time operation of nonrecursive filters at sampling rates up to about 5 MHz. Moreover the onedimensional C Z T has proven itself to be less costly

O 1982 North-Holland

D.R.R. Reddy, V. V. Rao / Two-dimensionalCZT

18

with CCD implementation and simpler to implement than its digital counterparts [4, 5]. As a natural extension to two-dimensional signals, the 2-D C Z T is developed and its CCD implementation is discussed. Further, the computational efficiency of the two-dimensional C Z T algorithm over the 2-D Fast Fourier Transform (2-D FFT) in evaluating discrete Fourier transform (DFT) values on a part of the two-dimensional hyperplane is established.

Let X ( k , l ) be the 2-D C Z T on the surtace described by eq. (2). Then M-1

N-1

X(k, l) = ~

~ x(m, n ) A - " w ' k B - " U "t.

m=O

n=0

(4) By using the relations

2kin = m2 + k E - ( k - m ) 2, 2In = n2 + 1 2 - ( l - n ) 2 in eq. (4), we obtain M-I

N-1

X(k, l) = •

2. Two-dimensional chirp z-transform

~ x(m, n ) A - " B - " W ~k:. . . .

m=O

Let x(m, n) be a 2-D sequence where m and n are integer variables. The 2-D z-transform for the sequence x(m, n) is given by [3]

¿k-m)2)~2

n=O

x U "'-+"2-~t-"m/2

(5a)

On rearranging and simplication, we have M-1

X(k, 1)= W k2/2 Z U12/2 m~--oo

m=0

tt = - c o N-I

For a finite sequence which is M x N points in extent, the z-transform is M-1

X(Z1, Z2) = ~ m=O

Z 2 = B U -I,

~ x(m,n)Z-1"Z2".

× A-'~ W m2/2W

n~O

(2) /=0,1 ..... L-1

A = Ao exp(j2rrao),

W = Wo exp(j2rr&o), U = Uo exp(j2w0o).

(5b)

M - 1

k=O, 1 . . . . . K - l ,

B = Bo exp(j2rrflo),

-(k -m)2/2

(1)

where K and L are arbitrary integers. A, B, W and U are arbitrary complex numbers defined by

Signal Processing

n=O

N-I

The chirp z-transform is an efficient algorithm for evaluating the z-transform of a finite duration sequence on certain general surfaces in the 2-D plane. Since a bi-circte is one of the allowable surfaces, the C Z T algorithm can be used to evaluate the 2-D D F T of the sequence efficiently. The C Z T eliminates many restrictions of the FFT. The C Z T algorithm can efficiently evaluate the D F T of a sequence, even when M and N are prime numbers. To evaluate the 2-D CZT, let

Z I = A W -k,

X ~ x(m,n)B-nU"2/2U -~t-'~2/2

(3)

= W k~/z Y~ y ( m , I ) W

-~-'"~/2

(5c)

m=O

where y(m, I ) = g(m, I ) A - m W m~/2, N-I

g(m, l)= U t2/2 ~. f(m, n ) U -(t-"~2/z,

(5d)

n=0

f(m, n) = x(m, n ) B - " U "2/2. It is observed that 2-D C Z T is a separable process. That is the one-dimensional algorithm can be used to operate on the variables m and n separately. 2-D C Z T requires that the 1-D C Z T be performed on each row of the input array and then on each column of the resulting array. The schematic diagram which summarizes the operations of 2-D C Z T is given in Fig. l(a). It is evident that the 2-D C Z T needs M one-dimensional CZTs to be operated on M rows of the input sequence x(m, n) and L one-dimensional CZTs to be performed on L columns of the resulting array. This is indicated in Fig. 1.

19

D.R.R. Reddy, V. V. Rao / Two-dimensional CZT r

1

i t l

I :,(,,,.n),

r

I

I

o o.o,o,,o ,,(,-,,.n~\^f

~

c = B - n U n 12

.,',,n~'U-

g(m.l)

,..:o I i=q 1,1 k£J

d = U - n 212

'~

"

I ! t" I

g,m,, .t.-!l stor=ge ]t=0I..L-1 ] d,,~,o, I I

l g(m,l)

__

kL.J

~=A-

e = U - I 12

,

W

coo.o,o.ou

I !

,,

lyO.,,.Oew ..... 21 ~.- J ,

12

q=W-

12

s=W

x~k.,~

/2

Fig. l(a). Operations in 2-D chirp z-transform.

vQryin 9 with

n

-

x (re,n) m=0t M-I

for the same N a n d L

cj(m.I ) rn=0 .M-1

n=0,1

-

I =0, .. L - I

.N- 1

v Q r y i n g wilh m fc¢ the same M Qnd K,

X(k.I} k= O,L.,K-I I = O,I.,..L-I

Fig, lib). Block diagram of 2-D chirp z-transform. It is to be noticed that one-dimensional CZTs are repeated M times with input sequence of length N, evaluating the CZTs at L points. Hence the constants c, d, e and the DFT of U -"s/'- (Fig. 1) can be computed once and stored. Similarly, the same applies to other set of L one-dimensional CZTs. The details of one-dimensional CZT algorithm on the assumption that an already existing FFT program is available to compute DFTs and IDFTs, is given in references [1, 3]. The additional freedoms over 2-D FFT and suitable application offered by 2-D CZT are summarized as follows: (a) The number of time samples does not have to equal the number of samples of z-transform. (b) M, N, L and K need not be composite numbers. (c) The angular spacings are arbitrary. (d) The evaluation of z-transform can be over any plane inside, outside and on bi-unit disc for any arbitrary resolution. In addition, the starting point is arbitrary, but this is also the case with FFT if the samples are multiplied by exponentials before transforming. (e) By evaluating the z-transform along a selected two-dimensional hyperplane, the transfer function of a linear system can be adjusted to sharpen or flatten the frequency response. The choice of the plane is dependent on which of the systems analytic regions should be emphasized.

(f) It can evaluate high resolution narrow frequency band at relatively low cost providing high resolution for a limited range of frequencies as CZT algorithm allows arbitrary selection of starting frequency and is independent of the number of time samples. (g) Interpolation between time samples of band limited function can be obtained with greater ease. Interpolation at arbitrary sampling intervals can be performed effectively by using the 2-D CZT wherein it is used efficiently to evaluate the 2-D

FF'r [3].

3. Functional diagram and CCD implementation of 2-D CZT as spectral analyzer The implementation of the algorithm is derived from Fig. l(a) as natural extension of onedimensional case [6]. In this derivation, the real and imaginary quantities are represented by r and i respectively. Let C = B - n U n2/2 = (Cr+jCi).

Hence f r( rn, n) = c~x~( m, n ) - cixi( rn, n ),

(6)

fi(m, n ) = G x i ( m , n ) + cix~(m, n ).

(7)

Vol. 4, No. 1, J a n u a r y 1982

D.R.R. Reddy, V. V. Rao / Two-dimensional C Z T

20

e, = cos wle/N,

Let

d = U -'2/2 = (dr + jdi). H e n c e (* r e p r e s e n t s c o n v o l u t i o n )

ar(m, 1) = f r ( m , n) * dr-fi(m, n) * d i,

(8)

ai(m, 1) = r a m , n) * di + fi(rn, n) * dr.

(9)

ei = - s i n

wl"/N,

pr = c o s ~ m ' / m ,

pi = - s i n rrrn21M,

qr = COS 7"rm2/m,

qi = sin ~rrn21M,

s~ = c o s

9

rrk-/M,

si = - s i n 7rkZlM.

S u b s t i t u t i n g the a b o v e v a l u e s in eqs. (6) to (17), we h a v e

Let

e = U t2/2 = (er + j e i ) .

fr(m, n)= xAm, n) cos wn2 / N

Hence

+ xi(m, n ) sin rrn 2/N,

gr(m, 1) = a~(m, l)er- ai(m, l)ei,

(10)

gi(m, l) = ar(m, l)ei + ai(m, l)er.

(11)

Let

fi(m, n) = - x A m , n) sin ,"rn2/ N + xi(m, n) cos wnZ/N,

a,.(m, l) = fr(m, n) * cos wn 2IN p = A - ' W ''2/z = (pr + jPi).

-fi(m, n) * sin "rrnZ/N,

Hence

ai(m, l) = f i ( m , n) * sin ,"rnZ/N

y~(m, l) = gr(m, l)p~- gi(m, l)pi,

(12)

yi(rn, l ) = g~(m, l)pi + gi(m, l)p,

(13)

Let

+ fi(rn, n) * cos wnZ/N,

g,.(rn, l) = a~(rn, l) cos vlZ / N + ai(m, l) sin 7rlZ/N,

q -~- W -m2/2

=

(qr+jqi)-

gi(m, l) = -a,.(m, l) sin wlZ/N

Hence

+ ai(m, l) cos 7rlZ/N,

br(k, l) = yr(m, l) * q ~ - yi(m, l)

* qi,

(14)

bi(k, l) = y~(m, 1) * qi + yi(m, l) * qr.

(15)

Let

y~(m, l) = gr(m, l) cos w m Z / M + gi(m, l) sin ~ m Z / m , yi(m, l) = -gr(m, l) sin v m Z / M

S = W k2/2 = (St + j s i ) .

+ gi(m, l) cos wm2/l"v[,

br( k, l) = y,.(m, l) * cos 7rrnZ/ M

Hence

X~(k, l)= br(k, l)s~-bi(k, l)si,

(16)

Xi( k, l)= bAk, l ) s i + bi( k, l)s~.

(17)

A s a s p e c i a l case, C Z T can be m a d e use of as a s p e c t r a l a n a l y z e r [6]. In such an i m p l e m e n t a t i o n , w e have A - - l , B = I , Wo=l and Uo=l. Also let 0o = - 1 / M

ci = - s i n rrn"/N,

dr=cos rrn2/N,

di=sin wn2/N,

Signal Processing

- yi(m, 1) * sin wrn21M, bi(k, I) = yr(m, l) * sin ,"rm21M

+ yi(m, l) * cos "rrmZ/M,

xr(k, l) = br(k, l) cos v k 2 / M + hi(k, l) sin ,rk2/M,

and/30 = -1/N. Then

c~ = cos "nn2/N,

(18)

(18)

xi(k, l) = -br(k, l) sin v k 2 / M

+ bi(k, l) cos rrk 2/M.

D.R.R. Reddy, V. V. Rao / Two-dimensional CZT

The functional diagram is shown in Fig. 2 including all the above values in real quantities. The functional diagram of Block of Fig. 2 is given in Fig. 3(a). The application of CCDs to the CZT algorithm requires the use of tapped and untapped delays [6] to perform the convolution. The CCD implementation of 1-D CZT (Block 1 of Fig. 2) is presented in Fig. 3(b). The operation and the CCD implementation of Block 2 of Fig. 2 is identical to that of Block 1. The execution of the set of L one-dimensional CZTs starts only after completion of the execution of the set of M onedimensional CZTs. 4.

Computational

~

m

~

= -~fl n = N:t

n

• 2w x e x p ( - j ff-~km),

--I

I

= 0,

xi(m.n )

m~..~A-1

b,,t-l~N-,

Y(k, l) = -1=oE ,,=oY y(m, n) exp

( .2~r 1)

-j~--~n

x exp(-j b~mk ) .

(21)

X ( k , l) is obtained by selecting the appropriate M N values of Y(k, l). The number of operations

needed for the two-dimensional FFT [1, 7] are: 0 . 5 a b M N log2 a b M N complex multiplications, and a b M N loga a b M N complex additions/subtractions.

Co$1rk 2IM

xr(k,= )

gr(m,I ) I =O,I..L- 1 STORAGE

I gi( m,I )

[n:O,L..N.-,[ ~ _ _ _ _

Jm:O,L.~M_, 1 =O,I..L -1

-SinTFI21N

(20)

elsewhere.

CoslTm2/M

[ xi(m_.n ) J

- Sin ITn21 N

MI<~m<-M=,

The DFT of y(m, n) with resolutions 2-~/aN and 2"rr/ bM is

i m =O,EM-1

BLOCK 1

(19)

Nt <~n <~N2,

CosVn2/N Cos I-I 12/N Set of M,1 -D | CZ Ts Xr(m.n)_ IXr(m-n), -

M=,

where a and b are integers, M1 and N~ are starting points and M2 and N2 are ending points such that N=Nz-Nt+I and M = M 2 - M ~ + I . X ( k , I ) can be evaluated either using the CZT or the FFT. In the comparison of these two methods, a radix-2 FFT is employed. Hence, a, b, N and M are all power of 2 so that X ( k , I) can be computed using a radix-2 FFT. The two-dimensional FFT is employed to evaluate the DFT of the following sequence. y(m,n)=x(m,n),

x(m,n)exp -j

m=0,1...M-1 'n=ql.'~N-l] .

M1, Ml + 1 . . . . .

I = N1, NI + I . . . . . N:

e f f i c i e n c y of 2 - D C Z T

It is well known [1, 7] that CZT can be efficiently used in evaluating DFT if number of samples is prime or not a power of some radix for which fast Fourier Transform (FFT) algorithms are existing. So far no investigation has been reported regarding the resolution for which CZT is more efficient in evaluating the DFT of N-samples when N is a power of some radix for which efficient FFT is available. In this section, a two-dimensional analysis is presented to the above problem. Let x(m, n) be the M N sequence. X(k,l)=

k =

21

DEVICE

k'= O...K-I I = O,..L-I

I=O...L-| [rn: 0,I...1

[M-, IBLOCK 2 J gitm,l)Jgr~m,t 11

xi( k,1 )

L-rL ,,I

Ik --'o..K-;

- ira=O, ..l

I I = O,..L-I

-Sin 1Tin2/M - Sin I"I"k2/M

Fig. 2. Functionaldiagram of 2-D chirp z-transform as spectral analyzer. Vol. 4, No. 1, January 1982

22

D.R.R. Reddy, V. V. Rao / Two-dimensional CZT

Cos~n2/N 1 ~

-S,n~'n21N

COST12I TN ~g6m~t

-~

It may be noted that computation involved is independent of the number of frequency coefficients required for the specified resolutions 2 r r / a N and 2~r/bM. The two-dimensional C Z T (Section 2) can be evaluated using eqs. (4) and (5). X ( k , l) of eq. (19) can be easily obtained on substituting the following in eq. (5). A = exp(j2rr/bM), B = exp(j2.rr/aN),

(23)

W = exp(j2rr/bM), Ct~$1Tn21N

Cosg I 21N

U = exp(j2rr/aN)

Fig. 3(a). Functional diagram of Block 1 of Fig. 2.

COSITn2/N

K=M

CosITt21N

Fig. 3(b). CCD implementation of Block I of Fig. 2.

A complex multiplication is equal to four real multiplications, one real addition and one real subtraction. A complex addition/subtraction is equal to two real additions/subtractions. As one addition is approximately equal to one subtraction, we have M U L = 2 a b M N log2 a b M N ,

and

L=N.

The 2-D C Z T is evaluated using convolutions with FFT for P and Q samples [1], where P and Q are chosen as the least values greater than ( N + L - 1) and ( M + K 1) respectively and are powers of 2. In this case P = 2 N and Q = 2M. In the case of one-dimensional C Z T the arithmetical operations needed are [1, 7]: (N + L + P + P log2 P) complex multiplications; 2Plog2 P complex additions/subtractions and 2 P real multiplications (or 2P real divisions). 2 P real multiplications (or 2P real divisions) are required to obtain the IDFT. The arithmetical operations needed for generation of constants c, d, e (Fig. 1) and for the evaluation of the D F T of U -"2/2 (Fig. 1), are not counted because these are computed as needed and prestored in the case of repeated CZTs with the same N and L as applicable in the evaluation of the 2-D CZT. Total number equal arithmetical operations needed when L = N and P = 2N, are: M U L = 2 0 N + 8 N log2 2N, A D D = 8 N + 12N log2 2N.

A D D = 3 a b M N log2 a b M N . If a and /3 are respectively the execution times for one real multiplication and one real addition, the time required by the 2-D FFT method is TOT1 = (2a + 3 ~ ) a b M N log2 a b M N . Signal Processing

(22)

The one-dimensional time

C Z T requires execution

T ( N ) = (20N + 8 N log,. 2N)o~

+ ( 8 N + 1 2 N log2 2N)/3.

D.R.R. Reddy, V. V. Rao / Two-dimensional C Z T

2-D C Z T with K = M execution time

and Q = 2 M

requires

T O T 2 = MT(N) + NT(M) = M N ( 4 0 + 8 log2 4MN)a + M N ( 1 6 + 12 log2 4MN)~.

(24)

2-D C Z T is superior computationally in evaluating

S(k, l) if (40 + 8 log2 4MN)a + (16 + 12 log2 4MN)fi < ab(log2 abMN)(2a + 3/3), i.e. log2 4MN(8a + 12/3) + (40a + 16fl) < ab(log2 abMN)(2a + 3fl).

(25)

The minimum values of a and b can be found such that eq. (25) is satisfied for various values of M and N'. For all values of a and b greater or equal to the minimum values, 2-D C Z T is superior to 2-D FFT in evaluating X(k, l) of eq. (19). For various values of M and N, the minimum values of a and b have been calculated and it is observed that the values a and b are different for different values of M and N for all values of a and ft. In the special case of M = N and the resolutions in both the directions being the same 2~r/aN (=2~r/bM), it is found that 2-D C Z T is having computational advantage for any resolution less than or equal to (2-rr/4N)= (2~r/4M) (i.e. minimum value of a or b is 4) for any value of N (or M ) which is a power of 2 for all values of a and/3. It may be observed that the one-dimensional case is a special case of 2-D one. x(/)=

~

x(n)exp

,

(26)

23

the one-dimensional FFT in evaluating x(1) on the arc of length 2,-r/4 where N is a power of 2. This is true if the one-dimensional C Z T is performed repeatedly for the same N and L. If the C Z T is not repeated the one-dimensional C Z T is more efficient than one-dimensional FFT in evaluating X(l) on the arc of length 2-rr/8 since the arithmetical operations necessary for generation of constants c, d and e (Fig. 1) and for evaluation of the D F T of U -"'-/2 (Fig. 1) has to be considered in this case. It may be recalled that one or two-dimensional C Z T can be advantageously used in evaluating D F T in all cases when a and b or M and N are not power of any radix for which efficient F F T is available.

5. Conclusions The concept of two-dimensional chirp z-transform is introduced and it is established that the process is separable and can be implemented with ML one-dimensional CZTs. CCD implementation of the 2-D C Z T is suggested. Conditions for the computational superiority of the 2-D C Z T over the 2-D FFT in evaluating D F T of a sequence with a specified resolution in the frequency domain are presented. Further it is observed that the 2-D C Z T is efficient in evaluating the 2-D D F T of M N samples with M = N (case of equal resolutions in both directions) when the required resolution is less than or equal to 2~r/4N. The effects of finite word length in the digital implementation of the C Z T algorithm are under investigation.

References

n=N 1

l = N1, N1 + 1 . . . . , N2.

x(1) is D F T of N samples evaluated on arc of unit circle of length 2rr/a with resolution 2rrlaN. Hence it can be inferred from the 2-D case that the one-dimensional C Z T is more efficient than

[1] L.R. Rabiner and B. Gold, Theory and Applications of Digital Signal Processing, Prentice-Hall, Englewood Cliffs, NJ, 1975, Ch. 6, pp. 356--437. [2] L.I. Blustein, Several Fourier transform algorithms, Sylvania Electron Syst. Tech. Note 49, Feb. 1968. [3] L.R. Rabiner et al., The chirp z-transform algorithm and its applications, Bell Syst. Tech., Vol. 48, No. 5, May-June 1969, pp. 1249-1292. Vol. 4, No. 1, January 1982

24

D.R.R. Reddy, V. V. Rao / Two-dimensional C Z T

[4] R.W. Erodersen et al., A 500 stage CCD transversal filter for spectral analysis, I E E E J. Solid State Circuits, Vol. 8C-11, Feb. 1976, pp. 75-84. [5] D.D. Buss et al., Communication applications of CCD transversal filters, presented at 1975, Nat. Telecomm. Con[., New Orleans, LA, 1975.

Signal Processing

[6] G.J. Mayer, The chirp z-transform - A CCD implementation, R C A Reeiew, Vol. 36, Dec. 1975, pp. 759-773. [7] A.V. Oppenheim et al., Digital Processing, Prentice-Hall, Englewood Cliffs, NJ, 1975, Ch. 6, pp. 284-336.