- Email: [email protected]
Discrete Gabor transforms with complexity O (NlogN)
Signal Processing 77 (1999) 159}170
Discrete Gabor transforms with complexity O(NlogN)夽 Sigang Qiu *, Feng Zhou, Phyllis E. Crandall Ambient Technologies Inc., NCDC, 110 Horizon Drive, Raleigh, NC 27615, USA PowerPC Microprocessor Validation, Department RYMA, Bldg. 667, Research Triangle Park, NC 27709, USA Los Alamos National Laboratory, P.O. Box 1663, ACL, MS B287, Los Alamos, NM 87545, USA Received 27 August 1997; received in revised form 11 February 1999
Abstract We develop fast computations of discrete Gabor transforms. We present algorithms for computations of both Gabor analysis and Gabor synthesis via the fast Fourier transform. The algorithms are performed with complexity O(N log N). In addition, we propose an algorithm for computing the dual Gabor wavelets with complexity less than O(N log N). 1999 Elsevier Science B.V. All rights reserved. Zusammenfassung Wir entwickeln schnelle Berechnungsmethoden fuK r diskrete Gabor-Transformationen. Wir praK sentieren Algorithmen fuK r die Berechnung der Gabor-Analyse sowie der Gabor-Synthese mittels der Schnellen Fouriertransformation (FFT). Diese Algorithmen weisen eine KomplexitaK t von O(N log N) auf. Weiters schlagen wir einen Algorithmus zur Berechnung des dualen Gabor-Wavelets mit einer KomplexitaK t von weniger als O(N log N) vor. 1999 Elsevier Science B.V. All rights reserved. Re2 sume2 Nous proposons dans cet article des meH thodes de calcul rapide de transformations de Gabor discre% tes. Nous preH sentons des algorithmes de calcul a` la fois pour l'analyse et la synthe`se de Gabor via la transformation de Fourier rapide (en anglais FFT). Ces algorithmes ont une complexite` O(N log N). De plus, nous proposons un algorithme pour le calcul d'ondelettes de Gabor duales ayant une complexiteH infeH rieure a` O(N log N). 1999 Elsevier Science B.V. All rights reserved. Keywords: Discrete Gabor transforms (DGT); Gabor analysis; Gabor synthesis; Fast Fourier transform
1. Introduction The interest in Gabor transforms arises from the fact that Gabor basis functions has good time}frequency localization [6,7] with appropriate Gabor window functions. Given a Gabor window u and a pair of lattice parameters (a,b), there are two basic and fundamental questions. 夽 The authors were with the University of Connecticut when this paper was submitted. The referee's insightful comments on the paper are acknowledged. * Corresponding author. Tel.: #1-919-870-5060; fax: #1-919-870-9898. E-mail address: [email protected] (S. Qiu)
0165-1684/99/$ - see front matter 1999 Elsevier Science B.V. All rights reserved. PII: S 0 1 6 5 - 1 6 8 4 ( 9 9 ) 0 0 0 3 0 - 4
160
S. Qiu et al. / Signal Processing 77 (1999) 159}170
E Gabor decomposition (analysis): For a signal x, how can we e$ciently "nd the Gabor coe$cients +c , L K such that x is represented (or best approximated) by the following form: ? \ @I \ x" c u , (1) L K L K L K where +u "M ¹ u, are Gabor basis (or Gabor elementary) functions. ¹ and M are the time L K K@ L? L K L? K@ translation operator and modulation operator. We call u a Gabor window function. a "N/a and bI "N/b are referred to as the dual lattice parameters. E Gabor reconstruction (synthesis): Given a set of coe$cients +c ,, how can we e$ciently determine L K a signal x with the form of Eq. (1). As Orr [11] observed, there are several approaches to the above questions. From Orr's summary [11, Table 1], we may say that matrix methods are comparatively expansive; however, there are many papers that describe matrix approaches to the topic. Variant matrix techniques have proven to be e!ective and clear. We refer the interested reader to some of the relevant papers [1}5,9,12}20] and references cited in each. We consider the critical Gabor transform in this paper. The aim is to show that both the Gabor analysis and synthesis can be performed with the same complexity as FFT. The paper is organized in the following way. Relevant notation is listed in Section 2. In Section 3 we investigate the structures of Gabor basic matrix. Then we are able to propose relevant algorithms in Section 4. We show that the relevant computations are performed with a couple of FFT. A necessary and su$cient condition for (u, a, b) to generate a Gabor frame is formulated. Numerical simulations are given in Section 5. We make our conclusions in Section 6.
2. Notation Signals are considered as row vectors in ",, and we identify signals with N-periodic complex sequences. That means, for a signal x"(x( j)),\3",, we treat x as an N-periodic sequence such that x( j#kN)"x( j) H for j"0,1,2,N!1, k39. We use bold upper case letters to denote matrices and bold lower case letters to denote row vectors. For a matrix A"(a ) , we de"ne vec(A) as the following pq-dimensional row vector: I J N"O vec(A)"(a , a ,2, a , ,a ,a , ,a ). (2) O\ 2 N\ N\ 2 N\ O\ Conversely, if a is a pq-dimensional vector as given in Eq. (2), we de"ne *ec\(a)"A to be the p;q matrix N O whose (k,l)th entry is a . For example, let a"(a , a , a , a , a , a ), then we have I J a a a a a . vec\(a)" a a and vec\(a)" a a a a a If p"q, we de"ne diag(A) as the following p-dimensional row vector, whose entries are the main diagonal elements of A.
diag(A)"(a , a ,2, a ). N\ N\ A denotes the conjugate transpose of A and A is the transpose of A. A*B is the matrix multiplication of two matrices A and B. A B is the tensor product [10] of A and B.
S. Qiu et al. / Signal Processing 77 (1999) 159}170
161
The r;r identity matrix is denoted as I . Following the notation from Davis' classical book on circulant P matrices [8], F denotes the Fourier matrix of order r [8]. F is de"ned [8, p. 32] as follows: P P 1 F" (s )P\ , P (r H I H I where s "uHI with u "e\p P and i"(!1. As it is known, F is a unitary matrix: H I P F\"F . P P It is also known that the tensor products I F and F I are unitary matrices, for any number r and r . P P P P In addition, we recall the de"nitions for time translation operator ¹ and modulation operator M I J [12}14,17]. For two integers k,l and a signal x, the actions of ¹ and M on x are de"ned as follows: I J ¹ (x)"(x(N!k), x(N#1!k),2, x(N!1), x(0), x(1),2,x(N!k!1)), (3) I M (x)"(x(0), uJx(1), uJx(2),2, u,\Jx(N!1)), where u"u "e\p ,. (4) J , 3. Gabor basic matrix 3.1. Matrix formulation For a Gabor window u and lattice parameters a and b, we de"ne the Gabor basic matrix as the following block matrix form: u Mu @ $ u MI @\@ $ M ¹ u K@ L? $
GAB(u, a, b)"
.
(5)
¹ u ?\? M¹ u @ ?\? $ M
¹ u @I \@ ? \?
Clearly, GAB(u, a, b) is an NI ;N matrix (NI "N/ab). M ¹ u is the (nbI #m#1)th row of GAB(u, a, b) for K@ L? n"0,1,2,a !1 and m"0,1,2,bI !1. Let c be the following NI -dimensional row vector: c"(c , c ,2, c I ,2, c ,c , ,c ). @\ ?\ ? \ 2 ? \ @I \ Clearly, c is the (nbI #m#1)th entry, where n"0,1,2,a !1 and m"0,1,2,bI !1. L K With the above notation, Eq. (1) is equivalent to the following matrix formulation: x"c*GAB(u, a, b).
(6)
162
S. Qiu et al. / Signal Processing 77 (1999) 159}170
3.2. Block circulant structure of Gabor basic matrices Throughout the rest of the paper, we always consider the critical sampling cases: ab"N. Then a "b and bI "a. The Gabor basic matrix GAB(u, a, b) is an N;N matrix. Theorem 1. Let V be the xrst a;N block matrix,
u
Mu @ $
V"
.
MI u @\@ For l"0,1,2,b!1, let V be the following a;a square matrix: J u(la) u(la#1) u(la#a!1) 2 wJ?u(la) wJ?>u(la#1) wJ?>?\u(la#a!1) 2 V" J $ $ $ $ wJ??\u(la)
wJ?>?\u(la#1)
2
wJ?>?\?\u(la#a!1)
where w"e\p @ ,"e\p ?. Then V is in the following block matrix form: V"(V V 2V ). @\ The Gabor basic matrix is formulated in the block-circulant form [8]:
V V GAB(u, a, b)" @\ $
V
2
$
2 \
V
V
V
V
(8)
Proof. By the de"nition of modulation operation and bI "a, u
Mu @ $
V"
u(0)
u(1)
u(0)
u@u(1)
$
$
2 \
u?\@u(1)
2
"
M u u(0) ?\@ where u"u "e\p ,. , Since w"e\p @ ,"u@,
V"
(7)
V @\ V @\ . $
2
,
u(N!1)
2
u,\@u(N!1) $
u,\?\@u(N!1)
u(0)
u(1)
2
u(N!1)
u(0)
wu(1)
w,\u(N!1)
$
$
2 \
u(0)
w?\u(1)
2
w,\?\u(N!1)
$
This implies Eq. (7). Eq. (8) can be easily derived from Eqs. (5) and (8). 䊐
.
,
S. Qiu et al. / Signal Processing 77 (1999) 159}170
163
Theorem 2. Let U"F I . Then GAB(u, a, b) is unitarily block-diagonalized by U: there exist a;a block @ ? matrices B for k"0,1,2, b!1 such that I U * GAB(u, a, b) * U "B,
(9)
where
B B
B"
.
\ B
@\
Proof. The proof is similar to the one presented in [12]. Let p be the basic circulant matrix of order b [8]. By Theorem 1, GAB(u, a, b)"I V #p V #2#p@\ V . @ @\
(10)
Using the properties of matrix}tensor products [10], @\ @\ U * GAB(u, a, b)*U " U*(pI V )*U " (F I )*(pI V )*(F I ) I @ ? I @ ? I I @\ @\ " (F I )*(pI V )*(F I )" (F *pI*F ) (I *V *I ) @ ? I @ ? @ @ ? I ? I I @\ " XI V . I I The last equality follows from [8] and
1
X I"
uI \
For k"0,1,2, b!1, set @\ B " uIJV . I J J
with u"ep @, for k"0,1,2, b!1.
uI@\
(11)
Eq. (9) follows immediately. 䊐 Next, we explore the structure of the block matrices B (k"0,1,2,b!1). We show that each block matrix I B is the product of the Fourier matrix F and a diagonal matrix D . I ? I
164
S. Qiu et al. / Signal Processing 77 (1999) 159}170
Theorem 3. For k"0,1,2,b!1, let D be the following diagonal matrices of order a: I
@\uIJu(la) J
D" I
@\uIJu(la#1) J
,
\
(12)
@\uIJu(la#a!1) J
where u"ep @. Then B "(aF *D . I ? I
(13)
Proof. By Eq. (11) and Theorem 1, we see that
@\uIJu(la) J @\uIJwJ?u(la) J $
@\uIJu(la#1) J @\uIJwJ?>u(la#1) J $
@\uIJwJ??\u(la) J
@\uIJwJ?>?\u(la#1) J
B" I
2 2 $ 2
@\uIJu(la#a!1) J @\uIJwJ?>?\u(la#a!1) J $
where w"e\p ? and ep @. Using w?"1, we derive that
@\uIJu(la) J @\uIJu(la) B " J I $ @\uIJu(la) J
@\uIJu(la#1) J @\uIJwu(la#1) J $
1
2
1
1
w
w?\
$
$
2 $
1
w?\
2
w?\?\
*
@\uIJu(la) J
2 $
@\uIJw?\u(la#1) J
1
"
2
$
@\uIJu(la#1) J
\
Corollary 1. Write U"F I and UI "I F . Then @ ? @ ? U*(GAB(u, a, b))*U "(a*UI *D
2
@\uIJwJ?>?\?\u(la#a!1) J
@\uIJu(la#a!1) J @\uIJw?\u(la#a!1) J $
@\uIJw?\?\u(la#a!1) J
"(aF *D . ? I
@\uIJu(la#a!1) J
䊐
,
S. Qiu et al. / Signal Processing 77 (1999) 159}170
165
and GAB(u, a, b)"(a*U *UI *D*U, where
D D
D"
.
\
D @\
Corollary 1 is straightforward from Theorems 2 and 3. Now let us form the row vectors d from the diagonal, the diagonal matrices D (k"0,1,2,b!1): I I d "diag (D ). Then we build the following b;a matrix D : I I B
d
d
. $
D" B
(14)
d @\
The next theorem shows how to calculate D by performing FFT. Hence, the diagonal elements of D can be B computed simply by FFT. Theorem 4. D "(bF **ec\(u). B @ @ ? Proof. By Eq. (12) and u"ep @,
@\ @\ @\ uIJu(la), uIJu(la#1),2, uIJu(la#a!1) . J J J
d" I Thus,
@\u(la) J @\uJu(la) J $
@\u(la#1) J @\uJu(la#1) J $
@\uJ@\u(la) J
@\uJ@\u(la#1) J
D" B
"
1
1
1
u
$
$
2 $
1
u@\
2
"(bF **ec\(u). @ @ ?
1
2
u@\ $
u@\@\ 䊐
*
2 2 $ 2
@\u(la#a!1) J @\uJu(la#a!1) J $
@\uJ@\u(la#a!1) J
u(0)
u(1)
2
u(a!1)
u(a)
u(a#1)
u(a#a!1)
$
$
2 $
u(N!a)
u(N!a#1)
2
u(N!1)
$
166
S. Qiu et al. / Signal Processing 77 (1999) 159}170
As a consequence, we formulate a simple condition for (u, a, b) to generate a Gabor frame. By Theorem 3 and Corollary 1, +M ¹ u, is a Gabor frame if and only if all the diagonal elements D are not zeros. By K@ L? Theorem 4, all the elements of D are not zeros. B Corollary 2. For a Gabor window u and lattice parameters a and b, a necessary and suzcient condition for the Gabor family +M ¹ u, to be a Gabor frame is that there will be no zero elements in F **ec\(u). K@ L? @ @ ? 4. Algorithms In this section, we present the algorithms for both Gabor analysis and Gabor synthesis. Let u be a Gabor window function. a and b are lattice parameters satisfying ab"N. Both algorithms for Gabor analysis and Gabor synthesis have complexity O(N log N), which perform the same as the fast Fourier transform (FFT) does. In addition, we present an algorithm for computing the dual Gabor window u . 4.1. Gabor analysis We present an algorithm to compute the Gabor coe$cients. The algorithm is easy to derive from Theorems 1}4. By Eq. (6) and Corollary 1, it is easy to see that 1 c" *x*U *D\*UI *U. (a Algorithm 1. For a signal x, the following steps are required to compute the Gabor coezcients. E Calculate x "(1/(a)*x*(F I ) with complexity aO(b log b)"O(N log b) via FFT. @ ? E Build the b;a matrix vec\(u) and calculate D "(bF **ec\(u) with complexity aO(b log b)"O(N log b) @ ? B @ @ ? via FFT. E Calculate the point-wise division x of vector x and the vector *ec(D ) with complexity O(N). B E Calculate x "x *(I F ) via FFT with complexity bO(a log a)"O(N log a). @ ? E Calculate Gabor coezcients c"x *(F I ) via FFT with complexity aO(b log b)"O(N log b). @ ? The total complexity for Gabor analysis is, therefore, O(N log N). 4.2. Gabor synthesis The algorithm for Gabor synthesis follows immediately from Theorems 1}4. By Eq. (6) and Corollary 1, x"(a*c*U *(UI )*D*U. Algorithm 2. Given Gabor coezcients C"(c ) . The Gabor synthesis is performed via the following L K ?"@ procedure: E Build the b;a matrix *ec\(u) and calculate D "(bF *vec\(u) with complexity aO(b log b)"O(N log b) @ ? B @ @ ? via FFT. E Form the vector c"*ec(C) and calculate c "(a*c*(F I ) with complexity aO(b log b)"O(N log b) via @ ? FFT, where c"*ec(C).
S. Qiu et al. / Signal Processing 77 (1999) 159}170
167
E Calculate c "c *(I F ) with complexity bO(a log a)"O(N log a) via FFT. @ ? E Calculate the point-wise multiplication c of vector c and vector *ec(D ) with complexity O(N). B E Calculate x"c *(F I ) with complexity aO(b log b)"O(N log b) via FFT. @ ? The total complexity for Gabor synthesis is O(N log N). 4.3. Dual Gabor window Algorithm 3. The following procedures are used to compute the dual Gabor window function u : E Build the b;a matrix *ec\(u) and calculate D "(bF **ec\(u) with complexity aO(b log b)"O(N log b) @ @ ? @ ? B via FFT. E Build the complex conjugate D and calculate the point-wise inverse u of vec(D )with complexity O(N). B B E Calculate the dual Gabor window u "(1/(a)*u *(F I ) with complexity aO(b log b)"O(N log b) via FFT. @ ? The total complexity for computing the dual Gabor window is O(N log b). Proof. The algorithm follows from the Wexler}Raz bi-orthogonality [21]: u *[GAB(u, a, b)] "e, where e"(1,0,2,0) is the N-dimensional unit vector. By Corollary 1, 1 u " *e*U *UI *(D )\*U. (a This leads to the algorithm. 䊐
5. Numerical simulations In this section, we present some of the numerical simulation results. We demonstrate the results computed by Algorithms 1}3. The Gabor window function and a Chirp signal is shown in Fig. 1. The signal length is N"256. In Fig. 2, we choose the lattice parameters a"16 and b"16. In Fig. 3, the lattice parameters are a"32 and b"8. In Fig. 4, the lattice parameters are taken as a"8 and b"32. The Gabor coe$cients and the dual Gabor u are illustrated. The Gabor reconstruction relative errors are of order O(10\).
6. Conclusion We observed the block-circulant structure of the Gabor basic matrices. We showed that the Gabor basic matrices are block- diagonalized by the unitary matrix U"F I , the tensor product of the Fourier matrix @ ? of order b and the identity matrix of order a. Then it is proved that the block-diagonal matrix is decomposed as the matrix product of I F and a diagonal matrix D. The diagonal elements of D are determined by u via @ ? FFT. The algorithms are proposed for Gabor analysis and Gabor synthesis and these algorithms are performed with complexity O(N log N).
168
S. Qiu et al. / Signal Processing 77 (1999) 159}170
Fig. 1. A Gabor window function u and a Chirp signal.
Fig. 2. Gabor coe$cients and dual Gabor window u .
S. Qiu et al. / Signal Processing 77 (1999) 159}170
Fig. 3. Gabor coe$cients and dual Gabor window u .
Fig. 4. Gabor coe$cients and dual Gabor window u .
169
170
S. Qiu et al. / Signal Processing 77 (1999) 159}170
References [1] L. Auslander, I.C. Gertner, R. Tolimieri, The discrete Zak transform application to time}frequency analysis and synthesis of nonstationary signals, IEEE Trans. Signal Process. 39 (4) (1991) 825}835. [2] R. Balart, Matrix reformulation of the Gabor transform, Opt. Eng. 31 (4) (June 1992) 1235}1242. [3] M.J. Bastiaans, The expansion of an optical signal into a discrete set of Gaussian beams, Optik 57 (1) (1980) 95}102. [4] M.J. Bastiaans, Gabor's signal expansion and its relation to sampling of the sliding-window spectrum, in: R.J. Marks II (Ed.), Advanced Topics in Shannon Sampling and Interpolation Theory, Springer, New York, 1993, pp. 1}35. [5] J.J. Benedetto, Gabor representations and wavelets, AMS Contemporary Mathematics Series, Vol. 91, AMS, Providence, RI, 1989, pp. 9}27. [6] I. Daubechies, The wavelet transform, time}frequency localization and signal analysis, IEEE Trans. Inform. Theory 36 (1989) 961}1005. [7] J.G. Daugman, Uncertainty relation for resolution in space, spatial frequency, and orientation optimized by two-dimensional visual cortical "lters, J. Opt. Soc. Amer. 2 (7) (1985) 1160}1169. [8] P.J. Davis, Circulant Matrices, A Wiley-Interscience Publication, New York, 1979. [9] T. Ebrahimi, M. Kunt, Image compression by Gabor expansion, Opt. Eng. 30 (7) (1991) 873}880. [10] P. Lancaster, M. Tismenetsky, The Theory of Matrices, Academic Press, Orlando, 1985. [11] R.S. Orr, The order of computation for "nite discrete Gabor transform, IEEE Trans. Signal Process. 41 (2) (1993) 122}130. [12] S. Qiu, Block-circulant Gabor-matrix structure and Gabor transforms, Opt. Eng. 34 (10) (October 1995) 2872}2878. [13] S. Qiu, Generalized dual Gabor atoms and best approximations by Gabor family, Signal Processing 49 (3) (1996) 167}186. [14] S. Qiu, Discrete Gabor transforms: the Gabor}Gram matrix approach, J. Fourier Anal. Appl. 4 (1998) 1}17. [15] S. Qian, D. Chen, Discrete Gabor transform, IEEE Trans. Signal Process. 41 (7) (1993) 2429}2438. [16] S. Qian, K. Chen, S. Li, Optimal biorthogonal sequence for "nite discrete-time Gabor expansion, Signal Processing 27 (2) (May 1992) 177}185. [17] S. Qiu, H.G. Feichtinger, Discrete Gabor structure and optimal representation, IEEE Trans. Signal Process. 43 (10) (1995) 2259}2268. [18] D. Stewart, L. Potter, S.C. Ahalt, Computationally attractive real Gabor transform, IEEE Trans. Signal Process. 43 (1) (1995) 77}84. [19] J. Yao, Complete Gabor transformation for signal representation, IEEE Trans. Image Process. 2 (2) (1993) 152}159. [20] J. Yao, P. Krolak, C. Steele, The generalized Gabor transform, IEEE Trans. Image Process. 4 (7) (1995) 978}988. [21] J. Wexler, S. Raz, Discrete Gabor expansion, Signal Processing 21 (1990) 207}220.
Discrete Gabor transforms with complexity O(NlogN)夽 Sigang Qiu *, Feng Zhou, Phyllis E. Crandall Ambient Technologies Inc., NCDC, 110 Horizon Drive, Raleigh, NC 27615, USA PowerPC Microprocessor Validation, Department RYMA, Bldg. 667, Research Triangle Park, NC 27709, USA Los Alamos National Laboratory, P.O. Box 1663, ACL, MS B287, Los Alamos, NM 87545, USA Received 27 August 1997; received in revised form 11 February 1999
Abstract We develop fast computations of discrete Gabor transforms. We present algorithms for computations of both Gabor analysis and Gabor synthesis via the fast Fourier transform. The algorithms are performed with complexity O(N log N). In addition, we propose an algorithm for computing the dual Gabor wavelets with complexity less than O(N log N). 1999 Elsevier Science B.V. All rights reserved. Zusammenfassung Wir entwickeln schnelle Berechnungsmethoden fuK r diskrete Gabor-Transformationen. Wir praK sentieren Algorithmen fuK r die Berechnung der Gabor-Analyse sowie der Gabor-Synthese mittels der Schnellen Fouriertransformation (FFT). Diese Algorithmen weisen eine KomplexitaK t von O(N log N) auf. Weiters schlagen wir einen Algorithmus zur Berechnung des dualen Gabor-Wavelets mit einer KomplexitaK t von weniger als O(N log N) vor. 1999 Elsevier Science B.V. All rights reserved. Re2 sume2 Nous proposons dans cet article des meH thodes de calcul rapide de transformations de Gabor discre% tes. Nous preH sentons des algorithmes de calcul a` la fois pour l'analyse et la synthe`se de Gabor via la transformation de Fourier rapide (en anglais FFT). Ces algorithmes ont une complexite` O(N log N). De plus, nous proposons un algorithme pour le calcul d'ondelettes de Gabor duales ayant une complexiteH infeH rieure a` O(N log N). 1999 Elsevier Science B.V. All rights reserved. Keywords: Discrete Gabor transforms (DGT); Gabor analysis; Gabor synthesis; Fast Fourier transform
1. Introduction The interest in Gabor transforms arises from the fact that Gabor basis functions has good time}frequency localization [6,7] with appropriate Gabor window functions. Given a Gabor window u and a pair of lattice parameters (a,b), there are two basic and fundamental questions. 夽 The authors were with the University of Connecticut when this paper was submitted. The referee's insightful comments on the paper are acknowledged. * Corresponding author. Tel.: #1-919-870-5060; fax: #1-919-870-9898. E-mail address: [email protected] (S. Qiu)
0165-1684/99/$ - see front matter 1999 Elsevier Science B.V. All rights reserved. PII: S 0 1 6 5 - 1 6 8 4 ( 9 9 ) 0 0 0 3 0 - 4
160
S. Qiu et al. / Signal Processing 77 (1999) 159}170
E Gabor decomposition (analysis): For a signal x, how can we e$ciently "nd the Gabor coe$cients +c , L K such that x is represented (or best approximated) by the following form: ? \ @I \ x" c u , (1) L K L K L K where +u "M ¹ u, are Gabor basis (or Gabor elementary) functions. ¹ and M are the time L K K@ L? L K L? K@ translation operator and modulation operator. We call u a Gabor window function. a "N/a and bI "N/b are referred to as the dual lattice parameters. E Gabor reconstruction (synthesis): Given a set of coe$cients +c ,, how can we e$ciently determine L K a signal x with the form of Eq. (1). As Orr [11] observed, there are several approaches to the above questions. From Orr's summary [11, Table 1], we may say that matrix methods are comparatively expansive; however, there are many papers that describe matrix approaches to the topic. Variant matrix techniques have proven to be e!ective and clear. We refer the interested reader to some of the relevant papers [1}5,9,12}20] and references cited in each. We consider the critical Gabor transform in this paper. The aim is to show that both the Gabor analysis and synthesis can be performed with the same complexity as FFT. The paper is organized in the following way. Relevant notation is listed in Section 2. In Section 3 we investigate the structures of Gabor basic matrix. Then we are able to propose relevant algorithms in Section 4. We show that the relevant computations are performed with a couple of FFT. A necessary and su$cient condition for (u, a, b) to generate a Gabor frame is formulated. Numerical simulations are given in Section 5. We make our conclusions in Section 6.
2. Notation Signals are considered as row vectors in ",, and we identify signals with N-periodic complex sequences. That means, for a signal x"(x( j)),\3",, we treat x as an N-periodic sequence such that x( j#kN)"x( j) H for j"0,1,2,N!1, k39. We use bold upper case letters to denote matrices and bold lower case letters to denote row vectors. For a matrix A"(a ) , we de"ne vec(A) as the following pq-dimensional row vector: I J N"O vec(A)"(a , a ,2, a , ,a ,a , ,a ). (2) O\ 2 N\ N\ 2 N\ O\ Conversely, if a is a pq-dimensional vector as given in Eq. (2), we de"ne *ec\(a)"A to be the p;q matrix N O whose (k,l)th entry is a . For example, let a"(a , a , a , a , a , a ), then we have I J a a a a a . vec\(a)" a a and vec\(a)" a a a a a If p"q, we de"ne diag(A) as the following p-dimensional row vector, whose entries are the main diagonal elements of A.
diag(A)"(a , a ,2, a ). N\ N\ A denotes the conjugate transpose of A and A is the transpose of A. A*B is the matrix multiplication of two matrices A and B. A B is the tensor product [10] of A and B.
S. Qiu et al. / Signal Processing 77 (1999) 159}170
161
The r;r identity matrix is denoted as I . Following the notation from Davis' classical book on circulant P matrices [8], F denotes the Fourier matrix of order r [8]. F is de"ned [8, p. 32] as follows: P P 1 F" (s )P\ , P (r H I H I where s "uHI with u "e\p P and i"(!1. As it is known, F is a unitary matrix: H I P F\"F . P P It is also known that the tensor products I F and F I are unitary matrices, for any number r and r . P P P P In addition, we recall the de"nitions for time translation operator ¹ and modulation operator M I J [12}14,17]. For two integers k,l and a signal x, the actions of ¹ and M on x are de"ned as follows: I J ¹ (x)"(x(N!k), x(N#1!k),2, x(N!1), x(0), x(1),2,x(N!k!1)), (3) I M (x)"(x(0), uJx(1), uJx(2),2, u,\Jx(N!1)), where u"u "e\p ,. (4) J , 3. Gabor basic matrix 3.1. Matrix formulation For a Gabor window u and lattice parameters a and b, we de"ne the Gabor basic matrix as the following block matrix form: u Mu @ $ u MI @\@ $ M ¹ u K@ L? $
GAB(u, a, b)"
.
(5)
¹ u ?\? M¹ u @ ?\? $ M
¹ u @I \@ ? \?
Clearly, GAB(u, a, b) is an NI ;N matrix (NI "N/ab). M ¹ u is the (nbI #m#1)th row of GAB(u, a, b) for K@ L? n"0,1,2,a !1 and m"0,1,2,bI !1. Let c be the following NI -dimensional row vector: c"(c , c ,2, c I ,2, c ,c , ,c ). @\ ?\ ? \ 2 ? \ @I \ Clearly, c is the (nbI #m#1)th entry, where n"0,1,2,a !1 and m"0,1,2,bI !1. L K With the above notation, Eq. (1) is equivalent to the following matrix formulation: x"c*GAB(u, a, b).
(6)
162
S. Qiu et al. / Signal Processing 77 (1999) 159}170
3.2. Block circulant structure of Gabor basic matrices Throughout the rest of the paper, we always consider the critical sampling cases: ab"N. Then a "b and bI "a. The Gabor basic matrix GAB(u, a, b) is an N;N matrix. Theorem 1. Let V be the xrst a;N block matrix,
u
Mu @ $
V"
.
MI u @\@ For l"0,1,2,b!1, let V be the following a;a square matrix: J u(la) u(la#1) u(la#a!1) 2 wJ?u(la) wJ?>u(la#1) wJ?>?\u(la#a!1) 2 V" J $ $ $ $ wJ??\u(la)
wJ?>?\u(la#1)
2
wJ?>?\?\u(la#a!1)
where w"e\p @ ,"e\p ?. Then V is in the following block matrix form: V"(V V 2V ). @\ The Gabor basic matrix is formulated in the block-circulant form [8]:
V V GAB(u, a, b)" @\ $
V
2
$
2 \
V
V
V
V
(8)
Proof. By the de"nition of modulation operation and bI "a, u
Mu @ $
V"
u(0)
u(1)
u(0)
u@u(1)
$
$
2 \
u?\@u(1)
2
"
M u u(0) ?\@ where u"u "e\p ,. , Since w"e\p @ ,"u@,
V"
(7)
V @\ V @\ . $
2
,
u(N!1)
2
u,\@u(N!1) $
u,\?\@u(N!1)
u(0)
u(1)
2
u(N!1)
u(0)
wu(1)
w,\u(N!1)
$
$
2 \
u(0)
w?\u(1)
2
w,\?\u(N!1)
$
This implies Eq. (7). Eq. (8) can be easily derived from Eqs. (5) and (8). 䊐
.
,
S. Qiu et al. / Signal Processing 77 (1999) 159}170
163
Theorem 2. Let U"F I . Then GAB(u, a, b) is unitarily block-diagonalized by U: there exist a;a block @ ? matrices B for k"0,1,2, b!1 such that I U * GAB(u, a, b) * U "B,
(9)
where
B B
B"
.
\ B
@\
Proof. The proof is similar to the one presented in [12]. Let p be the basic circulant matrix of order b [8]. By Theorem 1, GAB(u, a, b)"I V #p V #2#p@\ V . @ @\
(10)
Using the properties of matrix}tensor products [10], @\ @\ U * GAB(u, a, b)*U " U*(pI V )*U " (F I )*(pI V )*(F I ) I @ ? I @ ? I I @\ @\ " (F I )*(pI V )*(F I )" (F *pI*F ) (I *V *I ) @ ? I @ ? @ @ ? I ? I I @\ " XI V . I I The last equality follows from [8] and
1
X I"
uI \
For k"0,1,2, b!1, set @\ B " uIJV . I J J
with u"ep @, for k"0,1,2, b!1.
uI@\
(11)
Eq. (9) follows immediately. 䊐 Next, we explore the structure of the block matrices B (k"0,1,2,b!1). We show that each block matrix I B is the product of the Fourier matrix F and a diagonal matrix D . I ? I
164
S. Qiu et al. / Signal Processing 77 (1999) 159}170
Theorem 3. For k"0,1,2,b!1, let D be the following diagonal matrices of order a: I
@\uIJu(la) J
D" I
@\uIJu(la#1) J
,
\
(12)
@\uIJu(la#a!1) J
where u"ep @. Then B "(aF *D . I ? I
(13)
Proof. By Eq. (11) and Theorem 1, we see that
@\uIJu(la) J @\uIJwJ?u(la) J $
@\uIJu(la#1) J @\uIJwJ?>u(la#1) J $
@\uIJwJ??\u(la) J
@\uIJwJ?>?\u(la#1) J
B" I
2 2 $ 2
@\uIJu(la#a!1) J @\uIJwJ?>?\u(la#a!1) J $
where w"e\p ? and ep @. Using w?"1, we derive that
@\uIJu(la) J @\uIJu(la) B " J I $ @\uIJu(la) J
@\uIJu(la#1) J @\uIJwu(la#1) J $
1
2
1
1
w
w?\
$
$
2 $
1
w?\
2
w?\?\
*
@\uIJu(la) J
2 $
@\uIJw?\u(la#1) J
1
"
2
$
@\uIJu(la#1) J
\
Corollary 1. Write U"F I and UI "I F . Then @ ? @ ? U*(GAB(u, a, b))*U "(a*UI *D
2
@\uIJwJ?>?\?\u(la#a!1) J
@\uIJu(la#a!1) J @\uIJw?\u(la#a!1) J $
@\uIJw?\?\u(la#a!1) J
"(aF *D . ? I
@\uIJu(la#a!1) J
䊐
,
S. Qiu et al. / Signal Processing 77 (1999) 159}170
165
and GAB(u, a, b)"(a*U *UI *D*U, where
D D
D"
.
\
D @\
Corollary 1 is straightforward from Theorems 2 and 3. Now let us form the row vectors d from the diagonal, the diagonal matrices D (k"0,1,2,b!1): I I d "diag (D ). Then we build the following b;a matrix D : I I B
d
d
. $
D" B
(14)
d @\
The next theorem shows how to calculate D by performing FFT. Hence, the diagonal elements of D can be B computed simply by FFT. Theorem 4. D "(bF **ec\(u). B @ @ ? Proof. By Eq. (12) and u"ep @,
@\ @\ @\ uIJu(la), uIJu(la#1),2, uIJu(la#a!1) . J J J
d" I Thus,
@\u(la) J @\uJu(la) J $
@\u(la#1) J @\uJu(la#1) J $
@\uJ@\u(la) J
@\uJ@\u(la#1) J
D" B
"
1
1
1
u
$
$
2 $
1
u@\
2
"(bF **ec\(u). @ @ ?
1
2
u@\ $
u@\@\ 䊐
*
2 2 $ 2
@\u(la#a!1) J @\uJu(la#a!1) J $
@\uJ@\u(la#a!1) J
u(0)
u(1)
2
u(a!1)
u(a)
u(a#1)
u(a#a!1)
$
$
2 $
u(N!a)
u(N!a#1)
2
u(N!1)
$
166
S. Qiu et al. / Signal Processing 77 (1999) 159}170
As a consequence, we formulate a simple condition for (u, a, b) to generate a Gabor frame. By Theorem 3 and Corollary 1, +M ¹ u, is a Gabor frame if and only if all the diagonal elements D are not zeros. By K@ L? Theorem 4, all the elements of D are not zeros. B Corollary 2. For a Gabor window u and lattice parameters a and b, a necessary and suzcient condition for the Gabor family +M ¹ u, to be a Gabor frame is that there will be no zero elements in F **ec\(u). K@ L? @ @ ? 4. Algorithms In this section, we present the algorithms for both Gabor analysis and Gabor synthesis. Let u be a Gabor window function. a and b are lattice parameters satisfying ab"N. Both algorithms for Gabor analysis and Gabor synthesis have complexity O(N log N), which perform the same as the fast Fourier transform (FFT) does. In addition, we present an algorithm for computing the dual Gabor window u . 4.1. Gabor analysis We present an algorithm to compute the Gabor coe$cients. The algorithm is easy to derive from Theorems 1}4. By Eq. (6) and Corollary 1, it is easy to see that 1 c" *x*U *D\*UI *U. (a Algorithm 1. For a signal x, the following steps are required to compute the Gabor coezcients. E Calculate x "(1/(a)*x*(F I ) with complexity aO(b log b)"O(N log b) via FFT. @ ? E Build the b;a matrix vec\(u) and calculate D "(bF **ec\(u) with complexity aO(b log b)"O(N log b) @ ? B @ @ ? via FFT. E Calculate the point-wise division x of vector x and the vector *ec(D ) with complexity O(N). B E Calculate x "x *(I F ) via FFT with complexity bO(a log a)"O(N log a). @ ? E Calculate Gabor coezcients c"x *(F I ) via FFT with complexity aO(b log b)"O(N log b). @ ? The total complexity for Gabor analysis is, therefore, O(N log N). 4.2. Gabor synthesis The algorithm for Gabor synthesis follows immediately from Theorems 1}4. By Eq. (6) and Corollary 1, x"(a*c*U *(UI )*D*U. Algorithm 2. Given Gabor coezcients C"(c ) . The Gabor synthesis is performed via the following L K ?"@ procedure: E Build the b;a matrix *ec\(u) and calculate D "(bF *vec\(u) with complexity aO(b log b)"O(N log b) @ ? B @ @ ? via FFT. E Form the vector c"*ec(C) and calculate c "(a*c*(F I ) with complexity aO(b log b)"O(N log b) via @ ? FFT, where c"*ec(C).
S. Qiu et al. / Signal Processing 77 (1999) 159}170
167
E Calculate c "c *(I F ) with complexity bO(a log a)"O(N log a) via FFT. @ ? E Calculate the point-wise multiplication c of vector c and vector *ec(D ) with complexity O(N). B E Calculate x"c *(F I ) with complexity aO(b log b)"O(N log b) via FFT. @ ? The total complexity for Gabor synthesis is O(N log N). 4.3. Dual Gabor window Algorithm 3. The following procedures are used to compute the dual Gabor window function u : E Build the b;a matrix *ec\(u) and calculate D "(bF **ec\(u) with complexity aO(b log b)"O(N log b) @ @ ? @ ? B via FFT. E Build the complex conjugate D and calculate the point-wise inverse u of vec(D )with complexity O(N). B B E Calculate the dual Gabor window u "(1/(a)*u *(F I ) with complexity aO(b log b)"O(N log b) via FFT. @ ? The total complexity for computing the dual Gabor window is O(N log b). Proof. The algorithm follows from the Wexler}Raz bi-orthogonality [21]: u *[GAB(u, a, b)] "e, where e"(1,0,2,0) is the N-dimensional unit vector. By Corollary 1, 1 u " *e*U *UI *(D )\*U. (a This leads to the algorithm. 䊐
5. Numerical simulations In this section, we present some of the numerical simulation results. We demonstrate the results computed by Algorithms 1}3. The Gabor window function and a Chirp signal is shown in Fig. 1. The signal length is N"256. In Fig. 2, we choose the lattice parameters a"16 and b"16. In Fig. 3, the lattice parameters are a"32 and b"8. In Fig. 4, the lattice parameters are taken as a"8 and b"32. The Gabor coe$cients and the dual Gabor u are illustrated. The Gabor reconstruction relative errors are of order O(10\).
6. Conclusion We observed the block-circulant structure of the Gabor basic matrices. We showed that the Gabor basic matrices are block- diagonalized by the unitary matrix U"F I , the tensor product of the Fourier matrix @ ? of order b and the identity matrix of order a. Then it is proved that the block-diagonal matrix is decomposed as the matrix product of I F and a diagonal matrix D. The diagonal elements of D are determined by u via @ ? FFT. The algorithms are proposed for Gabor analysis and Gabor synthesis and these algorithms are performed with complexity O(N log N).
168
S. Qiu et al. / Signal Processing 77 (1999) 159}170
Fig. 1. A Gabor window function u and a Chirp signal.
Fig. 2. Gabor coe$cients and dual Gabor window u .
S. Qiu et al. / Signal Processing 77 (1999) 159}170
Fig. 3. Gabor coe$cients and dual Gabor window u .
Fig. 4. Gabor coe$cients and dual Gabor window u .
169
170
S. Qiu et al. / Signal Processing 77 (1999) 159}170
References [1] L. Auslander, I.C. Gertner, R. Tolimieri, The discrete Zak transform application to time}frequency analysis and synthesis of nonstationary signals, IEEE Trans. Signal Process. 39 (4) (1991) 825}835. [2] R. Balart, Matrix reformulation of the Gabor transform, Opt. Eng. 31 (4) (June 1992) 1235}1242. [3] M.J. Bastiaans, The expansion of an optical signal into a discrete set of Gaussian beams, Optik 57 (1) (1980) 95}102. [4] M.J. Bastiaans, Gabor's signal expansion and its relation to sampling of the sliding-window spectrum, in: R.J. Marks II (Ed.), Advanced Topics in Shannon Sampling and Interpolation Theory, Springer, New York, 1993, pp. 1}35. [5] J.J. Benedetto, Gabor representations and wavelets, AMS Contemporary Mathematics Series, Vol. 91, AMS, Providence, RI, 1989, pp. 9}27. [6] I. Daubechies, The wavelet transform, time}frequency localization and signal analysis, IEEE Trans. Inform. Theory 36 (1989) 961}1005. [7] J.G. Daugman, Uncertainty relation for resolution in space, spatial frequency, and orientation optimized by two-dimensional visual cortical "lters, J. Opt. Soc. Amer. 2 (7) (1985) 1160}1169. [8] P.J. Davis, Circulant Matrices, A Wiley-Interscience Publication, New York, 1979. [9] T. Ebrahimi, M. Kunt, Image compression by Gabor expansion, Opt. Eng. 30 (7) (1991) 873}880. [10] P. Lancaster, M. Tismenetsky, The Theory of Matrices, Academic Press, Orlando, 1985. [11] R.S. Orr, The order of computation for "nite discrete Gabor transform, IEEE Trans. Signal Process. 41 (2) (1993) 122}130. [12] S. Qiu, Block-circulant Gabor-matrix structure and Gabor transforms, Opt. Eng. 34 (10) (October 1995) 2872}2878. [13] S. Qiu, Generalized dual Gabor atoms and best approximations by Gabor family, Signal Processing 49 (3) (1996) 167}186. [14] S. Qiu, Discrete Gabor transforms: the Gabor}Gram matrix approach, J. Fourier Anal. Appl. 4 (1998) 1}17. [15] S. Qian, D. Chen, Discrete Gabor transform, IEEE Trans. Signal Process. 41 (7) (1993) 2429}2438. [16] S. Qian, K. Chen, S. Li, Optimal biorthogonal sequence for "nite discrete-time Gabor expansion, Signal Processing 27 (2) (May 1992) 177}185. [17] S. Qiu, H.G. Feichtinger, Discrete Gabor structure and optimal representation, IEEE Trans. Signal Process. 43 (10) (1995) 2259}2268. [18] D. Stewart, L. Potter, S.C. Ahalt, Computationally attractive real Gabor transform, IEEE Trans. Signal Process. 43 (1) (1995) 77}84. [19] J. Yao, Complete Gabor transformation for signal representation, IEEE Trans. Image Process. 2 (2) (1993) 152}159. [20] J. Yao, P. Krolak, C. Steele, The generalized Gabor transform, IEEE Trans. Image Process. 4 (7) (1995) 978}988. [21] J. Wexler, S. Raz, Discrete Gabor expansion, Signal Processing 21 (1990) 207}220.