- Email: [email protected]
Generalized dual Gabor atoms and best approximations by Gabor family
SIGNAL
PROCESSING Signal Processing
49 (1996)
167-186
Generalized dual Gabor atoms and best approximations by Gabor family Sigang Qiu* Department of Mathematics, Uniuersity of Connecticut U-9, 196 Auditorium Road. Storrs, CT 06269-3009, USA Received
13 January
1995; revised
14 December
1995
Abstract Let g be a Gabor window of length N and (a,b) be a pair of lattice constants. u and b are considered as time and frequency gaps for a TF-lattice with N*/ab elements. The corresponding time-frequency shifts of g form a so-called Gabor family. In this paper, we investigate the structural properties of the discrete Gabor transforms. We address the problem of finding the best approximation of a signal x E CN by linear combinations of a Gabor family. We consider critical sampling, oversampling and undersampling, and do not assume that the Gabor family is a flame. For the task we determine the (generalized) dual Gabor atom (GDGA). This amounts to determining the pseudoinverse of the Gabor matrix and can be solved by the conjugate-gradient (CG) algorithm with O(N) complexity for fixed lattice constants (a, b). We provide an easy practical criterion for checking whether a Gabor triple (g, a, b) generates a Gabor frame or not. We propose an efficient algorithm for estimating the Gabor frame bounds and an algorithm for determining tight Gabor atoms. Zusammenfassung Es seien g ein Gaborfenster der Lange N und (a, b) ein Parr von Gitterkonstanten; a und b werden als Zeit- und Frequenzhicken fir ein Zeit-Frequenz-Gitter mit N’/ub Elementen betrachtet. Die entsprechenden Zeit-Frequenz-Verschiebungen von g bilden eine sogenannte Gabor-Familie. In diesem Beitrag diskntieren wir die strukturellen Eigenschaften der diskreten Gabortransformationen. Wir sprechen das Problem an, wie die beste Approximation eines Signals x E CN durch lineare ijberlagerungen einer Gaborfamilie zu finden ist. Dabei betrachten wir kritische Abtastung wie Unter- und ijberlastung und nehmen nicht an, dalj die Gabor-Familie einen Rahmen bildet. Zu diesem Zweck bestimmen wir das (verallgemeinerte) duale Gaboratom (GDGA). Das ftihrt dazu, dag die Pseudoinverse der Gabormatrix bestimmt werden mug, was fiir feste Gitterkonstanten (a, b) mit dem “Conjugate-Gradient” (CC) Algorithmus mit einer Komplexitit O(N) zu l&en ist. Wir stellen ein einfaches praktisches Kriterium bereit, mit dem zu testen ist, ob ein Gabortripel (g, a, b) einen Gaborrahmen erzeugt oder nicht. Wir schlagen einen effizienten Algorithmus zur Schatzung der Gaborrahmen-Grenzen vor sowie ein Verfahren znr Bestimmung dichter Gaboratome.
Soit g une fen&e de Gabor de longueur N et (a, 6) une paire de constantes de treillis. a et b sont consider& comme etant les intervalles de temps et de frequence pour un treillis TF de N’/ub elements. Les deplacements temps-Mquence correspondants de g forment ce qu’on appelle une famille de Gabor. Dans cet article, nous discutons des proprietes structurelles des transfotmees de Gabor disc&es. Nous considerons le probleme de trouver la meilleure approximation d’un signal x E CN par * Tel.: INT++(860)-486-1280;
fax: INT++(860)-486-4238;
e-mail: [email protected].
0165-1684/96/$15.00 @ 1996 Elsevier Science B.V. All rights reserved PIISOl65-1684(96)00015-l
168
S. Qiu / Signal Processing
49 (1996)
167-186
des combinaisons lineaires d’une famille de Gabor. Nous etudions l’echantillomrage critique, le SW et le sous-tchantillonnage, et ne supposons pas que la famille de Gabor est une trame. Pour cette t&he, nous determinons l’atome de Gabor dual (generalise) (GDGA). Ceci revient a determiner le pseudo-inverse de la matrice de Gabor et peut etre resolu en utilisant l’algorithme du gradient conjugue (GC) de complexite O(N) pour des constantes de treillis fix&es (a, b). Nous donnons un critere pratique simple pour verifier si un triplet de Gabor (g, a, b) genere une trame de Gabor ou pas. Nous proposons un
algorithme efficace pour estimer les limites de la trame de Gabor et un algorithme pour determiner les atomes de Gabor. Keywords: Gabor expansion; (Generalized) dual Gabor atom (GDGA); 9%symmetry
1. Introduction The theory of Gabor expansions has been thoroughly studied by many scientists, mostly via the Zak transform, see e.g. [3] and the references cited there. Since Gabor elementary functions are nonorthogonal, the main difficulty of doing Gabor analysis is how to obtain suitable Gabor coefficients efficiently. As we know, these coefficients can be obtained by short-time Fourier transforms (STFT) with a so-called biorthogonal Gabor analysis window. It is the desire of finding efficient determinations of this dual Gabor window (atom) that has attracted considerable attention from both mathematicians and engineers. Actually, the oversampled and the critically sampled Gabor families which are frames have been studied essentially by the following two methods: - Via the Zak transform and the frame method: We refer to [l-3,7,14,18,27,30]. - Via matrix decompositions: We refer to [9,17,1925,28,29]. The main purpose of this paper is to consider discrete Gabor transforms generally. Specifically for an arbitrary Gabor triple (g,a,b), we present a general method to obtain the associated (generalized) dual Gabor atom (GDGA) g based on the CG(conjugategradient)-method, which can be implemented with O(N) complexity for fixed lattice constants (a, b). For a signal X, using the calculated & we are able to determine the best approximation of x by linear combinations of the Gabor family. We propose an algorithm to determine the tight Gabor atom gt based on the binomial series. In Section 2, we list some notation which have been widely used in our previous papers [21,23] and some preliminary results about discrete Gabor expansions. Section 3, the main part of this paper, will discuss the Rotated(g)-symmetric structure of the Gabor
matrix in detail and will give an easy way of computing the Gabor frame bounds. Section 3 also contains the corresponding numerical algorithms. In Section 4 we provide some numerical results. We draw our conclusions in Section 5.
2. Notation and preliminaries We use essentially the notation introduced in [21]. For the paper to be self-contained, we recall the main definitions. We use i= J--f to denote the imaginary unit and j, k, 1, a, b are always nonnegative integers. is considered as an N(i) A signal x=(x(k)}::: periodic row vector in the complex space CN. We use either Xk or x(k) to denote the (k + I)thcoordinateofasignalx(k=O,l,...,N-1). The inner product of two signals x,y E CN is
(ii)
&Y> = CL;’ xo’M.0 Signals are-treated as row vectors. For a matrix A E CNxN, it induces a linear mapping (operator) A : x H x * A. The adjoint operator corresponds to x c+ x * A’, where A’ is the conjugate transpose (or Hermitian transpose) of A. We use AT to denote the transpose of A. The time translation operator Tk corresponds to the cyclic rotation:
Tk(x041, .a. %-I) 3
:= (X-k,Xl__k, ... . x-l,&
.. . .
XN-k-1).
The modulation operator (or the frequency translation operator) is given by M
(x0,x1,
:=
*..,
Q-1)
. XN--Id--I) (X0,X~W,X~W2, ..)
>
S. Qiu/Signal
Processing 49 (1996)
where w=e (iii)
-2xi’fN for 1 E R.
(vi)
The following rotation operator (introduced in [21,23]) with the rotation number a acting on vectors is also used for the cyclic rotation:
167-186
169
satisfy that ub
r&(x,
T,g
a) := (XN_-a,XNma+l, ... , XN-13x0,
If y is a column
. . . . XN-a-1).
T(i- I jag Mb!? MbTaY
vector, we define rot@, a) :=
[rotdy’, a)]‘. The matrix rotation on matrix A with a rotation number a is obtained by rotating all the column vectors of A. Precisely if A E CPxq and A =
‘-W,
a,b) =
(2) MbT,,- I )a.4
(AI, AZ, ... , Aq), then
MmbTnag
rotm(A, a) := (rot(A,, a), rot(A2, a), ... . rot(A,, a)). For any number k = 1,2, ... . we define inductively rotmk(A, a) = rotm (rot&-‘(A, a),~). It
is
easy
to
check
that
rotm’(A, a) =
rotm(A, ku). (iv)
(v)
n=O m=O
where the Gabor elementary functions gnm are given as g,,,,, = M&,Tnag for n = 0,1, . . . ,6 - 1 and m = O,l,..., b - 1. The collection {gnm}nfl is called a Gabor family. {c,,}~~ are usually referred to as the Gabor coefficients. The Gabor operator S associated to (g, a, b) is defined as
=
d-l
6-l
C
Cb,
n=O
m=O
Then it is easy to see the following sentation for the Gabor operator: Sx=x*S
For a signal g E CN, we call (g,u, b) a Gabor triple if both a and b are positive integers and divisors of N. We call g a Gabor atom. The constants a and b are called time and frequency gaps. They form a pair of lattice constants. We call G = N/u and b”= N/b the dual lattice constants. A Gabor expansion of a signal x with respect to a Gabor triple (g,u, b) is a series expansion of the form
Sx
MC&_, y,T(d-I jag
snm)snm,
for
x E P.
(1)
Unless the contrary is explicitly stated, we assume from now on that lattice constants (a, b)
matrix repre-
forxECN,
where S = [GAH(g, a, b)]’ * [GAE%(g,a, b)] is called the Gabor matrix. In general, for any two Gabor triples (gi,u, b) and (gz,u, b), we call the matrix product [GAB(gi, a, b)]’ * [GAEi(gz, a, b)] the Gabor matrix associated to (gi,u, b) and (92, a, b). For a Gabor triple (g, a, b), the associated Gabor operator S and the Gabor matrix S may not be invertible. However, pinv(S) is uniquely defined, where pinv(S) denotes the pseudoinverse [26] (or the Moore-Penrose inverse [ 16, pp. 432-4401) matrix of S. Then we define the associated generalized dual Gabor atom (GDGA) as g := g * pinv(S). Since the Gabor matrix S is positive semidefinite [l 1, 16,261, the square root S1’2 is well defined [16, pp. 180-1841. We define the generalized tight Gabor atom (GTGA) as gt = g * pinv(S1j2). If a Gabor triple (g, a, b) generates a Gabor frame, the associated Gabor operator S is the usual Gabor frame operator (cf. [12, pp. 635-6361 or [6, pp. 58591). In this case, pinv(S) = inv(S) is the usual inverse of S. The associated GDGA g and GTGA gt are the dual and tight Gabor atoms in the usual sense (cf.
170
S. Qiu/ Signal Processing 49 (1996) 167-186
[27], [6] or [12]). In this case, we call g and gt the usual dual and tight Gabor atoms. The following proposition is used to obtain the best approximation of a signal x E CN by linear combinations of a Gabor family {gnm}n,nz. This is equivalent to determine the orthogonal projection of x onto the Gabor space: Span{g,,}, where Span{g,,} is the linear space spanned by the Gabor family {gnm}nfl [16, pp. 83-871. Proposition 1. Let J be the GDGA associated to (g, a, b). For a signal x E CN, the orthogonal projection of x onto Span{g,,} can be represented as
This yields that g,, E Span{g,,}. Then using (5), we derive that MmbTna(g * pinv(G)) * {G * pinv(G))
gnnm =
= MmbTna{g * pinv(G) * G} * pinv(G) = (MmbT,,g) * pinv( G) = gnm* pinv( G). This proves Eq. (6). Hence, i-1
6-1
7
x(X,
d-1 Km)gnm
6-1
= c
n=O m=O
x(X9
gnm * pinv(G))gnm
n=O m=O d-l
=
6-l
C X(X * pinv(G),gnm)gnm n=O m=O
In particular, ifx E Span{g,,}, 8-l
= {x * pinv( G)} * G.
n=O m=O
n=O m=O
6-l
8-l
(3)
then
We conclude from Eqs. (5) and (6) that if x E SPan{gnm1
6-1 X =
{X * pinv(G)}
* G
n=O m=O
n=O m=O
=
d-l
6-1
71
Fl(X,
n=O
m=O
inm>gnm,
and Let us sketch the proof of Proposition 1. For further details, we refer to [24]. Let G = [GAB(g,a, b)]’ * [GAEl(g,a, b)] be the Gabor matrix associated to (g, a, b). It is clear that G and thus pinv(G) are Hermitian symmetric. For a signal x E Span{g,,}, we see that x is in the row space of the Gabor matrix G (cf. [16, pp. 78-801, [26, pp. 137-1401). Hence, we have (see [16, pp. 428-4381, [26, pp. 140-1411) {x * G} * pinv(G) = {x * pinv(G)} * G =x
for X E Span{g,,}.
(5)
x = {x * G} * pinv( G)} 8-l = C
z-1 = T;7
=
(MmbTnag)
i.e. &,, = gym.
gnm)
(6)
In fact, by the commutative property of the Gabor operator with both T, and Mb, it is easy to see that
7:
tx3 gnm)!7nm.
Eq. (4) is proved. For the general case (3), since P(x) E Span{g,,}, using (4), we see that d-l = 71
6-1 x(P(X)9
= M,bT,,,
(g * @NW) (g
*
I?-1 6-1
Similarly,
pinv(G*) * G)
=MmbTna (g * pinv(G’))
* G.
gnm)tjnm
=
d-l
6-l
7:
FItxy
gnm)Lm
n=O m=O
= {x * G} * pinv(G) = {x * pinv(G)} * G
n=O m=O CTnnm= MmbTnn
* PWG))
6-l
n=O m=O
* PWG),
(gnm
n=O m=O
P(X) (g *WV(G))
x(X,
n=O m=O
In order to prove Eqs. (3) and (4), we need to show first that MmbTna
5-l
i-1
5-l
n=O m=O
S. Qiu / Signal Processing 49 (1996)
Eq. (3) is obtained. The sketch of the proof of Proposition 1 is complete.
167-186
171
More precisely
Remark. Since Eq. (6) holds without the assumption
if B is given (by Theorem A) as
Cl,1
0
..
0
c2,2
’
0 0
..
that {gnm I,,, is a frame, Eq. (6) is an extended version of the result presented in [4, Theorem 3.1 I], [6, pp.
0
831, [7, PP. 9711.
0 0
c6+l,1
0
Ck2.2
... ...
. cCZ,Cl 0 0
...
.. B=
3. Mathematical setup
0
As we have seen in [21,23], the Gabor matrix S associated to Gabor triples (gi,a,b) and (gz,a,b) is completely determined by its first N x a block matrix B. All the other a”- 1 block matrices are obtained from B by proper rotations. We state these descriptions in the following theorem. and Feichtinger [21,23]). Let (gl,a, b) and (gz,a, b) be two Gabor triples, S = [GAB(g,, a, b)]’ * [GAE!(gz, a, b)] be the associated Gabor matrix. Assume that B is thejirst N x a block matrix qf S. Then S is of the following block matrix ,form:
Theorem A (Qiu
C(b-l)6+l,l
0
0
forr=a”-1.
(7)
For k = 1,2, ... , a” - 1, rotn#(B, a) is exactly the (k + I)-th block matrix of size N x a. Moreover, the possibly nonzero entries of B are distributed only in the I-th diagonals for 1 = 0,416, f26,. . . , f(b - 1)b”. We call B the first block corresponding to Gabor triples (gi, a, b) and (92, a, b). If gi = g2 = g, we call B the first block matrix associated to (g, a, b). Extracting all the zero-entries of B, we obtain a compressed “nonzero”-block matrix B of size b x a. The Ith row-vector is exactly the Ith diagonal of B for 1 = 0, &, ,.. , (b - 1)J. We call B the (first) “nonzero”-block matrix corresponding to Gabor triples (gl,a, b) and (gz,a, 6). When gi = g2 = g, we call B the “nonzero”-block matrix associated to (g, a, b). By Theorem A and the definition of rotm, every (sub-)diagonal is a-periodic. We consider B as u-periodic in column-subscripts and b-periodic in row-subscripts.
0
:::
0
. ..
0
1)&+2,2
0
then the b x a “nonzero”-block
B=
CI,I
c2,2
ClT+1,1 .
%+2,2
. ctY+n,a
...
‘(b-
I$++a,a
matrix B is
...
C,,ll
...
C6+a,o
..
:
t'(b-I)&,,,
a), ... , rotm’(B, a)]
...
..
i S = [B, robn(B,
‘(b-
0
: ‘(b-1)6+2,2
“.
1.
C(b-I)i++a,ai
Since B is of small size b x a, B can be easily built up. We give a sufficient condition with r? for a Gabor triple (g, a, b) to generate a Gabor frame.
Corollary 1 (Sufficient condition). Let B = (Ui,j)bXa be the “nonzero”-block matrix associated to a Gabor triple (g, a, b), and assume that for any 1
then (g, a, b) generates a Gabor frame.
Proof. Condition
(8) implies that the Gabor matrix S associated to (g, a, b) is strictly diagonal dominant. Therefore, S is nonsingular and (g, a, b) does generate a frame. 0 In the following, we show that if gi = g2 = g then the Gabor matrix associated to (g, a, b) has the rotated-symmetry (%!-symmetry) structure.
S.
172
Qiu/ Signal Processing 49 (1996) 167-186
Theorem 1 (W-symmetry structure). Assume that gr = g2 = g, and let
In general, for Y = 46 + 1, s = (b - q)b” + 1 with q = 1, 2, ...) bo - 1, we can show similarly that
S = [B, rotm(B), ... , rotmr(B)],
di = d;&
wherer=a”-
1,
(9)
be the Gabor matrix associated to (g, a, b). Writing b,, =
i(b+l)
ifbisodd,
i(b + 2)
if b is even,
then the Gabor matrix S has at most abo independent nonzero elements. They are from the upper part of the N x a block matrix B. Proof. Using Theorem A, we show first that for r = qb”+l,s=(b-q)b”+lwithq=1,2, ....ba-1.the sth subdiagonal elements of B are completely determined by the rth subdiagonal elements of B via the rotation and the conjugation. Conversely, we can determine the rth subdiagonal from the sth subdiagonal. Without loss of generality, we take Y = 6 + 1 and s = (b - 1)b” + 1. Noticing that the diagonal elements of the block matrix are of periodicity N = bb”,for k = 1, 2, ... , a, the general elements of the sth subdiagonal can be written as
4 = S((bl)K+k),k i-1
=bCT,,g((bn=O
l)h+k-
l)Tnag(k-
1)
6-l
=bxT,,g(k-b-
l)T,,g(k-
1)
l)Tnag(k-
1)
n=O
for k = 1, 2, ... , a.
This means that the sth subdiagonal of B is compketely obtained via a rotation with the rotation number qb (or the remainder of qb divided by a) and the conjugation of the rth subdiagonal of B.
Hence, all the independent entries of S are located in the rth subdiagonal of B, for Y = qb + 1 with q = 0,l , ..*, bo - 1. The total number of independent entries are at most abo. The proof is complete. 0 Remark. Since abo c iab, the B-symmetry of the Gabor matrix has reduced the computational cost by half. In general, if gr # 92, the Gabor matrix corresponding to (gk, a, b) for k = 1,2 does not have the Wsymmetry. Using the introduced “nonzero”-block matrix, we rewrite Theorem 1 as the following equivalent form. Theorem 2. Let B be the “nonzero”-block matrix associated to (g,a, b). Assume that 8j is the jth row vector of B for j = 2, ... , bo. Then the (b - j + 2)-th row vector 8b_j+2 of 8 can be written as Bb-j+2 = rot(Bj, mod((j - l)b, a)), where mod(Cj - 1)b”, a) is the remainder of Cj - 1)b” divided by a. In particular if the signal length N is divisible by ab (especially for the critical cases), then
d-l
=bxT,,g(k-6n=O
bb_j+2 = 6
for j = 2, ... , bo.
If the Gabor atom is a real-valued signal, then
= d;_6
fib--j+2 = rot(LTj,mod(~ - l$,
Therefore, di = d;_6
for k = 1, 2, ... , a.
s = (b - 1)b” f 1 -th subdia( > gonal elements of B can be obtained via a suitable rotation and the conjugation from the (r = 6 + 1)-th subdiagonal elements of B. This yields that the
a))
for j = 2, ... , bo. Theorem 3. Let S be the Gabor matrix corresponding to a Gabor triple (g,a, b). Then for any k = O,l, ... . Sk is W-symmetric. More generally, if p(z) is any real polynomial, then p(S) also has the W-symmetry.
S. Qiu / Signal Processing 49 (1996)
For the proof of Theorem 3, we need the following interesting proposition.
Proposition 2. Let S k/2 be the linear operator realized by the matrix Ski2 for k = 0,1,2,. . . . Then Ski2 fork = 0,1,2,... commute with both the time translation operator T, and the frequency modulation operator Mb.
167-186
173
Proof of Theorem 3. Set h = S(k-1)/2g, then the Gabor operator Sh associated to (h, a, b) is defined as 6-l
&x =
d-l
c
x(x,h,,)h,,
m=O
n=O
X E 6?
On the other hand, by Proposition Sk1 &k-i)/=
2,
(S (S(k-‘Y=X))
Proof. If k = 0, the fact is trivial. We only need to consider the cases for k > 1. For k = 1, S is the associated Gabor operator. It is known that S commutes with both T, and Mb. Since S is positive semidefinite, Ski2 and Ski2 are well defined. Assume that U is the unitary matrix diagonalizing S, then there is a diagonal matrix D = diag(ai,az,...,aN) such that
for
6-I z-1
x ~(S(k-1)‘2x,Ynm~Ynm
= SW 1Y=
m=On=O 6-l
=
d-l
x
~(+,s(k-‘)~=gn,)S(k-‘):=g,m
m=O
n=O
6-l d-l
S=
U’*D*U,
(10)
where c113 1x23 . . . UN 2 0, diag (ai, ~12,.. . , UN) denotes the diagonal matrix with aj (j = 1,2,. . . , N) being its diagonal entries. From Eq. (lo), it is easy to see that
=
c
~(~,hmJhm.
m=O
n=O
Thus,
This implies that = (J’
e
diag
Sk = [GAJ3(h, a, b)]’ * [GAB(h, a, b)].
u:‘~,ui12,.. . , OLD* U.
(11)
>
By the polynomial interpolation theory [13], there is an interpolation polynomial L(x) such that L(uj) = $I’,
The l-norm
for j = 1,2,.. .,N.
Ski2 = U’* diag(L(ai),L(a&...,L(a,v))
[ll,
13, 161 of a matrix A = (a;,j),lx,
is defined as
From Eq. ( 11 ), we can easily check that
=U’*L(D)*U=L(S).
Therefore, Sk is &!-symmetric by Theorem 1. For any real polynomial p(z), p(S) also has the .C&symmetry. 0
n
*U
0
Then by the fact that S commutes with both T, and Mb, we deduce that Ski2 commutes with both the timetranslation operator T, and the modulation operator Mb. 0 With the same argument given above, we obtain the following corollary.
Corollary 2. Zf we define ,Ykj2 = (pinv(S))k’2 for k = 1,2,..., then Sk12for k = 0, fl, 4~2,. . . commute with both T, and Mb.
IPIll = ,y2mCladiI. i=l
It is easy to compute, especially for a small-size matrix. We show that the l-norm of the small “nonzero”-block matrix can be used to estimate the Gabor frame upper bounds. Suppose (g, a, b) generates a Gabor frame, S is the associated Gabor frame operator. Let 8 be the dual Gabor atom, 8 and B> be the b x a “nonzero”-block matrices associated to (g, a, b) and (6, a, b), respectively. Then
Theorem 4. (Gabor frame bounds).
S. Qiu/ Signal Processing 49 (1996)
174
Proof. Assume that S and S are the Gabor matrices associated to (g, a, b) and (g, II,b), respectively, then IlSll =
IPII
and IIS-‘II = IF’II
= IISII-
Since S = [B, rotm(B, a), ... . rotm’(B, a)] and B is obtained from B by filling in some proper zero-entries, it is clear that
Ilrota% a)111= IlBlll = 11~111. Hence,
ll~lll = IPIll = IPIll. Similarly,
II% = II&Ill,
i.e.,llS-‘Il~ = Il&ll~.
Because the norm (spectral norm) of S is majorized (see [13; 16, pp. 358-3621) by the l-norm of S, we see that
A
= ~~IIsIl~llsII1 W’ll1
= 11411.
Therefore,
167-186
(2) If u = 1, let 8, be the “nonzero”-block corresponding to (F(g), b,l), then the best Gubor upper bound is obtained us (14~ Ilo3/N, Proof. Since the Gabor matrix associated to (g,u, 1) is a nonnegative diagonal matrix [23, Corollary 81, the first result is obvious, while the second result is deduced from [23, Theorem 31. 0 Algorithm 1 (Nonzero-block and Gabor matrix). Let (g,u, b) be a Gabor triple and S be the associated Gabor matrix. Then we can construct S with only z(& ) M ;Nb multiplications and via proper rotations and the conjugation by the following steps: Step 1. We calculate the first bo row-vectors of the b x a “nonzero”-block matrix l?. The total computation needs only (ubo )C = Nbo M +Nb multiplications. The (k, j)-th general entry of the “nonzero’‘-block B can be computed as d-l ~~~=gCTnog((k-1)5+j-l)TnagV-l) n=O
for j = 1, 2, ... , a and k = 1, 2, ... , bo. Step 2. We obtain the remaining b - bo row-vectors in the lower part of B. If Bi denotes the jth u-elements row vector of B for j = 1,2, ... , b, then for j = 27 ... 3 bo,
Remark. (1) Let &,,, and ;Iminbe the maximal and minimal eigenvalues of the Gabor matrix S, then Aand &in are the best upper and lower Gabor frame bounds. Tables 1 and 2 show that the approximate Gabor frame bounds ub := l[~llt andLb := l/l(B>ll, calculated from Theorem 4 approximate well to I,, and kin. (2) With the simple and efficient formula presented for the Gabor frame bounds, the condition number of the Gabor matrix can be estimated as
which can be used to check the quality (that is, the stability of the associated Gabor expansion and synthesis) of a Gabor signal. Corollary 3. (1) If b = 1, then the “nonzero”-block corresponding to (g, a, 1) is a vector B = (tj)yE, E Cu. llB]l, =maxy=, tj is the best Gubor upper bound.
Bb-j+Z = rot(Bj,mod((‘j - l)b,
u)).
Step 3. Filling in some extra zero-entries into 8, we get the first N x a block matrix B of the Gabor matrix. Step 4. By the rotation operator rotm with the rotation number a, we can successively obtain the other 2 - 1 block matrices of the same size N x a. The Gabor matrix S is built up.
Algorithm 2 (Gabor matrix-vector multiplication). Let (g,a, b) be a Gabor triple, S and S be the associated Gabor operator and Gabor matrix, B be the associated “nonzero”-block matrix of size b x a and x be an arbitrary signal. Then the general elements of y = Sn = x * S := (Yo, Yl, ... . Yv_1) can be computed from B = (ck,l),,, by
S. Qiu/ Signal Processing 49 (1996)
167-186
175
Table 1 Comparison between the Gabor frame bounds ub (Lb,) and A,, (in,,“). We use ub (Lb) to denote the approximate Gabor upper (lower) bound obtained by Theorem 4. A,, (A,,) stands for the best Gabor frame upper (lower) bound. The Gabor atom shown in Fig. 3 is of a typical Gaussian shape Lattice
(ib
(a. b)
( uap*rw.
(12, 12) (IO, 10) (8, 8) (8, 6)
597.50260416615 717.00312961905 896.26049300375 897.82568699835
(0, b)
KDZ =
(12, 12)
570.215206083184 66.364320832716 11.420117443573 2.910703099535
(10, 10) (8, 8) (8, 6)
*mm 1
ub/Lb
Gxact
597.50260416597 717.00312961882 896.26049300366 897.82568699835
K =
I .047853832802 10.804045315651 78.480847279568 308.456636179028
1
1.047853832801 10.8040453 15657 78.4808447279583 308.456636179024 ub
-
Lb
ub
+
Lb
knmax/&b
570.215602083367 66.364320832654 11.420117443569 2.910703099535
0.91961152311793 0.78133544866770 0.5433 1293495694 0.26092308033431
0.99649869507642 0.9703 1069303042 0.83897092687839 0.48858301203230
Table 2 Comparison between the Gabor frame bounds ub (Lb) and A,,,= (&,). (ib (Lb) and &,, same meanings as those in Table 1. The exponential Gabor atom is shown in Fig. 4 Lattice
ub
ta,b)
( uapprox.
1max 1
(10, 10) (8, 8) (8, 6)
1922.76264640444 2306.32832577053 2882.42656567521 2886.37655929048
(a, b)
Km =
(12,
(12,
12)
12)
(10, 10) (8, 8) (8, 6)
Ub/Lb
813.739718040677 375.122874805591 168.773569670759 47.896900894645
(&act
Lb )
1919.88421461691 2303.70222154611 2880.82779694440 2886.37655929050
K =
~Lpprox
(&in)
have the
bin 1
2.36287181733500 6.14819431357204 17.07866090228586 60.26228222238036
vexact
)
2.37681424145194 7.11090203642794 17.08182182956302 60.26228222237835 Lb
a-v%
‘h
-
a+cb
ub
+Lb
&xxx/bin
807.755263803913 323.967649919045 168.648744009180 47.896900894647
wherej=(r-l)a+s--lforr=1,2,...,Zand s = 1, 2, , a. The derivation of this algorithm is similar to the proof of Algorithm 2 in [23]. With the precalculated “nonzero’‘-block 8, we only need Nb multiplications to perform Algorithm 2. We should also notice that if S is the Gabor matrix corresponding to (gk, a, b) for k = 1,2, i.e., S = [GAB(gi, a, b)]’ * [GAEl(gz,a, b)], then Algorithm 2 does also hold. In practical applications, we do not need to build the N x N Gabor matrix S. We only need to work with the small b x a “nonzero”block matrix B by Algorithm 2.
0.93226337332680 0.90180718516839 0.85705390609887 0.74749893572402
0.99754522830333 0.99468258876562 0.9882 1960329939 0.95909761225340
We now present the CG-method for obtaining the GDGA y” associated to a Gabor triple (g,a, b). The standard CG-method is stated in the following [ 11, p. 5231 Algorithm A (Conjugate gradient [ 111). If A is positive definite, then we have the following iterative algorithm for computing x so that x * A = y, for a vector y: Starting with x0 = 0, Q = y, p1 = t-0the iterative step is as follows, for k > 1, _ Set/-$=Oand/?k=*; _ Set pi = ?i, and pk = yk-_l + Pkpk-I;
S. Qiu/ Signal Processing
176
- setuk=&; - set - set
xk
=
xk__l
•,-
‘tkpk;
rk
=
rk-1
-
akpk
*A.
It is known that N is the upper bound of the number of iteration steps and XN is the exact solution to the linear system x * A = y if the algorithm is performed in infinite precision. In this case, we say that the CGalgorithm is convergent after N iterations. For a Gabor triple (g, a, b), let S be the associated Gabor matrix to (g,a,b). We show in the following theorem that the CG-method can always be applied and be convergent after at most ab iterations to obtain the associated GDGA j by solving the linear equation g * S = g. Using the Gabor matrix-vector multiplication algorithm 2, the product Pk *S for each CG-cycle can be implemented in only O(N) operations.
49 (1996)
167-186
at most ub-dimensional subspaces of CNxN and CN, respectively. If g is real-valued, by Theorem 1 and Theorem 3, for any 1 = 0,1, ... . the (r = (b - q)6 + 1)-th subdiagonal of B(l) can be obtained by (s = q6 + 1 )-th subdiagonal of B(‘) via only the rotation with the rotation number qb”.This yields that the dimension of 9 and consequently the dimension of Pg over R are at most ubo. 0 Proof of Algorithm 3 . Step 1. We show that the CG-procedure works and converges after at most ub iterative steps. Following the arguments in [ 11, pp. 5 16-5231, we only need to show that Ukis always defined. Equivalently, we need to show the following: ifrk-_l #O,
thenpk*S*pL#O.
(13)
Algorithm 3 (Generalized dual Gabor atom (GDGA)). Let S be the Gabor matrix corresponding to a Gabor triple (g, a, b). If we replace both A and y by S and g in Algorithm A, respectively, then the CG-algorithm is convergent after s (s < ub) iterations and the resulting x, = 6 is the solution for g * S = g. Moreover, g = x, is exactly the GDGA corresponding to (g, a, b), that is, i = g * pinv(S). In particular if g is real-valued, the algorithm is convergent after at most ubs iterations.
In fact, suppose that there exist pk and r&t such that Pk * S * pi = 0 but rk_1 # 0, then Pk * [GAB(g, a, b)]’ * [GAB(g, a, b)] * p; = 0, which yields that Pk * [GAB(g,u, b)]’ = 0. Since S = [GAB(g,u, b)]’ * [GAB(g, a, b)], g is in the row space of GAB(g, a, b). Inductively, it is easy to check that pk is in the row space of GAB(g,u, b). We deduce that (see [ 11, pp. 50]), pk = 0. This implies that
In order to derive Algorithm 3, we need the following lemma.
and
But rk_l~i_~ = 0 holds (see [ll, Theorem 10.2.2]), r,&_l = 0 follows. So we get a contradiction and Eq. (13) is true. Therefore, if rk_1 # 0, then we always have that pk * S * pk # 0. By [l 1, Theorem 10.2.2 and Theorem 10.2.31, after k - 1 steps of the CG-algorithm, the residuals
P’x = Span{x * s’ : I = 0, 1, ... }, for x E CN.
ra,rl, ... . rk__l E Span{g,g*S,
Then both 9’ and 9x (especially, ~9~) are at most ubdimensional subspaces of CN XNand CN, respectively. If g is a real-valued Gubor atom, then the dimensions of both 9’ und Ps over R are at most ubo.
are mutually orthogonal. By Lemma 1, the dimension of the subspace Span{g * S’, for 1 = 0,1, ... , } c CN is at most ub. Hence, there exists s (s Gab) such that the (s + 1 )-th residual r, = 0, which yields that x, * S = g. Write # = x,, then g * S = g. If g is a real-valued signal, the algorithm is convergent after ubo steps by Lemma 1, Step 2. We show that the solution g to the linear system 6 * S = g based on the CG-method is exactly the GDGA y = g * pinv(S).
Lemma 1. Write B = Span{& : I= O,l, ...}
Proof. For I = 0, 1, ... , let B(‘) be the first N x a block matrix of S’. By Theorem 3, S’ are all completely determined by B (‘1 which have at most ub nonzero entries. Because the nonzero entries of B(l) for 1 = 0, 1, ... ) are in the same positions, both 9 and 9’~ are
F-k-1 + bkpk--l
= 0.
... . g*Sk-t}
S. Qiu/ Signal Processing 49 (1996) 167-186
By the CG-iteration steps, we see that there exists a (real) polynomial p(t) with at most ab degree such that $=g*p(S)
and
g*S=g.
We show that p(S) = pinv(S) and therefore 6 = y * pinv(S). Since S and p(S) are Hermitian and commute with each other, with the definition of the pseudoinverse [ 11, 16,261, we only need to show that s * p(S)
* s = s.
Let S be the Gabor Since S commutes and the frequency 0, 1, ... , N - 1, we
(14) operator associated with (g, a, b). with both the time translation T, modulation A&, for any k, 1 = have
Yk.1* P(S) * S = S (P(S)&l)
= S MM,,,
177
bound for the number of CG-iteration steps is ab. Therefore, the total operations for the computation of the GDGA i is about 2ab(bN f4N). If the lattice constants (a, b) are fixed, then we obtain an O(N) CG-algorithm. - Jacobi-method: When the Gabor atom is given as a typical Gaussian function (for example, the Gabor atom shown in Fig. 1, the associated Gabor matrix is (strictly) diagonal dominant. We use the Jacobimethod [ 11, 131 instead to improve the implementation. In this case, the computations and decompositions of Gabor matrices can also be realized by the associated “nonzero”-block matrices. Fig. 1 shows the comparison of convergent rates between the Jacobi-method and the CG-method. - PCG-method: In practice, for a given Gabor triple (g, a, b), the Gabor atom g is well localized. Even if the associated Gabor matrix S = does not have good diagonal dominance, we can use the preconditioned conjugate gradient (PCG) algorithm [ 111 with the preconditioner (sj,k)N
= (S(P(S))& = (@ * S),,
= (9 * P(S) * S),,, gk,l.
=
Thus,
P = diag(sl,l,S*,2,...,SNsV)
gk,l * p(S) * S = gk,l
for
k, 1 = 0, 1, ... . N - 1. (15)
By the definition of the Gabor Eq. (2)), Eq. (15) yields that GAB(g,a,b)
xN
basic
* p(S) * S = GAB(g,a,b).
matrix
(see
(16)
Since S = [GAB(g, a, b)]’ * [GAB(g, a, b)], Eq. (14) follows immediately if we multiply both sides of Eq. (16) by [GAIQ,a, b)l’. The proof is complete.
0
Remark. - Complexity analysis: Algorithm
3 is always applied with Algorithm 2. By Algorithm 2, we do not need to build up the full N x N Gabor matrix S, we only need to compute and save the associated b x a “nonzero”-block matrix l?. This largely saves the storage space. Using Algorithm 2, the calculation of pk * S for every CG-iteration requires only Nb multiplications and Nb additions. Since L? can be precalculated with only bo(a + N) multiplications and boN addtions in total, Algorithm 3 can be performed fast. By Algorithm 3, the upper
to improve Algorithm 3. Equivalently, if we set SO=S*P-’ and go=g*P-‘, then we need to solve the equation S * SO = go by the CGalgorithm. Since SO = S * P-’ = (sj,k/sk,k)N: N. Without the assumption that (g, b”,G) generates a Gabor frame, we have extended the algorithms proposed in [23,24].
Corollary 4 (Gabor frame condition). Let i be the associated GDGA to a Gabor triple (g,a, b), i& = (ukj)hXa be the “nonzero”-block matrix associated
S. Qiu / Signal Processing 49 (1996)
178
167-186
Comparison of Convergent Rates of Jacobi- and CG-methods
Fig. 1. Convergent-rate comparison of the Jacobi- and the CG-methods; a Gabor atom and the associated dual Gabor atom. The Gabot atom is well concentrated and the associated Gabor matrix is diagonal dominant. The signal length is 128 and lattice constants are (8, 8).
with the Gabor triples (g,a, b) and (@,a, b). Then (g, a, b) generates a Gabor frame if and only if
‘k,j = “5, =
1 0
ifk=l, ifk # 1,
(17)
for k = 1,..., a and j = l,..., bo. Proof. Let Ss = [GAT3(g,a, b)]’ * [GAB(&a,b)]. Then (g, a, b) generates a Gabor frame if and only if SO is an identity matrix. The condition (17) follows. Now suppose that condition (17) is satisfied, we show that (g,a, b) generates a Gabor frame. Since there exists a real polynomial p(t) such that g = g * pinv(S) = g * p(S), following the similar argument to the proof of Theorem 3, we can verify that the Gabor matrix So can be identified with p(S) * S. By Theorem 3, SO satisfies the g-symmetry and the condition (17) implies that So is an identity matrix. Therefore, (g,a, b) generates a Gabor frame. Cl Remark. For an arbitrary Gabor triple (g,a, b), we can always check whether (g, a, b) generates a Gabor
frame or not via two steps: ( 1) Determine g by Algorithm 3. (2) Test only abo M $ab number of conditions (17). In the following, we show that the GDGA g obtained by Algorithm 3 is the closest one to the Gabor atom g in the sense of Qian and Chen [ 191. Theorem 5 can be seen as an extended discrete version of Proposition 3.2 in the recent work of Janssen [15]. Theorem 5. Assume that g and S are the GDGA and the Gabor matrix associated to a Gabor triple (g,a,b), then 5 = g * pinv(S) is the unique solution of the problem
minimize (I&
- &/Iover
with y * [GAB(g,b,a”)]’ = $e. Proof. Write &? = Span(g,,),
y E CN
(18)
then {gnm}is a frame of 4?. Let S and Sd be the Gabor operators associated with {gnm}in CN and in 4, respectively. Then S& is invertible. It is easy to check that S”fr = SI”4r
S. Qiu/ Signal Processing 49 (1996)
and
179
and
S~x=x*s
forxEd.
$
S,i’ has the following
matrix representation:
Si’x
for x E A.
= x * pinv(S)
gt := g * (pinv(S))‘j2
= g * pinv(S’!2)
is exactly the GTGA corresponding to (g, a, b). For any 1’ satisfying the biorthogonality y * [GAB(g, 6, a”)]’ = $e, we can express gt as d-l 5-l St = c C(gt, n=O m=O
Ynm)Snm.
On the other hand, gt E &? can also be expressed as
n=O
m=O
Since (g,a,b) generates a frame of .&’ and i is the dual Gabor atom associated to (g, a, b) in &‘, the frame property (cf. [6, Proposition 3.2.4; 12, Proposition 2.1.51 or originally [8, Lemma VIII]) yields that C-1
d-l 6-l =
IG II=0m=O
cc
d-l
6-l
=
6-1
= cc
Ikh
n=Om=O
d-1
&m)l’
6-l
Ih Ym)12= y x IbY mm)12
G :x n=O
G&l)12
m=O
n=O
= (59 9)
~llill lbll GIIYII Ilsll.
It follows that g = g * pinv(S) to the problem (18). 0
Hence, 6 = g * pinv(S) = S>‘g is the usual dual Gabor atom associated to (g,a, b) in A. Similar to Janssen’s argument [ 15, p. 171, we can show that @ is the unique solution to the problem (18). In fact, setting gt = S,>“2 g, then g, E _&’and
Ils”ll
167-186
m=O
11;‘.4Yl12~/lYl12>
where y.4 is the orthogonal projection of y onto A. Since y * [GAB(g, 6, Z)]’ = $e implies that (y, g) =
is the unique solution
Remark. By choosing a different initial vector ~0, we can modify Algorithm 3 to obtain desirable biorthogonal Gabor analysis atoms. In fact, let h be a signal, if we intend to obtain a biorthogonal Gabor atom i. which is close to h, then we set x0 = g - h * S and apply Algorithm 3. Assume go is the output, then g = go + h is what we want. Algorithm 4 (Tight Gabor atom). Given a Gabor triple (g, a, 6) which generates a frame, the following procedure is used to compute the associated tight Gabor atom g,: Step 1. By Algorithm 1, we obtain the “nonzero”block matrix B of size b x a and the Gabor frame upper bound ub by Theorem 4. Step 2. We calculate the dual Gabor atom and henceforth the lower Gabor frame bound Lb. then & := YE is associated to Step 3. Set r = &, the Gabor frame operator S, := KS which satisfies that
&~llS~lI
=rllSll<2Uh Lb
+
uh
This yields that
:= p< 1.
III -Srll<~
Step 4. Since /I1 - &(I 67
< 1,
00 s-112 r - c &(I - 0, k=O
where h(h-l)...(h-k+l)
!Xo= 1,
& =
for h = -i Hence,
and k = 1, 2,
abIN,
1 .2...k .
30
s--1/2
llnhd=2-2Rei%% =2(1- Nlhll llsll) ab
=
J;cuk(I
-
Sr)k.
k=O
Therefore, gt = s-g.
the tight Gabor atom can be obtained as
S. Qiu/ Signal Processing 49 (1996)
180
Iteratively, set St(O)= fig,
gt@+l) = J;a,+,
(gp
4. Numerical results
for n 20,
- S,gJ”‘) .
Then we have lim gt@) = gt. n-+oo The error estimation after n iterations is II&)
-
St11
In fact,
II2 J;
Ilgt(n)-&II=
ak(I-Sr)kg
k=n+l
=
&(I-$)“+’ (1
=
/I
gak+n+l k=O
167-186
(I-Sr)kg (1
E/,gll.
Remark. Even if (g,a,b) does not generate a Gabor frame, Algorithm 4 also works. Algorithm 5 (Best approximation). Assume that (g,a,b) is a Gabor triple and i is the associated GDGA, we obtain the best approximation of a signal x by linear combinations of the Gabor family via the following steps: Step 1. We obtain the “nonzero”-block of size b x a corresponding to the Gabor triples (g,a,b) and (&a, b). Step 2. Applying Algorithm 2 and Proposition 1, we obtain xapproX. In particular, if (g,a, b) generates a frame or a signal x is in the Gabor space spun{gnm}, then we can reconstruct a signal x via the above two steps. Remark. In the two-dimensional case, if the Gabor atom is given as a 2-D separable one, then we can calculate the orthogonal projection of an image onto the 2-D Gabor space by using Proposition 1 and the higher dimensional results presented in [23].
This section is devoted to numerical simulations. The performances of the previously proposed algorithms are illustrated. Given an arbitrary Gabor triple, the numerical constructions of the associated GDGA, GTGA and the best approximation of a given signal are shown. All the numerical experiments were carried out using MATLAB 4.0 on a Sun4 Station. As we have demonstrated, the algorithms proposed in this paper are in a very general sense. We do not assume that a given Gabor family is a frame. We do not restrict that the signal length is the power of 2. Signals can be either real valued or complex valued. Gabor atoms (Gabor window functions) can be chosen arbitrarily. Since Gabor atoms with Gaussian types are typical (they are minimizing the uncertainty principle [5]), most of our examples are experimented with Gabor atoms of Gaussian shapes. When a Gabor atom is well concentrated, the associated Gabor matrix is strongly diagonal dominant. We use the Jacobi-method instead of the CG-algorithm to speed up the convergence. In Fig. 1, the Gabor atom is of a Gaussian shape. The signal length is 128 and the lattice constants (a, b) are (8,8). Since Gabor window functions are usually with good localizations, the associated Gabor matrices have some diagonal dominance (even if they are not strictly diagonal dominant). In this case, the simple PCGmethod is employed to improve the ordinary CGalgorithm. Fig. 2 compares the convergent rates of the PCG, CG and the ordinary frame-method. Both PCG and CG algorithms work much faster, while PCG is better than CG. Meanwhile, if (g, a, b) does not generate a frame, the frame method does not work. However, the CG-method can always be applied to obtain the GDGA. Tables 1 and 2 illustrate the comparison between the best Gabor frame bounds and the estimated ones calculated according to Theorem 4. For Table 1, we take a Gabor atom with the Gaussian shape shown in Fig. 3. For Table 2, we choose a Gabor atom with the exponential shape shown in Fig. 4. Lattice constant pairs are (12,12), (10, lo), (8,8) and (8,6) for both Tables 1 and 2. It is clear that the calculated Gabor bounds ub and Lb approximate well to the best Gabor
S. Qiu/ Signal Processing 49 (19%)
167-186
181
frame-method
\ lo-l5 -
‘, CG
I
1o-2o -_ 10” Gabor atom
dual Gabor atom
I
r
I :
1. 0.8. 0.6. 0.4. 0.2. o-
-50
0
50
Fig. 2. Convergent-rate comparison of the PCG, CG and the frame methods. The solid line denotes the convergent rates for the PCG-method; the dash-dotted line is for the CC&method; the dotted line is for the frame method. Gabor atom is of the Gaussian type. Signal length is 120 and lattice constants are (8, 6).
bounds A,,, and &in. From the Gabor upper and lower bounds, we can predetermine the stability of the Gabor syntheses with the associated dual Gabor atoms shown in Figs. 3 and 4. In fact, the data in the column with Ub/Lb = K, [ 1 l] estimate the condition number of the associated Gabor matrices. If U&b is small, then the corresponding Gabor synthesis is stable. The last two columns of Tables 1 and 2 show clearly why the CG-method converges much better than the frame-method [6,7]. In fact, Theorem 10.2.5 in Golub [ 1 l] shows that the CC convergence rate is majorized by
g’.=
VG-6 Jrr;;+&=
Jlcoo-1 &+1’
while the frame-method &,-Lb .--==. ‘=‘Ub+&,
convergence
rate depends on
&.-1 &c + 1
The data in the last two columns of Tables 1 and 2 show that 01 is much smaller than a2. Therefore, the CG-method converges much faster.
Fig. 3 illustrates a Gabor atom and the associated dual and tight Gabor atoms. The signal length N is 480 and the lattice constants are (i) (12,12), (ii) (IO, lo), (iii) (8,8) and (iv) (8,6). For a random signal, the reconstruction relative errors with the dual and tight Gabor atoms are of order 10-15. Fig. 4 shows the exponential Gabor atom and the dual (tight) Gabor atoms corresponding to lattice constants: (i) (12, 12), (ii) (10, lo), (iii) (8, 8) and (iv) (8, 6). The signal length is 480. For a random signal, the reconstruction relative errors with the dual (tight) Gabor atoms are of order 10-14. Next we show some experiments with Gabor triples that do not generate Gabor frames. Given a Gabor triple (g,a,b) and a signal x with the same signal length as g. If x is in the Gabor space Span(g,,), then we get the reconstructed signal of x by Algorithm 5. If x is not in the Gabor space Span{g,,}, then we can only obtain the best approximation xapprox E Span{g,,} of the signal x. Fig. 5 illustrates the Gabor atom and the associated GDGA @and GTGA qt. The signal length N is 1200.
S. Qiu/ Signal Processing 49 (1996)
182
167-186
tight Gabor atom (i)
dual Gabor atom (i)
Gabor atom
I 1 0.5 0 FYIL
J
200
200
400
,
0 -5 L!!!!N
20
3 2 1 0 -1 LLI
5
200
0.06 0.04 0.02 0 0.08 •LI
10 0 200
400
tight Gabor atom (ii)
x IoA dual atom (iv)
x ,03 dual atom (iii)
x to-3 dual atom (ii)
0.02
•LI 200
400
tight Gabor atom (iii) 0.06 0.04
400
tight Gabor atom (iv) 0.04
I
I
0.02
0 200
400
0w 200
400
400
Fig. 3. A Gabor atom of the Gaussian shape and the associated dual (tight) Gabor atoms. The signal length N is 480 and the dual and tight atoms associated with different lattice constants pairs: (i) (12, 12); (ii) (10, 10); (iii) (8, 8) and (iv) (8, 6).
The Gabor atom with respect to the lattice constants (a, b) = (25, 16) does not generate a frame. In Fig. 6 we show (I) the original signal from the Gabor space Span{g,,} and the reconstructed signal, the (relative) reconstruction error is of order 10-14; (II) the original Chirp signal and the best approximate Chirp by the linear combinations of the Gabor family. Fig. 7 shows the original and the best approximate image (Lena) of size 128 x 128 with linear combinations of the 2-D Gabor family. A 2-D Gabor atom and the associated GDGA are shown in Fig. 8. The 2-D lattice constants (a, b) are taken to be a = (16,8) and b = (8,4). Fig. 9 gives another numerical simulation for the 2-D case. The original image given as a rectangular matrix of size 463 x 331 and the best approximate image by the linear combination of the 2-D Gabor family are illustrated. The 2-D Gabor atom of size 480 x 480 with the 2-D lattice constants (a,b) for
a = (30,24) and b = (15,10) Gabor frame.
does not generate a 2-D
5. Conclusion We have developed the W-symmetry feature of Gabor matrices. The B-symmetric property reduces the computational complexity by approximately half, for building up the Gabor matrix, its first N x a block matrix and the associated b x a “nonzero”-block matrix. The algorithms presented in [2 1,231 are improved. In practice, we do not need to construct the full N x N Gabor matrix S, the Gabor matrix-vector multiplications can be realized completely by the small b x a (ab d N) “nonzero’‘-block matrix 8. We save the storage space. On the other hand, for any Gabor triple (g, a, b), we have presented a general approach to compute the GDGA g based on the O(N) CG-method and an algorithm to calculate the tight Gabor atom based on the binomial series. We have shown that g is the closest biorthogonal window function in the sense of Qian
S. Qiu/ Signal Processing
49 (1996)
167-186
dual Gabor atom (i)
Gabor atom
tight Gabor atom (i)
O.1N -IN 200
200
400
0.1 0.05 0 LLI
400
200
400
x ,o-3 dual atom (iv)
x , o-3 dual atom (iii)
dual Gabor atom (ii)
183
dl
0.02 0.01 0
: 4
2
z
0
1: L&J
-0.01 LLI 200
‘W 400
200
400
200
tight Gabor atom (iii)
tight Gabor atom (ii) 0.1 []
400
tight Gabor atom (iv)
0.06
0.04
0.05 0.02
0.02 0.040 ;.h.,
0 LJLA 200
0 LAJ 200
400
400
200
Fig. 4. A Gabor atom of the exponential shape and the associated dual and tight Gabor atoms associated pairs: (i) (12, 12), (ii) (10, lo), (iii) (8, 8) and (iv) (8, 6). The signal length is 480.
400 with four different lattice constants
Gabor atom
loo
I
I
I
200
300
400
500
600
700
,
I
800
900
GDGA
x 1o-4
1 1000
1100
GTGA 0.047 0.03 0.02 0.01 0 LL
200
400
Fig. 5. A Gabor atom and the associated
600
800
1000
200
400
600
800
1000
GDGA and GTGA. The signal length N = 1200 and the lattice constants
(a,!~) = (25, 16)
S. Qiu/ Signal Processing
184
49 (1996)
(I). The signal from the Gabor space 101
(I). The reconstructed
i
J
-10’ 200
400
600
800
167-186
1000
signal
101
i
I
-1OL
200
(II). The original Chirp
400
600
800
(II). The best approximate
1000 Chirp
I
200
400
600
800
1000
1
200
400
600
800
Fig. 6. (I) A signal from the Gabor space and the reconstructed signal with GDGA S. (II) An original Chirp signal from the Gabor space. The Gabor atom 9 and the GDGA 6 are shown in Fig. 5. Original
Fig. 7. An original shown in Fig. 8.
image (Lena)
Image
(128 by 128)
of size 128x 128, and the orthogonal
and Chen [ 191. The algorithms presented in this paper have the following important applications: (1) to check whether a Gabor triple (g, a, b) generates a Gabor frame or not; to predetermine the quality of a Gabor triple (2) (g, a, b), i.e. the stability of Gabor synthesis with
(3)
(9, a, b); to obtain X E CN.
the best approximation
of a signal
Best Approximate
projection
1000 Chirp and the best approximate
Image
onto the Gabor space where the 2-D Gabor atom is
Acknowledgements The author acknowledges fruitful discussions with Prof. Hans G. Feichtinger. The author would like to thank the hospitality of the Department of Mathematits of the University of Vienna where the author began to prepare the paper. Special thanks to Prof. K. Grochenig for carefully reading the paper and making good suggestions. The author would like to thank
S. Qiul Signal Processing 49 (1996) 167-186 2D Gabor atom
185
2D generalized x 10-3 3
0.8
0.6
0.2
0 12 120
Fig. 8. A 2-D Gabor atom and the associated
GDGA.
2D Gabor atom (contour plot)
I 50
100
150
200
Original image (463 by 331)
Fig. 9. A 2-D Gabor atom, an original
250
300
350
400
Best Approximate
450 Image
image and the best approximate
image
dual Gabor
186
S. Qiu / Signal Processing 49 (1996)
Mr. K. Marinelli for setting up the computer enabling to do numerical experiments.
facility
References [l] L. Auslander, I. Gertner and R. Tolimieri, “The discrete Zak transform application to time-frequency analysis and synthesis of nonstationary signals”, IEEE Trans. Signal Process., Vol. 39, 1991, pp. 825-835. [2] M.J. Bastiaans, “Gabor’s signal expansion of a signal into Gaussian elementary signals”, Proc. IEEE, Vol. 68, No. 4, 1980, pp. 538-539. [3] 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. l-35. [4] J.J. Benedetto and D.F. Walnut, “Gabor frames for L2 and related spaces”, in: J. Benedetto and M. Frazier, eds., WA VELETS Mathematics and Application, CRC Press, Boca Raton, FL, 1994, pp. 91-162. [5] J.J. Benedetto, “Frame decompositions, sampling, and uncertainty principle inequalities”, in: J. Benedetto and M. Frazier, eds., WAVELETS Mathematics and Application, CRC Press, Boca Raton, FL, 1994, pp. 247-304. [6] I. Daubechies, Ten Lectures on Wavelets, CBMS-NSF Regional Conf Series in Applied Mathematics, Vol. 61, Philadelphia. SIAM, 1992. [7] I. Daubechies, “The wavelet transform, time-frequency localization, and signal analysis”, IEEE Trans. Inform. Theory, Vol. 36, 1990, pp. 961-1005. [8] R.J. Duffine and A.C. Schaeffer, “A class of nonharmonic Fourier series”, Trans. Amer. Math. Sot., Vol. 72, 1952, pp. 341-366. [9] T. Ebrahimi and M. Kunt, “Image compression by Gabor expansion”, Optical Engrg., Vol. 30, No. 7, 1991, pp. 873880. [lo] H.G. Feichtinger and K. Grdchenig, “Theory and practice of irregular sampling”, in: J. Benedetto and M. Frazier, eds., WAVELETS Mathematics and Application, CRC Press, Boca Raton, FL, 1994, pp. 305-363. [l l] G.H. Golub and CF. Van Loan, Matrix Computations, The Johns Hopkins University Press, Baltimore, 1993. [12] C.E. Heil and D.F. Walnut, “Continuous and discrete wavelet transforms”, SIAM Rev., Vol. 31, No. 4, 1989, pp. 628-666. [13] E. Isaacson and H.B. Keller, Analysis for Numerical Methoo!s, Wiley, New York, 1966.
167-186
[14] A.J.E.M. Janssen, “The Zak transform: A signal transform for sampled time-continuous signals”, Philips J. Res., Vol. 43, 1988, pp. 2369. [15] A.J.E.M. Janssen, “Duality and biorthogonality for WeylHeisenberg frames”, J. Fourier Anal. Appl., Vol. 1, No. 3, 1995, pp. 403436. [16] P. Lancaster and M. Tismenetsky, The Theory of Matrices, Academic Press, New York, 1985. [17] S. Li and D.M. Healy Jr, “On pseudo frame decompositions and a class of discrete Gabor expansions”, IEEE Trans. Signal Process., Submitted. [18] R.S. On; “The order of computation for finite discrete Gabor transform”, IEEE Trans. Signal Process., Vol. 41, No. 2, 1993, pp. 122-130. [19] S. Qian and D. Chen, “Discrete Gabor Transform”, IEEE Trans. Signal Process., Vol. 41, No. 7, 1993, pp. 2429-2438. [20] S. Qian and D. Chen, “Optimal biorthogonal analysis window functions for discrete Gabor transform”, IEEE Trans. Signal Process., Vol. 42, No. 7, 1994, pp. 694697. [21] S. Qiu and H.G. Feichtinger, “The structure of the Gabor matrix and efficient numerical algorithms for discrete Gabor expansions”, in: K. Katsaggelos, ed., Visual Communications and Image Processing ‘94, SPZE-Proc., Vol. 2308, 1994, pp. 11461157. [22] S. Qiu and H.G. Feichtinger, “Gabor-type matrices and discrete huge Gabor transforms”, IEEE Znternat. Conf Acoust. Speech Signal Process.-95, Vol. 2, Detroit, MI, May 1995, pp. 1089-1092. [23] S. Qiu and H.G. Feichtinger, “Discrete Gabor structure and optimal representation”, IEEE Trans. Signal Process., Vol. 43, No. 10, 1995, pp. 2258-2268. [24] S. Qiu, “On the best approximation by undersampled Gabor families”, IEEE Trans. Signal Process., Preprint, May 1994. [25] D. Stewart, L. Potter and SC. Ahalt, “Computationally attractive real Gabor transform”, IEEE Trans. Signal Process., Vol. 43, No. 1, 1995, pp. 77-84. [26] G. Strang, Linear Algebra and its Applications, Harcourt Brace Jovanovich Publishers, San Diego, 1986. [27] R. Tolimieri and R.S. Orr, Poisson summation, the ambiguity function and the theory of Weyl-Heisenberg frames, J. Fourier Anal. App., Vol. 1, No. 3, 1995, pp. 233-248. [28] I. Wexler and S. Raz, “Discrete Gabor expansion”, Signal Processing, Vol. 21, 1990, pp. 207-220. [29] J. Yao, “Complete Gabor transformation for signal representation”, IEEE Trans. Image Process., Vol. 2, No. 2, 1993, pp. 152-159. [30] M. Zibulski and Y.Y. Zeevi, “Oversampling in the Gabor scheme”, IEEE Trans. Signal Process., Vol. 41, No. 8, 1993, pp. 2670-2687.
PROCESSING Signal Processing
49 (1996)
167-186
Generalized dual Gabor atoms and best approximations by Gabor family Sigang Qiu* Department of Mathematics, Uniuersity of Connecticut U-9, 196 Auditorium Road. Storrs, CT 06269-3009, USA Received
13 January
1995; revised
14 December
1995
Abstract Let g be a Gabor window of length N and (a,b) be a pair of lattice constants. u and b are considered as time and frequency gaps for a TF-lattice with N*/ab elements. The corresponding time-frequency shifts of g form a so-called Gabor family. In this paper, we investigate the structural properties of the discrete Gabor transforms. We address the problem of finding the best approximation of a signal x E CN by linear combinations of a Gabor family. We consider critical sampling, oversampling and undersampling, and do not assume that the Gabor family is a flame. For the task we determine the (generalized) dual Gabor atom (GDGA). This amounts to determining the pseudoinverse of the Gabor matrix and can be solved by the conjugate-gradient (CG) algorithm with O(N) complexity for fixed lattice constants (a, b). We provide an easy practical criterion for checking whether a Gabor triple (g, a, b) generates a Gabor frame or not. We propose an efficient algorithm for estimating the Gabor frame bounds and an algorithm for determining tight Gabor atoms. Zusammenfassung Es seien g ein Gaborfenster der Lange N und (a, b) ein Parr von Gitterkonstanten; a und b werden als Zeit- und Frequenzhicken fir ein Zeit-Frequenz-Gitter mit N’/ub Elementen betrachtet. Die entsprechenden Zeit-Frequenz-Verschiebungen von g bilden eine sogenannte Gabor-Familie. In diesem Beitrag diskntieren wir die strukturellen Eigenschaften der diskreten Gabortransformationen. Wir sprechen das Problem an, wie die beste Approximation eines Signals x E CN durch lineare ijberlagerungen einer Gaborfamilie zu finden ist. Dabei betrachten wir kritische Abtastung wie Unter- und ijberlastung und nehmen nicht an, dalj die Gabor-Familie einen Rahmen bildet. Zu diesem Zweck bestimmen wir das (verallgemeinerte) duale Gaboratom (GDGA). Das ftihrt dazu, dag die Pseudoinverse der Gabormatrix bestimmt werden mug, was fiir feste Gitterkonstanten (a, b) mit dem “Conjugate-Gradient” (CC) Algorithmus mit einer Komplexitit O(N) zu l&en ist. Wir stellen ein einfaches praktisches Kriterium bereit, mit dem zu testen ist, ob ein Gabortripel (g, a, b) einen Gaborrahmen erzeugt oder nicht. Wir schlagen einen effizienten Algorithmus zur Schatzung der Gaborrahmen-Grenzen vor sowie ein Verfahren znr Bestimmung dichter Gaboratome.
Soit g une fen&e de Gabor de longueur N et (a, 6) une paire de constantes de treillis. a et b sont consider& comme etant les intervalles de temps et de frequence pour un treillis TF de N’/ub elements. Les deplacements temps-Mquence correspondants de g forment ce qu’on appelle une famille de Gabor. Dans cet article, nous discutons des proprietes structurelles des transfotmees de Gabor disc&es. Nous considerons le probleme de trouver la meilleure approximation d’un signal x E CN par * Tel.: INT++(860)-486-1280;
fax: INT++(860)-486-4238;
e-mail: [email protected].
0165-1684/96/$15.00 @ 1996 Elsevier Science B.V. All rights reserved PIISOl65-1684(96)00015-l
168
S. Qiu / Signal Processing
49 (1996)
167-186
des combinaisons lineaires d’une famille de Gabor. Nous etudions l’echantillomrage critique, le SW et le sous-tchantillonnage, et ne supposons pas que la famille de Gabor est une trame. Pour cette t&he, nous determinons l’atome de Gabor dual (generalise) (GDGA). Ceci revient a determiner le pseudo-inverse de la matrice de Gabor et peut etre resolu en utilisant l’algorithme du gradient conjugue (GC) de complexite O(N) pour des constantes de treillis fix&es (a, b). Nous donnons un critere pratique simple pour verifier si un triplet de Gabor (g, a, b) genere une trame de Gabor ou pas. Nous proposons un
algorithme efficace pour estimer les limites de la trame de Gabor et un algorithme pour determiner les atomes de Gabor. Keywords: Gabor expansion; (Generalized) dual Gabor atom (GDGA); 9%symmetry
1. Introduction The theory of Gabor expansions has been thoroughly studied by many scientists, mostly via the Zak transform, see e.g. [3] and the references cited there. Since Gabor elementary functions are nonorthogonal, the main difficulty of doing Gabor analysis is how to obtain suitable Gabor coefficients efficiently. As we know, these coefficients can be obtained by short-time Fourier transforms (STFT) with a so-called biorthogonal Gabor analysis window. It is the desire of finding efficient determinations of this dual Gabor window (atom) that has attracted considerable attention from both mathematicians and engineers. Actually, the oversampled and the critically sampled Gabor families which are frames have been studied essentially by the following two methods: - Via the Zak transform and the frame method: We refer to [l-3,7,14,18,27,30]. - Via matrix decompositions: We refer to [9,17,1925,28,29]. The main purpose of this paper is to consider discrete Gabor transforms generally. Specifically for an arbitrary Gabor triple (g,a,b), we present a general method to obtain the associated (generalized) dual Gabor atom (GDGA) g based on the CG(conjugategradient)-method, which can be implemented with O(N) complexity for fixed lattice constants (a, b). For a signal X, using the calculated & we are able to determine the best approximation of x by linear combinations of the Gabor family. We propose an algorithm to determine the tight Gabor atom gt based on the binomial series. In Section 2, we list some notation which have been widely used in our previous papers [21,23] and some preliminary results about discrete Gabor expansions. Section 3, the main part of this paper, will discuss the Rotated(g)-symmetric structure of the Gabor
matrix in detail and will give an easy way of computing the Gabor frame bounds. Section 3 also contains the corresponding numerical algorithms. In Section 4 we provide some numerical results. We draw our conclusions in Section 5.
2. Notation and preliminaries We use essentially the notation introduced in [21]. For the paper to be self-contained, we recall the main definitions. We use i= J--f to denote the imaginary unit and j, k, 1, a, b are always nonnegative integers. is considered as an N(i) A signal x=(x(k)}::: periodic row vector in the complex space CN. We use either Xk or x(k) to denote the (k + I)thcoordinateofasignalx(k=O,l,...,N-1). The inner product of two signals x,y E CN is
(ii)
&Y> = CL;’ xo’M.0 Signals are-treated as row vectors. For a matrix A E CNxN, it induces a linear mapping (operator) A : x H x * A. The adjoint operator corresponds to x c+ x * A’, where A’ is the conjugate transpose (or Hermitian transpose) of A. We use AT to denote the transpose of A. The time translation operator Tk corresponds to the cyclic rotation:
Tk(x041, .a. %-I) 3
:= (X-k,Xl__k, ... . x-l,&
.. . .
XN-k-1).
The modulation operator (or the frequency translation operator) is given by M
(x0,x1,
:=
*..,
Q-1)
. XN--Id--I) (X0,X~W,X~W2, ..)
>
S. Qiu/Signal
Processing 49 (1996)
where w=e (iii)
-2xi’fN for 1 E R.
(vi)
The following rotation operator (introduced in [21,23]) with the rotation number a acting on vectors is also used for the cyclic rotation:
167-186
169
satisfy that ub
r&(x,
T,g
a) := (XN_-a,XNma+l, ... , XN-13x0,
If y is a column
. . . . XN-a-1).
T(i- I jag Mb!? MbTaY
vector, we define rot@, a) :=
[rotdy’, a)]‘. The matrix rotation on matrix A with a rotation number a is obtained by rotating all the column vectors of A. Precisely if A E CPxq and A =
‘-W,
a,b) =
(2) MbT,,- I )a.4
(AI, AZ, ... , Aq), then
MmbTnag
rotm(A, a) := (rot(A,, a), rot(A2, a), ... . rot(A,, a)). For any number k = 1,2, ... . we define inductively rotmk(A, a) = rotm (rot&-‘(A, a),~). It
is
easy
to
check
that
rotm’(A, a) =
rotm(A, ku). (iv)
(v)
n=O m=O
where the Gabor elementary functions gnm are given as g,,,,, = M&,Tnag for n = 0,1, . . . ,6 - 1 and m = O,l,..., b - 1. The collection {gnm}nfl is called a Gabor family. {c,,}~~ are usually referred to as the Gabor coefficients. The Gabor operator S associated to (g, a, b) is defined as
=
d-l
6-l
C
Cb,
n=O
m=O
Then it is easy to see the following sentation for the Gabor operator: Sx=x*S
For a signal g E CN, we call (g,u, b) a Gabor triple if both a and b are positive integers and divisors of N. We call g a Gabor atom. The constants a and b are called time and frequency gaps. They form a pair of lattice constants. We call G = N/u and b”= N/b the dual lattice constants. A Gabor expansion of a signal x with respect to a Gabor triple (g,u, b) is a series expansion of the form
Sx
MC&_, y,T(d-I jag
snm)snm,
for
x E P.
(1)
Unless the contrary is explicitly stated, we assume from now on that lattice constants (a, b)
matrix repre-
forxECN,
where S = [GAH(g, a, b)]’ * [GAE%(g,a, b)] is called the Gabor matrix. In general, for any two Gabor triples (gi,u, b) and (gz,u, b), we call the matrix product [GAB(gi, a, b)]’ * [GAEi(gz, a, b)] the Gabor matrix associated to (gi,u, b) and (92, a, b). For a Gabor triple (g, a, b), the associated Gabor operator S and the Gabor matrix S may not be invertible. However, pinv(S) is uniquely defined, where pinv(S) denotes the pseudoinverse [26] (or the Moore-Penrose inverse [ 16, pp. 432-4401) matrix of S. Then we define the associated generalized dual Gabor atom (GDGA) as g := g * pinv(S). Since the Gabor matrix S is positive semidefinite [l 1, 16,261, the square root S1’2 is well defined [16, pp. 180-1841. We define the generalized tight Gabor atom (GTGA) as gt = g * pinv(S1j2). If a Gabor triple (g, a, b) generates a Gabor frame, the associated Gabor operator S is the usual Gabor frame operator (cf. [12, pp. 635-6361 or [6, pp. 58591). In this case, pinv(S) = inv(S) is the usual inverse of S. The associated GDGA g and GTGA gt are the dual and tight Gabor atoms in the usual sense (cf.
170
S. Qiu/ Signal Processing 49 (1996) 167-186
[27], [6] or [12]). In this case, we call g and gt the usual dual and tight Gabor atoms. The following proposition is used to obtain the best approximation of a signal x E CN by linear combinations of a Gabor family {gnm}n,nz. This is equivalent to determine the orthogonal projection of x onto the Gabor space: Span{g,,}, where Span{g,,} is the linear space spanned by the Gabor family {gnm}nfl [16, pp. 83-871. Proposition 1. Let J be the GDGA associated to (g, a, b). For a signal x E CN, the orthogonal projection of x onto Span{g,,} can be represented as
This yields that g,, E Span{g,,}. Then using (5), we derive that MmbTna(g * pinv(G)) * {G * pinv(G))
gnnm =
= MmbTna{g * pinv(G) * G} * pinv(G) = (MmbT,,g) * pinv( G) = gnm* pinv( G). This proves Eq. (6). Hence, i-1
6-1
7
x(X,
d-1 Km)gnm
6-1
= c
n=O m=O
x(X9
gnm * pinv(G))gnm
n=O m=O d-l
=
6-l
C X(X * pinv(G),gnm)gnm n=O m=O
In particular, ifx E Span{g,,}, 8-l
= {x * pinv( G)} * G.
n=O m=O
n=O m=O
6-l
8-l
(3)
then
We conclude from Eqs. (5) and (6) that if x E SPan{gnm1
6-1 X =
{X * pinv(G)}
* G
n=O m=O
n=O m=O
=
d-l
6-1
71
Fl(X,
n=O
m=O
inm>gnm,
and Let us sketch the proof of Proposition 1. For further details, we refer to [24]. Let G = [GAB(g,a, b)]’ * [GAEl(g,a, b)] be the Gabor matrix associated to (g, a, b). It is clear that G and thus pinv(G) are Hermitian symmetric. For a signal x E Span{g,,}, we see that x is in the row space of the Gabor matrix G (cf. [16, pp. 78-801, [26, pp. 137-1401). Hence, we have (see [16, pp. 428-4381, [26, pp. 140-1411) {x * G} * pinv(G) = {x * pinv(G)} * G =x
for X E Span{g,,}.
(5)
x = {x * G} * pinv( G)} 8-l = C
z-1 = T;7
=
(MmbTnag)
i.e. &,, = gym.
gnm)
(6)
In fact, by the commutative property of the Gabor operator with both T, and Mb, it is easy to see that
7:
tx3 gnm)!7nm.
Eq. (4) is proved. For the general case (3), since P(x) E Span{g,,}, using (4), we see that d-l = 71
6-1 x(P(X)9
= M,bT,,,
(g * @NW) (g
*
I?-1 6-1
Similarly,
pinv(G*) * G)
=MmbTna (g * pinv(G’))
* G.
gnm)tjnm
=
d-l
6-l
7:
FItxy
gnm)Lm
n=O m=O
= {x * G} * pinv(G) = {x * pinv(G)} * G
n=O m=O CTnnm= MmbTnn
* PWG))
6-l
n=O m=O
* PWG),
(gnm
n=O m=O
P(X) (g *WV(G))
x(X,
n=O m=O
In order to prove Eqs. (3) and (4), we need to show first that MmbTna
5-l
i-1
5-l
n=O m=O
S. Qiu / Signal Processing 49 (1996)
Eq. (3) is obtained. The sketch of the proof of Proposition 1 is complete.
167-186
171
More precisely
Remark. Since Eq. (6) holds without the assumption
if B is given (by Theorem A) as
Cl,1
0
..
0
c2,2
’
0 0
..
that {gnm I,,, is a frame, Eq. (6) is an extended version of the result presented in [4, Theorem 3.1 I], [6, pp.
0
831, [7, PP. 9711.
0 0
c6+l,1
0
Ck2.2
... ...
. cCZ,Cl 0 0
...
.. B=
3. Mathematical setup
0
As we have seen in [21,23], the Gabor matrix S associated to Gabor triples (gi,a,b) and (gz,a,b) is completely determined by its first N x a block matrix B. All the other a”- 1 block matrices are obtained from B by proper rotations. We state these descriptions in the following theorem. and Feichtinger [21,23]). Let (gl,a, b) and (gz,a, b) be two Gabor triples, S = [GAB(g,, a, b)]’ * [GAE!(gz, a, b)] be the associated Gabor matrix. Assume that B is thejirst N x a block matrix qf S. Then S is of the following block matrix ,form:
Theorem A (Qiu
C(b-l)6+l,l
0
0
forr=a”-1.
(7)
For k = 1,2, ... , a” - 1, rotn#(B, a) is exactly the (k + I)-th block matrix of size N x a. Moreover, the possibly nonzero entries of B are distributed only in the I-th diagonals for 1 = 0,416, f26,. . . , f(b - 1)b”. We call B the first block corresponding to Gabor triples (gi, a, b) and (92, a, b). If gi = g2 = g, we call B the first block matrix associated to (g, a, b). Extracting all the zero-entries of B, we obtain a compressed “nonzero”-block matrix B of size b x a. The Ith row-vector is exactly the Ith diagonal of B for 1 = 0, &, ,.. , (b - 1)J. We call B the (first) “nonzero”-block matrix corresponding to Gabor triples (gl,a, b) and (gz,a, 6). When gi = g2 = g, we call B the “nonzero”-block matrix associated to (g, a, b). By Theorem A and the definition of rotm, every (sub-)diagonal is a-periodic. We consider B as u-periodic in column-subscripts and b-periodic in row-subscripts.
0
:::
0
. ..
0
1)&+2,2
0
then the b x a “nonzero”-block
B=
CI,I
c2,2
ClT+1,1 .
%+2,2
. ctY+n,a
...
‘(b-
I$++a,a
matrix B is
...
C,,ll
...
C6+a,o
..
:
t'(b-I)&,,,
a), ... , rotm’(B, a)]
...
..
i S = [B, robn(B,
‘(b-
0
: ‘(b-1)6+2,2
“.
1.
C(b-I)i++a,ai
Since B is of small size b x a, B can be easily built up. We give a sufficient condition with r? for a Gabor triple (g, a, b) to generate a Gabor frame.
Corollary 1 (Sufficient condition). Let B = (Ui,j)bXa be the “nonzero”-block matrix associated to a Gabor triple (g, a, b), and assume that for any 1
then (g, a, b) generates a Gabor frame.
Proof. Condition
(8) implies that the Gabor matrix S associated to (g, a, b) is strictly diagonal dominant. Therefore, S is nonsingular and (g, a, b) does generate a frame. 0 In the following, we show that if gi = g2 = g then the Gabor matrix associated to (g, a, b) has the rotated-symmetry (%!-symmetry) structure.
S.
172
Qiu/ Signal Processing 49 (1996) 167-186
Theorem 1 (W-symmetry structure). Assume that gr = g2 = g, and let
In general, for Y = 46 + 1, s = (b - q)b” + 1 with q = 1, 2, ...) bo - 1, we can show similarly that
S = [B, rotm(B), ... , rotmr(B)],
di = d;&
wherer=a”-
1,
(9)
be the Gabor matrix associated to (g, a, b). Writing b,, =
i(b+l)
ifbisodd,
i(b + 2)
if b is even,
then the Gabor matrix S has at most abo independent nonzero elements. They are from the upper part of the N x a block matrix B. Proof. Using Theorem A, we show first that for r = qb”+l,s=(b-q)b”+lwithq=1,2, ....ba-1.the sth subdiagonal elements of B are completely determined by the rth subdiagonal elements of B via the rotation and the conjugation. Conversely, we can determine the rth subdiagonal from the sth subdiagonal. Without loss of generality, we take Y = 6 + 1 and s = (b - 1)b” + 1. Noticing that the diagonal elements of the block matrix are of periodicity N = bb”,for k = 1, 2, ... , a, the general elements of the sth subdiagonal can be written as
4 = S((bl)K+k),k i-1
=bCT,,g((bn=O
l)h+k-
l)Tnag(k-
1)
6-l
=bxT,,g(k-b-
l)T,,g(k-
1)
l)Tnag(k-
1)
n=O
for k = 1, 2, ... , a.
This means that the sth subdiagonal of B is compketely obtained via a rotation with the rotation number qb (or the remainder of qb divided by a) and the conjugation of the rth subdiagonal of B.
Hence, all the independent entries of S are located in the rth subdiagonal of B, for Y = qb + 1 with q = 0,l , ..*, bo - 1. The total number of independent entries are at most abo. The proof is complete. 0 Remark. Since abo c iab, the B-symmetry of the Gabor matrix has reduced the computational cost by half. In general, if gr # 92, the Gabor matrix corresponding to (gk, a, b) for k = 1,2 does not have the Wsymmetry. Using the introduced “nonzero”-block matrix, we rewrite Theorem 1 as the following equivalent form. Theorem 2. Let B be the “nonzero”-block matrix associated to (g,a, b). Assume that 8j is the jth row vector of B for j = 2, ... , bo. Then the (b - j + 2)-th row vector 8b_j+2 of 8 can be written as Bb-j+2 = rot(Bj, mod((j - l)b, a)), where mod(Cj - 1)b”, a) is the remainder of Cj - 1)b” divided by a. In particular if the signal length N is divisible by ab (especially for the critical cases), then
d-l
=bxT,,g(k-6n=O
bb_j+2 = 6
for j = 2, ... , bo.
If the Gabor atom is a real-valued signal, then
= d;_6
fib--j+2 = rot(LTj,mod(~ - l$,
Therefore, di = d;_6
for k = 1, 2, ... , a.
s = (b - 1)b” f 1 -th subdia( > gonal elements of B can be obtained via a suitable rotation and the conjugation from the (r = 6 + 1)-th subdiagonal elements of B. This yields that the
a))
for j = 2, ... , bo. Theorem 3. Let S be the Gabor matrix corresponding to a Gabor triple (g,a, b). Then for any k = O,l, ... . Sk is W-symmetric. More generally, if p(z) is any real polynomial, then p(S) also has the W-symmetry.
S. Qiu / Signal Processing 49 (1996)
For the proof of Theorem 3, we need the following interesting proposition.
Proposition 2. Let S k/2 be the linear operator realized by the matrix Ski2 for k = 0,1,2,. . . . Then Ski2 fork = 0,1,2,... commute with both the time translation operator T, and the frequency modulation operator Mb.
167-186
173
Proof of Theorem 3. Set h = S(k-1)/2g, then the Gabor operator Sh associated to (h, a, b) is defined as 6-l
&x =
d-l
c
x(x,h,,)h,,
m=O
n=O
X E 6?
On the other hand, by Proposition Sk1 &k-i)/=
2,
(S (S(k-‘Y=X))
Proof. If k = 0, the fact is trivial. We only need to consider the cases for k > 1. For k = 1, S is the associated Gabor operator. It is known that S commutes with both T, and Mb. Since S is positive semidefinite, Ski2 and Ski2 are well defined. Assume that U is the unitary matrix diagonalizing S, then there is a diagonal matrix D = diag(ai,az,...,aN) such that
for
6-I z-1
x ~(S(k-1)‘2x,Ynm~Ynm
= SW 1Y=
m=On=O 6-l
=
d-l
x
~(+,s(k-‘)~=gn,)S(k-‘):=g,m
m=O
n=O
6-l d-l
S=
U’*D*U,
(10)
where c113 1x23 . . . UN 2 0, diag (ai, ~12,.. . , UN) denotes the diagonal matrix with aj (j = 1,2,. . . , N) being its diagonal entries. From Eq. (lo), it is easy to see that
=
c
~(~,hmJhm.
m=O
n=O
Thus,
This implies that = (J’
e
diag
Sk = [GAJ3(h, a, b)]’ * [GAB(h, a, b)].
u:‘~,ui12,.. . , OLD* U.
(11)
>
By the polynomial interpolation theory [13], there is an interpolation polynomial L(x) such that L(uj) = $I’,
The l-norm
for j = 1,2,.. .,N.
Ski2 = U’* diag(L(ai),L(a&...,L(a,v))
[ll,
13, 161 of a matrix A = (a;,j),lx,
is defined as
From Eq. ( 11 ), we can easily check that
=U’*L(D)*U=L(S).
Therefore, Sk is &!-symmetric by Theorem 1. For any real polynomial p(z), p(S) also has the .C&symmetry. 0
n
*U
0
Then by the fact that S commutes with both T, and Mb, we deduce that Ski2 commutes with both the timetranslation operator T, and the modulation operator Mb. 0 With the same argument given above, we obtain the following corollary.
Corollary 2. Zf we define ,Ykj2 = (pinv(S))k’2 for k = 1,2,..., then Sk12for k = 0, fl, 4~2,. . . commute with both T, and Mb.
IPIll = ,y2mCladiI. i=l
It is easy to compute, especially for a small-size matrix. We show that the l-norm of the small “nonzero”-block matrix can be used to estimate the Gabor frame upper bounds. Suppose (g, a, b) generates a Gabor frame, S is the associated Gabor frame operator. Let 8 be the dual Gabor atom, 8 and B> be the b x a “nonzero”-block matrices associated to (g, a, b) and (6, a, b), respectively. Then
Theorem 4. (Gabor frame bounds).
S. Qiu/ Signal Processing 49 (1996)
174
Proof. Assume that S and S are the Gabor matrices associated to (g, a, b) and (g, II,b), respectively, then IlSll =
IPII
and IIS-‘II = IF’II
= IISII-
Since S = [B, rotm(B, a), ... . rotm’(B, a)] and B is obtained from B by filling in some proper zero-entries, it is clear that
Ilrota% a)111= IlBlll = 11~111. Hence,
ll~lll = IPIll = IPIll. Similarly,
II% = II&Ill,
i.e.,llS-‘Il~ = Il&ll~.
Because the norm (spectral norm) of S is majorized (see [13; 16, pp. 358-3621) by the l-norm of S, we see that
A
= ~~IIsIl~llsII1 W’ll1
= 11411.
Therefore,
167-186
(2) If u = 1, let 8, be the “nonzero”-block corresponding to (F(g), b,l), then the best Gubor upper bound is obtained us (14~ Ilo3/N, Proof. Since the Gabor matrix associated to (g,u, 1) is a nonnegative diagonal matrix [23, Corollary 81, the first result is obvious, while the second result is deduced from [23, Theorem 31. 0 Algorithm 1 (Nonzero-block and Gabor matrix). Let (g,u, b) be a Gabor triple and S be the associated Gabor matrix. Then we can construct S with only z(& ) M ;Nb multiplications and via proper rotations and the conjugation by the following steps: Step 1. We calculate the first bo row-vectors of the b x a “nonzero”-block matrix l?. The total computation needs only (ubo )C = Nbo M +Nb multiplications. The (k, j)-th general entry of the “nonzero’‘-block B can be computed as d-l ~~~=gCTnog((k-1)5+j-l)TnagV-l) n=O
for j = 1, 2, ... , a and k = 1, 2, ... , bo. Step 2. We obtain the remaining b - bo row-vectors in the lower part of B. If Bi denotes the jth u-elements row vector of B for j = 1,2, ... , b, then for j = 27 ... 3 bo,
Remark. (1) Let &,,, and ;Iminbe the maximal and minimal eigenvalues of the Gabor matrix S, then Aand &in are the best upper and lower Gabor frame bounds. Tables 1 and 2 show that the approximate Gabor frame bounds ub := l[~llt andLb := l/l(B>ll, calculated from Theorem 4 approximate well to I,, and kin. (2) With the simple and efficient formula presented for the Gabor frame bounds, the condition number of the Gabor matrix can be estimated as
which can be used to check the quality (that is, the stability of the associated Gabor expansion and synthesis) of a Gabor signal. Corollary 3. (1) If b = 1, then the “nonzero”-block corresponding to (g, a, 1) is a vector B = (tj)yE, E Cu. llB]l, =maxy=, tj is the best Gubor upper bound.
Bb-j+Z = rot(Bj,mod((‘j - l)b,
u)).
Step 3. Filling in some extra zero-entries into 8, we get the first N x a block matrix B of the Gabor matrix. Step 4. By the rotation operator rotm with the rotation number a, we can successively obtain the other 2 - 1 block matrices of the same size N x a. The Gabor matrix S is built up.
Algorithm 2 (Gabor matrix-vector multiplication). Let (g,a, b) be a Gabor triple, S and S be the associated Gabor operator and Gabor matrix, B be the associated “nonzero”-block matrix of size b x a and x be an arbitrary signal. Then the general elements of y = Sn = x * S := (Yo, Yl, ... . Yv_1) can be computed from B = (ck,l),,, by
S. Qiu/ Signal Processing 49 (1996)
167-186
175
Table 1 Comparison between the Gabor frame bounds ub (Lb,) and A,, (in,,“). We use ub (Lb) to denote the approximate Gabor upper (lower) bound obtained by Theorem 4. A,, (A,,) stands for the best Gabor frame upper (lower) bound. The Gabor atom shown in Fig. 3 is of a typical Gaussian shape Lattice
(ib
(a. b)
( uap*rw.
(12, 12) (IO, 10) (8, 8) (8, 6)
597.50260416615 717.00312961905 896.26049300375 897.82568699835
(0, b)
KDZ =
(12, 12)
570.215206083184 66.364320832716 11.420117443573 2.910703099535
(10, 10) (8, 8) (8, 6)
*mm 1
ub/Lb
Gxact
597.50260416597 717.00312961882 896.26049300366 897.82568699835
K =
I .047853832802 10.804045315651 78.480847279568 308.456636179028
1
1.047853832801 10.8040453 15657 78.4808447279583 308.456636179024 ub
-
Lb
ub
+
Lb
knmax/&b
570.215602083367 66.364320832654 11.420117443569 2.910703099535
0.91961152311793 0.78133544866770 0.5433 1293495694 0.26092308033431
0.99649869507642 0.9703 1069303042 0.83897092687839 0.48858301203230
Table 2 Comparison between the Gabor frame bounds ub (Lb) and A,,,= (&,). (ib (Lb) and &,, same meanings as those in Table 1. The exponential Gabor atom is shown in Fig. 4 Lattice
ub
ta,b)
( uapprox.
1max 1
(10, 10) (8, 8) (8, 6)
1922.76264640444 2306.32832577053 2882.42656567521 2886.37655929048
(a, b)
Km =
(12,
(12,
12)
12)
(10, 10) (8, 8) (8, 6)
Ub/Lb
813.739718040677 375.122874805591 168.773569670759 47.896900894645
(&act
Lb )
1919.88421461691 2303.70222154611 2880.82779694440 2886.37655929050
K =
~Lpprox
(&in)
have the
bin 1
2.36287181733500 6.14819431357204 17.07866090228586 60.26228222238036
vexact
)
2.37681424145194 7.11090203642794 17.08182182956302 60.26228222237835 Lb
a-v%
‘h
-
a+cb
ub
+Lb
&xxx/bin
807.755263803913 323.967649919045 168.648744009180 47.896900894647
wherej=(r-l)a+s--lforr=1,2,...,Zand s = 1, 2, , a. The derivation of this algorithm is similar to the proof of Algorithm 2 in [23]. With the precalculated “nonzero’‘-block 8, we only need Nb multiplications to perform Algorithm 2. We should also notice that if S is the Gabor matrix corresponding to (gk, a, b) for k = 1,2, i.e., S = [GAB(gi, a, b)]’ * [GAEl(gz,a, b)], then Algorithm 2 does also hold. In practical applications, we do not need to build the N x N Gabor matrix S. We only need to work with the small b x a “nonzero”block matrix B by Algorithm 2.
0.93226337332680 0.90180718516839 0.85705390609887 0.74749893572402
0.99754522830333 0.99468258876562 0.9882 1960329939 0.95909761225340
We now present the CG-method for obtaining the GDGA y” associated to a Gabor triple (g,a, b). The standard CG-method is stated in the following [ 11, p. 5231 Algorithm A (Conjugate gradient [ 111). If A is positive definite, then we have the following iterative algorithm for computing x so that x * A = y, for a vector y: Starting with x0 = 0, Q = y, p1 = t-0the iterative step is as follows, for k > 1, _ Set/-$=Oand/?k=*; _ Set pi = ?i, and pk = yk-_l + Pkpk-I;
S. Qiu/ Signal Processing
176
- setuk=&; - set - set
xk
=
xk__l
•,-
‘tkpk;
rk
=
rk-1
-
akpk
*A.
It is known that N is the upper bound of the number of iteration steps and XN is the exact solution to the linear system x * A = y if the algorithm is performed in infinite precision. In this case, we say that the CGalgorithm is convergent after N iterations. For a Gabor triple (g, a, b), let S be the associated Gabor matrix to (g,a,b). We show in the following theorem that the CG-method can always be applied and be convergent after at most ab iterations to obtain the associated GDGA j by solving the linear equation g * S = g. Using the Gabor matrix-vector multiplication algorithm 2, the product Pk *S for each CG-cycle can be implemented in only O(N) operations.
49 (1996)
167-186
at most ub-dimensional subspaces of CNxN and CN, respectively. If g is real-valued, by Theorem 1 and Theorem 3, for any 1 = 0,1, ... . the (r = (b - q)6 + 1)-th subdiagonal of B(l) can be obtained by (s = q6 + 1 )-th subdiagonal of B(‘) via only the rotation with the rotation number qb”.This yields that the dimension of 9 and consequently the dimension of Pg over R are at most ubo. 0 Proof of Algorithm 3 . Step 1. We show that the CG-procedure works and converges after at most ub iterative steps. Following the arguments in [ 11, pp. 5 16-5231, we only need to show that Ukis always defined. Equivalently, we need to show the following: ifrk-_l #O,
thenpk*S*pL#O.
(13)
Algorithm 3 (Generalized dual Gabor atom (GDGA)). Let S be the Gabor matrix corresponding to a Gabor triple (g, a, b). If we replace both A and y by S and g in Algorithm A, respectively, then the CG-algorithm is convergent after s (s < ub) iterations and the resulting x, = 6 is the solution for g * S = g. Moreover, g = x, is exactly the GDGA corresponding to (g, a, b), that is, i = g * pinv(S). In particular if g is real-valued, the algorithm is convergent after at most ubs iterations.
In fact, suppose that there exist pk and r&t such that Pk * S * pi = 0 but rk_1 # 0, then Pk * [GAB(g, a, b)]’ * [GAB(g, a, b)] * p; = 0, which yields that Pk * [GAB(g,u, b)]’ = 0. Since S = [GAB(g,u, b)]’ * [GAB(g, a, b)], g is in the row space of GAB(g, a, b). Inductively, it is easy to check that pk is in the row space of GAB(g,u, b). We deduce that (see [ 11, pp. 50]), pk = 0. This implies that
In order to derive Algorithm 3, we need the following lemma.
and
But rk_l~i_~ = 0 holds (see [ll, Theorem 10.2.2]), r,&_l = 0 follows. So we get a contradiction and Eq. (13) is true. Therefore, if rk_1 # 0, then we always have that pk * S * pk # 0. By [l 1, Theorem 10.2.2 and Theorem 10.2.31, after k - 1 steps of the CG-algorithm, the residuals
P’x = Span{x * s’ : I = 0, 1, ... }, for x E CN.
ra,rl, ... . rk__l E Span{g,g*S,
Then both 9’ and 9x (especially, ~9~) are at most ubdimensional subspaces of CN XNand CN, respectively. If g is a real-valued Gubor atom, then the dimensions of both 9’ und Ps over R are at most ubo.
are mutually orthogonal. By Lemma 1, the dimension of the subspace Span{g * S’, for 1 = 0,1, ... , } c CN is at most ub. Hence, there exists s (s Gab) such that the (s + 1 )-th residual r, = 0, which yields that x, * S = g. Write # = x,, then g * S = g. If g is a real-valued signal, the algorithm is convergent after ubo steps by Lemma 1, Step 2. We show that the solution g to the linear system 6 * S = g based on the CG-method is exactly the GDGA y = g * pinv(S).
Lemma 1. Write B = Span{& : I= O,l, ...}
Proof. For I = 0, 1, ... , let B(‘) be the first N x a block matrix of S’. By Theorem 3, S’ are all completely determined by B (‘1 which have at most ub nonzero entries. Because the nonzero entries of B(l) for 1 = 0, 1, ... ) are in the same positions, both 9 and 9’~ are
F-k-1 + bkpk--l
= 0.
... . g*Sk-t}
S. Qiu/ Signal Processing 49 (1996) 167-186
By the CG-iteration steps, we see that there exists a (real) polynomial p(t) with at most ab degree such that $=g*p(S)
and
g*S=g.
We show that p(S) = pinv(S) and therefore 6 = y * pinv(S). Since S and p(S) are Hermitian and commute with each other, with the definition of the pseudoinverse [ 11, 16,261, we only need to show that s * p(S)
* s = s.
Let S be the Gabor Since S commutes and the frequency 0, 1, ... , N - 1, we
(14) operator associated with (g, a, b). with both the time translation T, modulation A&, for any k, 1 = have
Yk.1* P(S) * S = S (P(S)&l)
= S MM,,,
177
bound for the number of CG-iteration steps is ab. Therefore, the total operations for the computation of the GDGA i is about 2ab(bN f4N). If the lattice constants (a, b) are fixed, then we obtain an O(N) CG-algorithm. - Jacobi-method: When the Gabor atom is given as a typical Gaussian function (for example, the Gabor atom shown in Fig. 1, the associated Gabor matrix is (strictly) diagonal dominant. We use the Jacobimethod [ 11, 131 instead to improve the implementation. In this case, the computations and decompositions of Gabor matrices can also be realized by the associated “nonzero”-block matrices. Fig. 1 shows the comparison of convergent rates between the Jacobi-method and the CG-method. - PCG-method: In practice, for a given Gabor triple (g, a, b), the Gabor atom g is well localized. Even if the associated Gabor matrix S = does not have good diagonal dominance, we can use the preconditioned conjugate gradient (PCG) algorithm [ 111 with the preconditioner (sj,k)N
= (S(P(S))& = (@ * S),,
= (9 * P(S) * S),,, gk,l.
=
Thus,
P = diag(sl,l,S*,2,...,SNsV)
gk,l * p(S) * S = gk,l
for
k, 1 = 0, 1, ... . N - 1. (15)
By the definition of the Gabor Eq. (2)), Eq. (15) yields that GAB(g,a,b)
xN
basic
* p(S) * S = GAB(g,a,b).
matrix
(see
(16)
Since S = [GAB(g, a, b)]’ * [GAB(g, a, b)], Eq. (14) follows immediately if we multiply both sides of Eq. (16) by [GAIQ,a, b)l’. The proof is complete.
0
Remark. - Complexity analysis: Algorithm
3 is always applied with Algorithm 2. By Algorithm 2, we do not need to build up the full N x N Gabor matrix S, we only need to compute and save the associated b x a “nonzero”-block matrix l?. This largely saves the storage space. Using Algorithm 2, the calculation of pk * S for every CG-iteration requires only Nb multiplications and Nb additions. Since L? can be precalculated with only bo(a + N) multiplications and boN addtions in total, Algorithm 3 can be performed fast. By Algorithm 3, the upper
to improve Algorithm 3. Equivalently, if we set SO=S*P-’ and go=g*P-‘, then we need to solve the equation S * SO = go by the CGalgorithm. Since SO = S * P-’ = (sj,k/sk,k)N:
Corollary 4 (Gabor frame condition). Let i be the associated GDGA to a Gabor triple (g,a, b), i& = (ukj)hXa be the “nonzero”-block matrix associated
S. Qiu / Signal Processing 49 (1996)
178
167-186
Comparison of Convergent Rates of Jacobi- and CG-methods
Fig. 1. Convergent-rate comparison of the Jacobi- and the CG-methods; a Gabor atom and the associated dual Gabor atom. The Gabot atom is well concentrated and the associated Gabor matrix is diagonal dominant. The signal length is 128 and lattice constants are (8, 8).
with the Gabor triples (g,a, b) and (@,a, b). Then (g, a, b) generates a Gabor frame if and only if
‘k,j = “5, =
1 0
ifk=l, ifk # 1,
(17)
for k = 1,..., a and j = l,..., bo. Proof. Let Ss = [GAT3(g,a, b)]’ * [GAB(&a,b)]. Then (g, a, b) generates a Gabor frame if and only if SO is an identity matrix. The condition (17) follows. Now suppose that condition (17) is satisfied, we show that (g,a, b) generates a Gabor frame. Since there exists a real polynomial p(t) such that g = g * pinv(S) = g * p(S), following the similar argument to the proof of Theorem 3, we can verify that the Gabor matrix So can be identified with p(S) * S. By Theorem 3, SO satisfies the g-symmetry and the condition (17) implies that So is an identity matrix. Therefore, (g,a, b) generates a Gabor frame. Cl Remark. For an arbitrary Gabor triple (g,a, b), we can always check whether (g, a, b) generates a Gabor
frame or not via two steps: ( 1) Determine g by Algorithm 3. (2) Test only abo M $ab number of conditions (17). In the following, we show that the GDGA g obtained by Algorithm 3 is the closest one to the Gabor atom g in the sense of Qian and Chen [ 191. Theorem 5 can be seen as an extended discrete version of Proposition 3.2 in the recent work of Janssen [15]. Theorem 5. Assume that g and S are the GDGA and the Gabor matrix associated to a Gabor triple (g,a,b), then 5 = g * pinv(S) is the unique solution of the problem
minimize (I&
- &/Iover
with y * [GAB(g,b,a”)]’ = $e. Proof. Write &? = Span(g,,),
y E CN
(18)
then {gnm}is a frame of 4?. Let S and Sd be the Gabor operators associated with {gnm}in CN and in 4, respectively. Then S& is invertible. It is easy to check that S”fr = SI”4r
S. Qiu/ Signal Processing 49 (1996)
and
179
and
S~x=x*s
forxEd.
$
S,i’ has the following
matrix representation:
Si’x
for x E A.
= x * pinv(S)
gt := g * (pinv(S))‘j2
= g * pinv(S’!2)
is exactly the GTGA corresponding to (g, a, b). For any 1’ satisfying the biorthogonality y * [GAB(g, 6, a”)]’ = $e, we can express gt as d-l 5-l St = c C(gt, n=O m=O
Ynm)Snm.
On the other hand, gt E &? can also be expressed as
n=O
m=O
Since (g,a,b) generates a frame of .&’ and i is the dual Gabor atom associated to (g, a, b) in &‘, the frame property (cf. [6, Proposition 3.2.4; 12, Proposition 2.1.51 or originally [8, Lemma VIII]) yields that C-1
d-l 6-l =
IG II=0m=O
cc
d-l
6-l
=
6-1
= cc
Ikh
n=Om=O
d-1
&m)l’
6-l
Ih Ym)12= y x IbY mm)12
G :x n=O
G&l)12
m=O
n=O
= (59 9)
~llill lbll GIIYII Ilsll.
It follows that g = g * pinv(S) to the problem (18). 0
Hence, 6 = g * pinv(S) = S>‘g is the usual dual Gabor atom associated to (g,a, b) in A. Similar to Janssen’s argument [ 15, p. 171, we can show that @ is the unique solution to the problem (18). In fact, setting gt = S,>“2 g, then g, E _&’and
Ils”ll
167-186
m=O
11;‘.4Yl12~/lYl12>
where y.4 is the orthogonal projection of y onto A. Since y * [GAB(g, 6, Z)]’ = $e implies that (y, g) =
is the unique solution
Remark. By choosing a different initial vector ~0, we can modify Algorithm 3 to obtain desirable biorthogonal Gabor analysis atoms. In fact, let h be a signal, if we intend to obtain a biorthogonal Gabor atom i. which is close to h, then we set x0 = g - h * S and apply Algorithm 3. Assume go is the output, then g = go + h is what we want. Algorithm 4 (Tight Gabor atom). Given a Gabor triple (g, a, 6) which generates a frame, the following procedure is used to compute the associated tight Gabor atom g,: Step 1. By Algorithm 1, we obtain the “nonzero”block matrix B of size b x a and the Gabor frame upper bound ub by Theorem 4. Step 2. We calculate the dual Gabor atom and henceforth the lower Gabor frame bound Lb. then & := YE is associated to Step 3. Set r = &, the Gabor frame operator S, := KS which satisfies that
&~llS~lI
=rllSll<2Uh Lb
+
uh
This yields that
:= p< 1.
III -Srll<~
Step 4. Since /I1 - &(I 67
< 1,
00 s-112 r - c &(I - 0, k=O
where h(h-l)...(h-k+l)
!Xo= 1,
& =
for h = -i Hence,
and k = 1, 2,
abIN,
1 .2...k .
30
s--1/2
llnhd=2-2Rei%% =2(1- Nlhll llsll) ab
=
J;cuk(I
-
Sr)k.
k=O
Therefore, gt = s-g.
the tight Gabor atom can be obtained as
S. Qiu/ Signal Processing 49 (1996)
180
Iteratively, set St(O)= fig,
gt@+l) = J;a,+,
(gp
4. Numerical results
for n 20,
- S,gJ”‘) .
Then we have lim gt@) = gt. n-+oo The error estimation after n iterations is II&)
-
St11
In fact,
II2 J;
Ilgt(n)-&II=
ak(I-Sr)kg
k=n+l
=
&(I-$)“+’ (1
=
/I
gak+n+l k=O
167-186
(I-Sr)kg (1
E/,gll.
Remark. Even if (g,a,b) does not generate a Gabor frame, Algorithm 4 also works. Algorithm 5 (Best approximation). Assume that (g,a,b) is a Gabor triple and i is the associated GDGA, we obtain the best approximation of a signal x by linear combinations of the Gabor family via the following steps: Step 1. We obtain the “nonzero”-block of size b x a corresponding to the Gabor triples (g,a,b) and (&a, b). Step 2. Applying Algorithm 2 and Proposition 1, we obtain xapproX. In particular, if (g,a, b) generates a frame or a signal x is in the Gabor space spun{gnm}, then we can reconstruct a signal x via the above two steps. Remark. In the two-dimensional case, if the Gabor atom is given as a 2-D separable one, then we can calculate the orthogonal projection of an image onto the 2-D Gabor space by using Proposition 1 and the higher dimensional results presented in [23].
This section is devoted to numerical simulations. The performances of the previously proposed algorithms are illustrated. Given an arbitrary Gabor triple, the numerical constructions of the associated GDGA, GTGA and the best approximation of a given signal are shown. All the numerical experiments were carried out using MATLAB 4.0 on a Sun4 Station. As we have demonstrated, the algorithms proposed in this paper are in a very general sense. We do not assume that a given Gabor family is a frame. We do not restrict that the signal length is the power of 2. Signals can be either real valued or complex valued. Gabor atoms (Gabor window functions) can be chosen arbitrarily. Since Gabor atoms with Gaussian types are typical (they are minimizing the uncertainty principle [5]), most of our examples are experimented with Gabor atoms of Gaussian shapes. When a Gabor atom is well concentrated, the associated Gabor matrix is strongly diagonal dominant. We use the Jacobi-method instead of the CG-algorithm to speed up the convergence. In Fig. 1, the Gabor atom is of a Gaussian shape. The signal length is 128 and the lattice constants (a, b) are (8,8). Since Gabor window functions are usually with good localizations, the associated Gabor matrices have some diagonal dominance (even if they are not strictly diagonal dominant). In this case, the simple PCGmethod is employed to improve the ordinary CGalgorithm. Fig. 2 compares the convergent rates of the PCG, CG and the ordinary frame-method. Both PCG and CG algorithms work much faster, while PCG is better than CG. Meanwhile, if (g, a, b) does not generate a frame, the frame method does not work. However, the CG-method can always be applied to obtain the GDGA. Tables 1 and 2 illustrate the comparison between the best Gabor frame bounds and the estimated ones calculated according to Theorem 4. For Table 1, we take a Gabor atom with the Gaussian shape shown in Fig. 3. For Table 2, we choose a Gabor atom with the exponential shape shown in Fig. 4. Lattice constant pairs are (12,12), (10, lo), (8,8) and (8,6) for both Tables 1 and 2. It is clear that the calculated Gabor bounds ub and Lb approximate well to the best Gabor
S. Qiu/ Signal Processing 49 (19%)
167-186
181
frame-method
\ lo-l5 -
‘, CG
I
1o-2o -_ 10” Gabor atom
dual Gabor atom
I
r
I :
1. 0.8. 0.6. 0.4. 0.2. o-
-50
0
50
Fig. 2. Convergent-rate comparison of the PCG, CG and the frame methods. The solid line denotes the convergent rates for the PCG-method; the dash-dotted line is for the CC&method; the dotted line is for the frame method. Gabor atom is of the Gaussian type. Signal length is 120 and lattice constants are (8, 6).
bounds A,,, and &in. From the Gabor upper and lower bounds, we can predetermine the stability of the Gabor syntheses with the associated dual Gabor atoms shown in Figs. 3 and 4. In fact, the data in the column with Ub/Lb = K, [ 1 l] estimate the condition number of the associated Gabor matrices. If U&b is small, then the corresponding Gabor synthesis is stable. The last two columns of Tables 1 and 2 show clearly why the CG-method converges much better than the frame-method [6,7]. In fact, Theorem 10.2.5 in Golub [ 1 l] shows that the CC convergence rate is majorized by
g’.=
VG-6 Jrr;;+&=
Jlcoo-1 &+1’
while the frame-method &,-Lb .--==. ‘=‘Ub+&,
convergence
rate depends on
&.-1 &c + 1
The data in the last two columns of Tables 1 and 2 show that 01 is much smaller than a2. Therefore, the CG-method converges much faster.
Fig. 3 illustrates a Gabor atom and the associated dual and tight Gabor atoms. The signal length N is 480 and the lattice constants are (i) (12,12), (ii) (IO, lo), (iii) (8,8) and (iv) (8,6). For a random signal, the reconstruction relative errors with the dual and tight Gabor atoms are of order 10-15. Fig. 4 shows the exponential Gabor atom and the dual (tight) Gabor atoms corresponding to lattice constants: (i) (12, 12), (ii) (10, lo), (iii) (8, 8) and (iv) (8, 6). The signal length is 480. For a random signal, the reconstruction relative errors with the dual (tight) Gabor atoms are of order 10-14. Next we show some experiments with Gabor triples that do not generate Gabor frames. Given a Gabor triple (g,a,b) and a signal x with the same signal length as g. If x is in the Gabor space Span(g,,), then we get the reconstructed signal of x by Algorithm 5. If x is not in the Gabor space Span{g,,}, then we can only obtain the best approximation xapprox E Span{g,,} of the signal x. Fig. 5 illustrates the Gabor atom and the associated GDGA @and GTGA qt. The signal length N is 1200.
S. Qiu/ Signal Processing 49 (1996)
182
167-186
tight Gabor atom (i)
dual Gabor atom (i)
Gabor atom
I 1 0.5 0 FYIL
J
200
200
400
,
0 -5 L!!!!N
20
3 2 1 0 -1 LLI
5
200
0.06 0.04 0.02 0 0.08 •LI
10 0 200
400
tight Gabor atom (ii)
x IoA dual atom (iv)
x ,03 dual atom (iii)
x to-3 dual atom (ii)
0.02
•LI 200
400
tight Gabor atom (iii) 0.06 0.04
400
tight Gabor atom (iv) 0.04
I
I
0.02
0 200
400
0w 200
400
400
Fig. 3. A Gabor atom of the Gaussian shape and the associated dual (tight) Gabor atoms. The signal length N is 480 and the dual and tight atoms associated with different lattice constants pairs: (i) (12, 12); (ii) (10, 10); (iii) (8, 8) and (iv) (8, 6).
The Gabor atom with respect to the lattice constants (a, b) = (25, 16) does not generate a frame. In Fig. 6 we show (I) the original signal from the Gabor space Span{g,,} and the reconstructed signal, the (relative) reconstruction error is of order 10-14; (II) the original Chirp signal and the best approximate Chirp by the linear combinations of the Gabor family. Fig. 7 shows the original and the best approximate image (Lena) of size 128 x 128 with linear combinations of the 2-D Gabor family. A 2-D Gabor atom and the associated GDGA are shown in Fig. 8. The 2-D lattice constants (a, b) are taken to be a = (16,8) and b = (8,4). Fig. 9 gives another numerical simulation for the 2-D case. The original image given as a rectangular matrix of size 463 x 331 and the best approximate image by the linear combination of the 2-D Gabor family are illustrated. The 2-D Gabor atom of size 480 x 480 with the 2-D lattice constants (a,b) for
a = (30,24) and b = (15,10) Gabor frame.
does not generate a 2-D
5. Conclusion We have developed the W-symmetry feature of Gabor matrices. The B-symmetric property reduces the computational complexity by approximately half, for building up the Gabor matrix, its first N x a block matrix and the associated b x a “nonzero”-block matrix. The algorithms presented in [2 1,231 are improved. In practice, we do not need to construct the full N x N Gabor matrix S, the Gabor matrix-vector multiplications can be realized completely by the small b x a (ab d N) “nonzero’‘-block matrix 8. We save the storage space. On the other hand, for any Gabor triple (g, a, b), we have presented a general approach to compute the GDGA g based on the O(N) CG-method and an algorithm to calculate the tight Gabor atom based on the binomial series. We have shown that g is the closest biorthogonal window function in the sense of Qian
S. Qiu/ Signal Processing
49 (1996)
167-186
dual Gabor atom (i)
Gabor atom
tight Gabor atom (i)
O.1N -IN 200
200
400
0.1 0.05 0 LLI
400
200
400
x ,o-3 dual atom (iv)
x , o-3 dual atom (iii)
dual Gabor atom (ii)
183
dl
0.02 0.01 0
: 4
2
z
0
1: L&J
-0.01 LLI 200
‘W 400
200
400
200
tight Gabor atom (iii)
tight Gabor atom (ii) 0.1 []
400
tight Gabor atom (iv)
0.06
0.04
0.05 0.02
0.02 0.040 ;.h.,
0 LJLA 200
0 LAJ 200
400
400
200
Fig. 4. A Gabor atom of the exponential shape and the associated dual and tight Gabor atoms associated pairs: (i) (12, 12), (ii) (10, lo), (iii) (8, 8) and (iv) (8, 6). The signal length is 480.
400 with four different lattice constants
Gabor atom
loo
I
I
I
200
300
400
500
600
700
,
I
800
900
GDGA
x 1o-4
1 1000
1100
GTGA 0.047 0.03 0.02 0.01 0 LL
200
400
Fig. 5. A Gabor atom and the associated
600
800
1000
200
400
600
800
1000
GDGA and GTGA. The signal length N = 1200 and the lattice constants
(a,!~) = (25, 16)
S. Qiu/ Signal Processing
184
49 (1996)
(I). The signal from the Gabor space 101
(I). The reconstructed
i
J
-10’ 200
400
600
800
167-186
1000
signal
101
i
I
-1OL
200
(II). The original Chirp
400
600
800
(II). The best approximate
1000 Chirp
I
200
400
600
800
1000
1
200
400
600
800
Fig. 6. (I) A signal from the Gabor space and the reconstructed signal with GDGA S. (II) An original Chirp signal from the Gabor space. The Gabor atom 9 and the GDGA 6 are shown in Fig. 5. Original
Fig. 7. An original shown in Fig. 8.
image (Lena)
Image
(128 by 128)
of size 128x 128, and the orthogonal
and Chen [ 191. The algorithms presented in this paper have the following important applications: (1) to check whether a Gabor triple (g, a, b) generates a Gabor frame or not; to predetermine the quality of a Gabor triple (2) (g, a, b), i.e. the stability of Gabor synthesis with
(3)
(9, a, b); to obtain X E CN.
the best approximation
of a signal
Best Approximate
projection
1000 Chirp and the best approximate
Image
onto the Gabor space where the 2-D Gabor atom is
Acknowledgements The author acknowledges fruitful discussions with Prof. Hans G. Feichtinger. The author would like to thank the hospitality of the Department of Mathematits of the University of Vienna where the author began to prepare the paper. Special thanks to Prof. K. Grochenig for carefully reading the paper and making good suggestions. The author would like to thank
S. Qiul Signal Processing 49 (1996) 167-186 2D Gabor atom
185
2D generalized x 10-3 3
0.8
0.6
0.2
0 12 120
Fig. 8. A 2-D Gabor atom and the associated
GDGA.
2D Gabor atom (contour plot)
I 50
100
150
200
Original image (463 by 331)
Fig. 9. A 2-D Gabor atom, an original
250
300
350
400
Best Approximate
450 Image
image and the best approximate
image
dual Gabor
186
S. Qiu / Signal Processing 49 (1996)
Mr. K. Marinelli for setting up the computer enabling to do numerical experiments.
facility
References [l] L. Auslander, I. Gertner and R. Tolimieri, “The discrete Zak transform application to time-frequency analysis and synthesis of nonstationary signals”, IEEE Trans. Signal Process., Vol. 39, 1991, pp. 825-835. [2] M.J. Bastiaans, “Gabor’s signal expansion of a signal into Gaussian elementary signals”, Proc. IEEE, Vol. 68, No. 4, 1980, pp. 538-539. [3] 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. l-35. [4] J.J. Benedetto and D.F. Walnut, “Gabor frames for L2 and related spaces”, in: J. Benedetto and M. Frazier, eds., WA VELETS Mathematics and Application, CRC Press, Boca Raton, FL, 1994, pp. 91-162. [5] J.J. Benedetto, “Frame decompositions, sampling, and uncertainty principle inequalities”, in: J. Benedetto and M. Frazier, eds., WAVELETS Mathematics and Application, CRC Press, Boca Raton, FL, 1994, pp. 247-304. [6] I. Daubechies, Ten Lectures on Wavelets, CBMS-NSF Regional Conf Series in Applied Mathematics, Vol. 61, Philadelphia. SIAM, 1992. [7] I. Daubechies, “The wavelet transform, time-frequency localization, and signal analysis”, IEEE Trans. Inform. Theory, Vol. 36, 1990, pp. 961-1005. [8] R.J. Duffine and A.C. Schaeffer, “A class of nonharmonic Fourier series”, Trans. Amer. Math. Sot., Vol. 72, 1952, pp. 341-366. [9] T. Ebrahimi and M. Kunt, “Image compression by Gabor expansion”, Optical Engrg., Vol. 30, No. 7, 1991, pp. 873880. [lo] H.G. Feichtinger and K. Grdchenig, “Theory and practice of irregular sampling”, in: J. Benedetto and M. Frazier, eds., WAVELETS Mathematics and Application, CRC Press, Boca Raton, FL, 1994, pp. 305-363. [l l] G.H. Golub and CF. Van Loan, Matrix Computations, The Johns Hopkins University Press, Baltimore, 1993. [12] C.E. Heil and D.F. Walnut, “Continuous and discrete wavelet transforms”, SIAM Rev., Vol. 31, No. 4, 1989, pp. 628-666. [13] E. Isaacson and H.B. Keller, Analysis for Numerical Methoo!s, Wiley, New York, 1966.
167-186
[14] A.J.E.M. Janssen, “The Zak transform: A signal transform for sampled time-continuous signals”, Philips J. Res., Vol. 43, 1988, pp. 2369. [15] A.J.E.M. Janssen, “Duality and biorthogonality for WeylHeisenberg frames”, J. Fourier Anal. Appl., Vol. 1, No. 3, 1995, pp. 403436. [16] P. Lancaster and M. Tismenetsky, The Theory of Matrices, Academic Press, New York, 1985. [17] S. Li and D.M. Healy Jr, “On pseudo frame decompositions and a class of discrete Gabor expansions”, IEEE Trans. Signal Process., Submitted. [18] R.S. On; “The order of computation for finite discrete Gabor transform”, IEEE Trans. Signal Process., Vol. 41, No. 2, 1993, pp. 122-130. [19] S. Qian and D. Chen, “Discrete Gabor Transform”, IEEE Trans. Signal Process., Vol. 41, No. 7, 1993, pp. 2429-2438. [20] S. Qian and D. Chen, “Optimal biorthogonal analysis window functions for discrete Gabor transform”, IEEE Trans. Signal Process., Vol. 42, No. 7, 1994, pp. 694697. [21] S. Qiu and H.G. Feichtinger, “The structure of the Gabor matrix and efficient numerical algorithms for discrete Gabor expansions”, in: K. Katsaggelos, ed., Visual Communications and Image Processing ‘94, SPZE-Proc., Vol. 2308, 1994, pp. 11461157. [22] S. Qiu and H.G. Feichtinger, “Gabor-type matrices and discrete huge Gabor transforms”, IEEE Znternat. Conf Acoust. Speech Signal Process.-95, Vol. 2, Detroit, MI, May 1995, pp. 1089-1092. [23] S. Qiu and H.G. Feichtinger, “Discrete Gabor structure and optimal representation”, IEEE Trans. Signal Process., Vol. 43, No. 10, 1995, pp. 2258-2268. [24] S. Qiu, “On the best approximation by undersampled Gabor families”, IEEE Trans. Signal Process., Preprint, May 1994. [25] D. Stewart, L. Potter and SC. Ahalt, “Computationally attractive real Gabor transform”, IEEE Trans. Signal Process., Vol. 43, No. 1, 1995, pp. 77-84. [26] G. Strang, Linear Algebra and its Applications, Harcourt Brace Jovanovich Publishers, San Diego, 1986. [27] R. Tolimieri and R.S. Orr, Poisson summation, the ambiguity function and the theory of Weyl-Heisenberg frames, J. Fourier Anal. App., Vol. 1, No. 3, 1995, pp. 233-248. [28] I. Wexler and S. Raz, “Discrete Gabor expansion”, Signal Processing, Vol. 21, 1990, pp. 207-220. [29] J. Yao, “Complete Gabor transformation for signal representation”, IEEE Trans. Image Process., Vol. 2, No. 2, 1993, pp. 152-159. [30] M. Zibulski and Y.Y. Zeevi, “Oversampling in the Gabor scheme”, IEEE Trans. Signal Process., Vol. 41, No. 8, 1993, pp. 2670-2687.