- Email: [email protected]
Parallel lattice structure of block time-recursive generalized Gabor transforms
SIGNAL
PROCESSING ELSEVIER
Signal Processing 57 (1997) 195-203
Parallel lattice structure of block time-recursive generalized Gabor transforms’ Chao LIP, Sanjay Joshib, Joel M. Morrisb9* “Department of Computer & Information Sciences, Towson State University, Baltimore, MD 21204, USA bComputer Science and Electrical Engineering Department, University of Maryland Baltimore County, Baltimore MD 21250, USA
Received 19 April 1996; revised 3 October 1996
Abstract The problems of efficient computation of generalized Gabor transforms are considered in this work. Block timerecursive Gabor transforms are defined for both complex and real data sequences, fast algorithms are derived through two different approaches, one by the finite Zak transform, the other by biorthogonal functions. All the block timerecursive algorithms for the computation of generalized Gabor transforms can be unified by a parallel lattice structure. The unified lattice structure computes the transformed data sequentially with data sliding into/out-of a block of size N at a time. Gabor transform computation is based on the most recent L-point data, where L = MN. This model accommodates timely processing of the sequentially received data in digital signal transmission. Implementation on parallel computer system will be discussed. 0 1997 Elsevier Science B.V. Zusammenfassung In diesem Beitrag wird das Problem der etiienten Berechnung der verallgemeinerten Gabor Transformation aufgegriffen. Blockorientierte zeitrekursive Gabor Transformationen werden fiir komplexe und reelle Datenfolgen definiert. Zwei verschiedene Ansltze fiir schnelle Algorithmen werden abgeleitet, von denen der eine auf der Zak Transformation und der andere auf biorthogonalen Funktionen basiert. All diese blockorientierten zeitrekursiven Algorithmen zur Berechnung der verallgemeinerten Gabor Transformation lassen sich auf eine einheitliche parallele Lattice (Kreuzglied-)Struktur zuriickfiihren. Diese parallele Lattice Struktur berechnet die transformierten Daten sequentiell, indem jeweils BlGcke von N Daten gleichzeitig verarbeitet werden. Die Berechnung der Gabor Transformation wird jeweils fiir die letzte aktuellste Datenfolge der LInge L durchgeftihrt, wobei L = MN. Dieses Model1 paljt die zeitaufwendige Verarbeitung sequentiell empfangener Daten an die Erfordernisse der digitalen Signaliibertragung an. Die Implementierung auf einem Parallelrechnersystem wird diskutiert. 0 1997 Elsevier Science B.V. Rbumi! Ce travail examine les probl6mes du calcul efficient des transform&es de Gabor gCnCralis6es. Les transformtes de Gabor par blocs rCcursives dans le temps sont dCfinies pour des stquences de don&es complexes aussi bien que rCelles; des algorithmes rapides sont ensuite d&iv&s $ partir de deux approches diffkrentes, l’une par la transform&e de Zak finie,
*Corresponding ’ This research
author. Tel.: 1410455 3503; fax: 1410455 3969, e-mail: [email protected]. was supported by an ONR Grant # NOOO14-89-51210.
0165-1684/97/$17.00 0 PIISO165-1684(96)00195-8
1997 Elsevier Science B.V. All rights reserved.
196
C. Lu et al. / Signal Processing 57 (1997) 195-203
l’autre par les fonctions bi-orthogonales. Tous les algorithmes de calcul par blocs et rkcursifs dans le temps de transform&s de Gabor ginCralisCespeuvent etre unifiCsdans une structure parallkle en treillis. Cette structure en treillis unifike calcule les donnires transformbes SCquenciellementg l’aide de donnies glissant dans/hors d’un bloc de taille N g chaque instant du temps. Ce calcul de la transform&e de Gabor est basi: sur les donntes de L points les plus ricentes, oti L = MN. Ce modkle s’adapte au traitement opportun de donnkes regues &quenciellement dans les transmissions numkriques des signaux. L’implantation sur un ordinateur parallble sera discutke. 0 1997 Elsevier Science B.V. Keywords: Gabor transform; Time-recursive algorithms; Zak transform; FFT; Lattice structure
1. Introduction
A theoretical foundation for designing algorithms for computing Gabor coefficients at critical sampling was established in [3,4]. The Zak transform played a central role in establishing clear and easily computable conditions for the existence of the Gabor expansion and for the stability of computation. In [19], the Zak transform provided the framework for computing biorthogonals for rationally oversampled Weyl-Heisenberg systemforming frames. A similar approach for integer oversampled Weyl-Heisenberg systems can be found in [lo]. Computational complexity was studied in [lS] for the finite discrete Gabor transform, and several algorithms were compared in the work. Fast algorithms for computing Gabor coefficients, at critical sampling, integer oversampling as well as rationally oversampling, were given in [2, 10, 13, IS] for discrete and periodic signal x(k) with periodicity L such that x(k) = x(k + 15). Implementation of these algorithms were carried out on RISC processors and parallel computing environments, and timing results were reported in those papers. The main computational task in those algorithms was a, or a set of two-dimensional FFTs. Real Gabor transforms were defined in [12, 181, and computationally attractive algorithms based on the discrete cosine transform were given in [18]. Time-recursive structures have been developed for computing the discrete Fourier transform (DFT) [14], discrete cosine (DCT)/sine (DST)/ Hartley (DHT) transforms [9], as well as the timevarying Fourier transform (TVFT) [l, 6,7]. All of the time-recursive algorithms proposed in these papers share a common characteristic: when the
moving-window moves forward one step, a new sample is included while excluding the oldest sample. In [l 11, we proposed a block time-recursive algorithm for the computation of the Gabor transform in the critical-sampling case, which is based on the most recent L-point input data where L = MN. Data is sliding in/out a block of size N at a time. The time-recursive algorithm presented in [ 1l] was derived from the finite Zak transform. The main computational task is a set of 1-D FFTs followed by an error correction stage to the previously computed Gabor coefficients. In this paper, block time-recursive Gabor transforms (BTRGTs) are defined for both complex and real data sequences, and fast algorithms are derived through two different approaches: one by the finite Zak transform, and the other by biorthogonal functions. All the block time-recursive algorithms for the computation of generalized Gabor transforms can be unified by a parallel lattice structure. The unified lattice structure computes the transformed data sequentially with data sliding into/out-of in blocks of size N at a time. The building blocks of the parallel lattice structures are similar for both the complex and the real Gabor transforms regardless of which approach is used. They consist of a set of 1-D FFTs or 1-D DCTs on the difference of the incoming block and outgoing block of data, and the transformed data is then added to the previously computed Gabor transform. These algorithms decouple the transformed data components, and interprocessor communication is needed only at the last addition stage with N-point data transmitted to the nearest neighbor. Hence, they are very attractive for parallel implementation for real-time signal processing.
197
C. Lu et al. / Signal Processing 57 (1997) 195-203
The rest of the paper is organized as follows. In Section 2, the block time-recursive complex Gabor transform is defined and, in Section 3, the finite Zak transform approach is used to derive a fast algorithm. In Section 4, the block time-recursive real Gabor transform is defined, and a fast algorithm based on the DCT is given. In Section 5, the block time-recursive Gabor transforms (BTRGTs) are unified into an unified parallel structure and, in Section 6, parallel implementation of the algorithm on the IBM SP2 will be discussed. Finally, in Section 7 the work will be summarized.
2. Block time-recursive complex Gabor transform For an integer L = MN > 0, the generalized Gabor expansion of a complex-valued, finite, discrete-time signal x(k) is M-l N-l x(k) = 1 1 a(m, n)g(k - mN)ejzXknlN, m=o n=O
&=
[x0,
xl,
X2N-I
. . . ,XN-l,XN,
,X2N,
...
... ,
XL-IIT
,
(6)
= cxo, XI, .” >;YM-IIT, F is a block-diagonal
matrix with N-point FFT matrices on its diagonal,
1
[F(N) F=
I
F(N)
(7)
..
f’(N)
and r is a block circulant matrix with each block being a diagonal matrix. The diagonal elements from each block in the first row of r comprise the biorthogonal function
(1)
where g(k) is a discrete-time periodic, window sequence with period L = MN. Under certain conditions [2] the expansion coefficients u(m, n) can be found by
(8)
1 l-1
l-2
.. .
with each r,,, an N x N diagonal matrix
1
ry(mN) a(m, n) = $ Lil x(k)y(k
- mN)e-jzAknlN.
(2)
k-0
rm=
y(mN+ 1)
The discrete-time periodic sequence y(k) is the biorthogonal function to the window function g(k) satisfying L-l ,zo
g(k + mN)y(k)e-j2”k”‘N
= i&,6,.
(3)
Eq. (2) is called the generalized complex Gabor transform. The matrix formulation of Eq. (2) is given in [ 181 by 4 = FM,
y(mN
(9)
An algorithm was provided in [18] for the calculation of the biorthogonal function. If we write Eq. (4) in detail, we have go = F(N)[T,x, f1
=w)~r,-,~,
+ i-lx1 +
G..
+
rM_l~Y_l],
+ TOT1+ ... +rM-2+11,
(4)
where
gM-1 =F(N)[r,$+r2X1
4 = [a(O, O),~(0, l), . . . ,a(O, N - l), u(l,O), . . . ,
2$f-llT,
+ ..*+roxM_l]. (10)
u((M- 1, N - l)]’ = Lao, t1, ...
+ N - 1)
(5)
We define the Gabor transform at block time t in order to derive a block time-recursive relation
C. Lu et al. 1 Signal Processing 57 (1997) 195-203
198
between the updated Gabor transform and the previous computed Gabor transform. This is consistent with the new block of data sliding in and the old block of data sliding out, the block size is N, and the transform size is L = MN. Thus,
Denote the FZT of Gaussian function g as G( i, 1). An algorithm for the computation of Gabor coefficients of critically sampled signals of size L is given as follows [lo]: Step 1. Compute the FTZ X(i, 1) of signal x. Case I. G never vanishes in IMM,N. Step 2. Compute
al(t) = fxW~M-I$
+ ~0$+1 + ...
(14)
+ rM-$*+M- 11,
(11) whereZM,,=((i,1):0di
+
roxt+,-II,
and at block time t + 1, the Gabor transform can be expressed as ao(r + 1) = al(t) + F(NVM-r(Xt+M
- x,)9
ar(r + 1) = a&) + F(JWM-,(x*+,
- X,)>
C
i=O
I=0
O
CM-l@
+
1)
=
$I@)
+
F(W~o($+,
-
-
xtx
-I
i, I)
=
G(i’‘)’
(15)
G,441, 6, I) E k 6, I) E IM,N.
0,
$1.
Eq. (12) is called the block time-recursive Gabor transform (BTRGT). The Gabor coefficients A(t + 1) at block time t + 1 are given by the pre;iously computed A(t) at block time t with an error correction term - F(N)ri(Xt+M - X,), i = 0, 1, 2, , . . , A4 - 1. The same formulation can be derived through the finite Zak transform (FZT).
3. Block time-recursive via the FZT
O
XG, 1)
p( F(w~l(xt+M
P(i, l)wLMmlwNni,
Case II. Let 5 be the zero set of G in IM,N. If the signal x has its FZT X vanishing on i, then Step 2. Compute
(12) aM-z(r + I) = aM-l(t) +
M-l
a(m, n) = C
Step 3. Compute P( i, 1): N-l
the 2D M x N inverse FFT of
M-l
a(m, n) = C
C
i=O
I=0
O
P(i, l)~i”~w;~‘,
O
The block time-recursive Gabor (BTRGT) can be derived as follows: The FZT of x at time t,
Gabor transform
(17)
transform
M-l
X(i, 1, t) =
Define
c x(i + mN + tN)w$‘. m=o
(18)
The FZT of x at time t + 1: M-l
X(i, 1) =
1
x(i + mN)wg,
X(i, 1,t + 1) = f
m=O
i,lEZ,
wM=e
x(i + mN + tN)w$‘“-”
m=l - jZn/M
.
(13)
X(i, 1) is called the jinite Zak transform (FZT) of x.
=
X(i, 1, t)wk”I + [x(i + MN + tN) - x(i + tN)] w;‘.
(19)
199
C. Lu et al. / Signal Processing 57 (1997) 195-203
Compute P(i, 1,t + 1) =
X(i, 1,t + 1) G(i, 0
X(L 1,r) =mWM
-1
+ [x(i + MN + tN) - x(i + tN)]
wiMI
(20)
2
G(i, 1)
I
then the Gabor transform of x at t + 1, a(m, n, t + 1)
I
Err(m, n, t)
a(m +
1, n, 1) + Err(m, n, t)
I
c
o(m, n, t + 1)
rxti + MN
T.2 F.2
+zi=o
A-‘
+ tm - x(i + tmi
G(i, ;)
I=0
Fig. 1. Implementation flow diagram of Gabor coefficient computation on infinite data.
~
‘-
(21)
Step 2. Compute N-point FFTs on E,(m, i, t): N-l
Set
Err(m, n, t) = c M-l
e&n,
i) =
c
2!!!-
lzo
G(i,
E,(m, i, t)wNin,
i=O
-I -h
w 1)
M
(22)
’
Odm
(25)
Step 3. Add the error term Err(m, n, t) to Gabor coefficients at time t:
Then a(m, It, t + 1)
a(m, n, t + 1) = a(m + 1, n, t) + Err(m, n, t), N-l
=
a(m
+
1, n, t)
+
c
[x(i
+
AffN
+
O
tw
(26)
i=O -
x(i + tN)] e,(m, i)wivi”
(23)
When we write Eq.(23) in detail, it has the same format as Eq. (12), the block time-recursive complex Gabor transform. Procedure to compute the BTRGT at block time t+ 1: 1. Pre-compute e&m, i) and store it in memory. 2. Compute the first set Gabor coejicients at t = 0: a(m, n, 0). 3. Time-recursively compute: Step I. Compute
The implementation flow chart of the time-recursive algorithm is given in Fig. 1.
4. Block time-recursive real Gabor transform For an integer L = MN > 0, the generalized real Gabor expansion of the real finite discrete-time signal x(k) is defined as
x cos((2k + l)nn/2N)
E,(m, i, t) = [x(i + MN + tN) - x(i + tN)]e,(m, i), -!-
O
andO
(24)
‘3
1
b(m, O)g(k - mN) . (27)
200
C. Lu et al. / Signal Processing 57 (1997) 195-203
The matrix formulation for the computation expansion coefficients b(m, n) is given by [ 183
of
If we write Eq. (28) in detail, we have + Jrlx, _ + ... + JfM_1~,+1],
LJ,,= C(N)[f,x,, B = C&,
(28)
where Jj = [b(O, O),b(O, l), . . . ,b(0,iv - l), 41, O),
b, = C(N)JIJfM_lxo
+ foxI
+ Lf-23-11, ‘.
>
b((M - 1, N - l)lT
bypI
= Cbo, b1, .‘. >~hf-IIT,
= C(N)J[J&
+ f2g
-t ...
(29) + b$-
x=
+ ...
(35)
11.
[x,~x,,...,xN-l,xN,....X2N-l. XZN, '.. ,XL-IIT
(30)
= C~O,~I, .'.,$f-llT.
The matrix C is block-diagonal with C(N) the N-point discrete cosine transform (DCT) matrices on its main diagonal (N being an even integer). Thus,
WV C(N)J
C=
L
C(N)J
and the biorthogonal Jf,
f, JpMel
PO
12
(31)
We now define the real Gabor transform at block time t in order to derive a block time-recursive relation between the updated Gabor transform and the previous computed Gabor transform. This is consistent with the new block of data sliding in and the old block of data sliding out, the block size is N, and the transform size is L = MN. Thus, I b,(r) = C(N)CQ, + Jr,xt+l + ...
+ Jhf-&+M-11, !v~l(t) = C(N)JIJf’M-lxc
_!
+ &,+l
+ ...
+ LZ~t+M-Il,
function matrix p is given as ... ...
JpM+ fMm2
&l(t)
= C(N)J[J&,
+ &xt+l
+ ...
^
(32)
+ roxt+&f-1 I,
(36)
and at block time t + 1, the real Gabor transform can be given as with each f, is an N x N diagonal matrix
r?(mN) ?(mN
%, =
+
bo(r + 1) = WVC(C(N)J)-‘~I@)
+ Jhe
1)
I(X~+M - ?,)I,
b,(t+ 1) = C(N)J[C+(N)&(t) f(mN + N - 1)
+ L2(3+h4
- g1,
(33)
and J is the N x N row exchange matrix
ro ...
0
b~-z(t + 1)=
1-l
C(N)C(C(W-
lb-
+ Jh(xt+, (34)
10) - x,)1,
bMml(t + 1) = C(N)JIC-l(N)b,(t) +
rotxt+hl - x,)1.
(37)
C. Lu et al. / Signal Processing 57 (1997) 195-203
It can be proved that l-1
C(N)JCl(N)
=
0
0 -1
...
0
...
0
1
.
Lo
0
...
(38)
- 1 -I
Set s = C(N)JC’(N);
(39)
then we have b,(t - 1) = ~~1,(t) + C(N)JLf-1(x,+,
- x,)3
b,(t + 1) = @b,(t) + C(N)JL,(:,+,
- x,), Fig. 2. Parallel lattice structure of block time-recursive Gabor transform.
&$.f-,(t + 1) =
&f-,(t) + wvJ~&,+,
t)M-l@ + 1) = &o(t) + C(W~o(~*+,
- xt), - xc). (40)
Eq. (40) is called the block time-recursive real Gabor transform. The Gabor coefficients lj(t + 1) at block time t + 1 are given by the previously computed B(t) at block time t with an error correction term C(N)JTi(x,+,$$ - x,), i = 1,2, . . . , M - 1. Eq. (40) has a similar format to the block complex Gabor transform in time-recursive Eq. (12).
5. Parallel lattice structure of block time-recursive Gabor transforms Both the block time-recursive complex (Eq. (12)) and real (Eq. (40)) Gabor transforms have the same lattice structure. Their implementation by a parallel lattice array consisting of M lattice modules is illustrated in Fig. 2. It should be noted that data processing at each stage is by an array size of N. The latest sliding-in data block x,+~ of size N subtracts the sliding-out data block xt, and then the difference is multiplied by the biorthogonal coefficients Ti at each module. An N-point FFT (or DCT for real) is calculated at each module followed by an addition of shifted previously computed Gabor coefficients. An array of size N needs to be transferred to the nearest-neighbor module in circular order.
To obtain a block of N transformed data, and update the other M - 1 blocks, each module of the lattice array carries N-point multiplications, N-point additions, and an N-point FFT (or DCT). The N-point FFT requires an order of N log N additions and an order of N log N multiplications for the Cooley-Tukey FFT algorithm, or an order of N multiplications for the Good-Thomas prime factor algorithm combined with the Winograd algorithm for each relatively prime factor. If we assume the Cooley-Tukey FFT algorithm for the N-point FFT computation, then the computational complexity is of the order of MN + MNlog N.
(41)
For the block time-recursive real Gabor transform, the parallel lattice structure in Fig. 2 will include a rotation of J and a trivial multiplication of 3 at each lattice module; but they do not contribute any computational complexity. The N-point FFT will be replaced by an N-point DCT, which has similar computational complexity to an N-point real FFT.
6. Parallel computer implementation Assume that a distributed-memory parallel system has M processors. A parallel implementation
C. Lu et al. 1 Signal Processing 57 (1997) 195-203
202 Table 1 Timing results on a 4 node IBM SP2 Size L
M
N
Time (ms)
16 384 65 536 262 144 589 824 1048 576 1327 104
128 256 512 768 1024 1152
128 256 512 168 1024 1152
3 14 101 229 441 834
and local; interprocessor communication is only required, therefore, at the last error correction stage of size N. The algorithms are suitable for processing sequential input data, and it is very attractive for parallel implementation. Timing results show that real-time signal processing is possible.
References Dl
of the block time-recursive algorithm can be described by the following procedures: - Broadcast to each processor the latest input block of size N. - Compute the difference between the sliding-in and sliding-out data blocks. - Perform point-wise multiplications by the precalculated e&r, i). - Compute N-point FFTs or DCTs. - Transfer data to the next-neighbor processor. - Perform point-wise additions with the previously calculated Gabor coefficients. Timing results of the time-recursive algorithm on a 4-node IBM SP2 configuration are given in Table 1, where L is the Gabor transform size, N is the sliding in/out block size for each update, and L = MN. The timing results given in Table 1 are the CPU time for each recursion. The simulation programs use double precision complex data and complex FFTs. Optimized FFT routines from the scientific libaries ESSL (SP-2) were used to compute the 1D FFTs. For real data each stage of the computation will be real or Hermitian, and computational complexity should be reduced approximately by 50%.
7. Conclusions In this paper, the block time-recursive Gabor transforms (BTRGTs) for both complex and real infinite signals are defined with window size L = MN. They are unified by a parallel lattice structure. Fast algorithms for their computation are presented with 1-D FFTs or DCTs at the center stage. Every stage of the computation is modular
M.G. Amin, “A new approach to recursive Fourier transform”, Proc. IEEE, Vol. 75, No. 11, November 1987, pp. 1537-1538. 121M. An, A. Brodzik, G. Kechriotis, C. Lu and R. Tolimieri, “Weyl-Heisenberg systems and the finite Zak transform”, Signal Processing, Submitted. r31 L. Auslander, I. Gertner and R. Tolimieri, “The finite Zak transforms and the finite Fourier transforms”, in: IMA on Radar and Sonar, Part II, Vol. 39, Springer, New York, 1991, pp. 21-36. M L. Auslander, I. Gertner and R. Tolimieri, “The discrete Zak transform application to time-frequency analysis and synthesis of non-stationary signals”, IEEE Trans. Signal Process., Vol. 39, No. 4, April 1991, pp. 825-835. WI M.J. Bastiaans, “Gabor signal expansion and degree of freedom of a signal”, Opt. Acta, Vol. 29, April 1982, pp. 8255835. C61 R.R. Bitmead, “On recursive discrete Fourier transform”, IEEE Trans. Acoust. Speech Signal Process., Vol. 30, No. 2, April 1982, pp. 319-322. c71 W. Chen, N. Kehtarnavaz and T.W. Spencer, “An efficient recursive algorithm for time-varying Fourier transform”, IEEE Trans. Signal Process., Vol. 41, No. 7, July 1993, pp. 2488-2490. VI D. Gabor, “Theory of communications”, J. IEE, Vol. 93, No. III, November 1946, pp. 429-457. c91 K.J.R. Liu and C.T. Chiu, “Unified parallel lattice structures for time-recursive discrete cosine/sine/Hartley transforms”, IEEE Trans. Signal Process., Vol. 41, No. 3, March 1993, pp. 1357-1377. UOI C. Lu, G. Kechriotis, M. An and R. Tolimieri, “The computation of Weyl-Heisenberg coefficients for critically sampled and oversampled signals”, Proc. ICSPAT’94, Dallas, TX, USA, 18-21 October 1994, pp. 824-829. Cl11 C. Lu, G. Kechriotis, M. An and R. Tolimieri, “A timerecursive algorithm for the computation of Gabor coefficients”, Proc. ICSPAT’95, Boston, MA, USA, 24-26 October 1995, pp. 642-646. 121 Y. Lu and J.M. Morris, “On novel critically-sampled Gabor transforms”, Proc. 29th Ann. Conf. on Inform. Science and Systems, The Johns Hopkins University, Baltimore, MD, March 1995, p. 95. C131 J.M. Morris and H. Xie, “Fast algorithms for generalized discrete Gabor expansion”, Signal Processing, Vol. 39, No. 3, September 1994, pp. 317-331.
C. Lu et al. 1 Signal Processing 57 (1997) 195-203 [14] S.H. Nawab and T.F. Quatieri, “Short-time Fourier transforms”, in: J.S. Lim and A.V. Oppenheim, eds., Advanced Topics in Signal Processing, Chapter 6, Prentice-Hall,
Englewood Cliffs, NJ, 1988. [lS] R.S. Orr, “The order of computation for finite discrete Gabor transforms”, IEEE Trans. Signal Process., Vol. 41, No. 1, January 1993, pp. 122-130. [16] S. Qian and J.M. Morris, “Wigner Distribution decomposition and cross-term deleted representation”, Signal Processing, Vol. 27, No. 2, May 1992, pp. 125-144.
203
[17] K.R. Rao and P. Yip, Discrete Cosine Transform: Algorithms, Advantages, Applications, Academic Press, New York, 1990. [18] D.F. Stewart and L.C. Potter, “Computationally attractive real Gabor transforms”, IEEE Trans. Signal Process., Vol. 43, No. 1, January 1995, pp. 77-83. [19] M. Zibulski and Y.Y. Zeevi, “Oversampling in the Gabor scheme”, IEEE Trans. Signal Process., Vol. 41, No. 8, August 1993, pp. 2679-2687.
PROCESSING ELSEVIER
Signal Processing 57 (1997) 195-203
Parallel lattice structure of block time-recursive generalized Gabor transforms’ Chao LIP, Sanjay Joshib, Joel M. Morrisb9* “Department of Computer & Information Sciences, Towson State University, Baltimore, MD 21204, USA bComputer Science and Electrical Engineering Department, University of Maryland Baltimore County, Baltimore MD 21250, USA
Received 19 April 1996; revised 3 October 1996
Abstract The problems of efficient computation of generalized Gabor transforms are considered in this work. Block timerecursive Gabor transforms are defined for both complex and real data sequences, fast algorithms are derived through two different approaches, one by the finite Zak transform, the other by biorthogonal functions. All the block timerecursive algorithms for the computation of generalized Gabor transforms can be unified by a parallel lattice structure. The unified lattice structure computes the transformed data sequentially with data sliding into/out-of a block of size N at a time. Gabor transform computation is based on the most recent L-point data, where L = MN. This model accommodates timely processing of the sequentially received data in digital signal transmission. Implementation on parallel computer system will be discussed. 0 1997 Elsevier Science B.V. Zusammenfassung In diesem Beitrag wird das Problem der etiienten Berechnung der verallgemeinerten Gabor Transformation aufgegriffen. Blockorientierte zeitrekursive Gabor Transformationen werden fiir komplexe und reelle Datenfolgen definiert. Zwei verschiedene Ansltze fiir schnelle Algorithmen werden abgeleitet, von denen der eine auf der Zak Transformation und der andere auf biorthogonalen Funktionen basiert. All diese blockorientierten zeitrekursiven Algorithmen zur Berechnung der verallgemeinerten Gabor Transformation lassen sich auf eine einheitliche parallele Lattice (Kreuzglied-)Struktur zuriickfiihren. Diese parallele Lattice Struktur berechnet die transformierten Daten sequentiell, indem jeweils BlGcke von N Daten gleichzeitig verarbeitet werden. Die Berechnung der Gabor Transformation wird jeweils fiir die letzte aktuellste Datenfolge der LInge L durchgeftihrt, wobei L = MN. Dieses Model1 paljt die zeitaufwendige Verarbeitung sequentiell empfangener Daten an die Erfordernisse der digitalen Signaliibertragung an. Die Implementierung auf einem Parallelrechnersystem wird diskutiert. 0 1997 Elsevier Science B.V. Rbumi! Ce travail examine les probl6mes du calcul efficient des transform&es de Gabor gCnCralis6es. Les transformtes de Gabor par blocs rCcursives dans le temps sont dCfinies pour des stquences de don&es complexes aussi bien que rCelles; des algorithmes rapides sont ensuite d&iv&s $ partir de deux approches diffkrentes, l’une par la transform&e de Zak finie,
*Corresponding ’ This research
author. Tel.: 1410455 3503; fax: 1410455 3969, e-mail: [email protected]. was supported by an ONR Grant # NOOO14-89-51210.
0165-1684/97/$17.00 0 PIISO165-1684(96)00195-8
1997 Elsevier Science B.V. All rights reserved.
196
C. Lu et al. / Signal Processing 57 (1997) 195-203
l’autre par les fonctions bi-orthogonales. Tous les algorithmes de calcul par blocs et rkcursifs dans le temps de transform&s de Gabor ginCralisCespeuvent etre unifiCsdans une structure parallkle en treillis. Cette structure en treillis unifike calcule les donnires transformbes SCquenciellementg l’aide de donnies glissant dans/hors d’un bloc de taille N g chaque instant du temps. Ce calcul de la transform&e de Gabor est basi: sur les donntes de L points les plus ricentes, oti L = MN. Ce modkle s’adapte au traitement opportun de donnkes regues &quenciellement dans les transmissions numkriques des signaux. L’implantation sur un ordinateur parallble sera discutke. 0 1997 Elsevier Science B.V. Keywords: Gabor transform; Time-recursive algorithms; Zak transform; FFT; Lattice structure
1. Introduction
A theoretical foundation for designing algorithms for computing Gabor coefficients at critical sampling was established in [3,4]. The Zak transform played a central role in establishing clear and easily computable conditions for the existence of the Gabor expansion and for the stability of computation. In [19], the Zak transform provided the framework for computing biorthogonals for rationally oversampled Weyl-Heisenberg systemforming frames. A similar approach for integer oversampled Weyl-Heisenberg systems can be found in [lo]. Computational complexity was studied in [lS] for the finite discrete Gabor transform, and several algorithms were compared in the work. Fast algorithms for computing Gabor coefficients, at critical sampling, integer oversampling as well as rationally oversampling, were given in [2, 10, 13, IS] for discrete and periodic signal x(k) with periodicity L such that x(k) = x(k + 15). Implementation of these algorithms were carried out on RISC processors and parallel computing environments, and timing results were reported in those papers. The main computational task in those algorithms was a, or a set of two-dimensional FFTs. Real Gabor transforms were defined in [12, 181, and computationally attractive algorithms based on the discrete cosine transform were given in [18]. Time-recursive structures have been developed for computing the discrete Fourier transform (DFT) [14], discrete cosine (DCT)/sine (DST)/ Hartley (DHT) transforms [9], as well as the timevarying Fourier transform (TVFT) [l, 6,7]. All of the time-recursive algorithms proposed in these papers share a common characteristic: when the
moving-window moves forward one step, a new sample is included while excluding the oldest sample. In [l 11, we proposed a block time-recursive algorithm for the computation of the Gabor transform in the critical-sampling case, which is based on the most recent L-point input data where L = MN. Data is sliding in/out a block of size N at a time. The time-recursive algorithm presented in [ 1l] was derived from the finite Zak transform. The main computational task is a set of 1-D FFTs followed by an error correction stage to the previously computed Gabor coefficients. In this paper, block time-recursive Gabor transforms (BTRGTs) are defined for both complex and real data sequences, and fast algorithms are derived through two different approaches: one by the finite Zak transform, and the other by biorthogonal functions. All the block time-recursive algorithms for the computation of generalized Gabor transforms can be unified by a parallel lattice structure. The unified lattice structure computes the transformed data sequentially with data sliding into/out-of in blocks of size N at a time. The building blocks of the parallel lattice structures are similar for both the complex and the real Gabor transforms regardless of which approach is used. They consist of a set of 1-D FFTs or 1-D DCTs on the difference of the incoming block and outgoing block of data, and the transformed data is then added to the previously computed Gabor transform. These algorithms decouple the transformed data components, and interprocessor communication is needed only at the last addition stage with N-point data transmitted to the nearest neighbor. Hence, they are very attractive for parallel implementation for real-time signal processing.
197
C. Lu et al. / Signal Processing 57 (1997) 195-203
The rest of the paper is organized as follows. In Section 2, the block time-recursive complex Gabor transform is defined and, in Section 3, the finite Zak transform approach is used to derive a fast algorithm. In Section 4, the block time-recursive real Gabor transform is defined, and a fast algorithm based on the DCT is given. In Section 5, the block time-recursive Gabor transforms (BTRGTs) are unified into an unified parallel structure and, in Section 6, parallel implementation of the algorithm on the IBM SP2 will be discussed. Finally, in Section 7 the work will be summarized.
2. Block time-recursive complex Gabor transform For an integer L = MN > 0, the generalized Gabor expansion of a complex-valued, finite, discrete-time signal x(k) is M-l N-l x(k) = 1 1 a(m, n)g(k - mN)ejzXknlN, m=o n=O
&=
[x0,
xl,
X2N-I
. . . ,XN-l,XN,
,X2N,
...
... ,
XL-IIT
,
(6)
= cxo, XI, .” >;YM-IIT, F is a block-diagonal
matrix with N-point FFT matrices on its diagonal,
1
[F(N) F=
I
F(N)
(7)
..
f’(N)
and r is a block circulant matrix with each block being a diagonal matrix. The diagonal elements from each block in the first row of r comprise the biorthogonal function
(1)
where g(k) is a discrete-time periodic, window sequence with period L = MN. Under certain conditions [2] the expansion coefficients u(m, n) can be found by
(8)
1 l-1
l-2
.. .
with each r,,, an N x N diagonal matrix
1
ry(mN) a(m, n) = $ Lil x(k)y(k
- mN)e-jzAknlN.
(2)
k-0
rm=
y(mN+ 1)
The discrete-time periodic sequence y(k) is the biorthogonal function to the window function g(k) satisfying L-l ,zo
g(k + mN)y(k)e-j2”k”‘N
= i&,6,.
(3)
Eq. (2) is called the generalized complex Gabor transform. The matrix formulation of Eq. (2) is given in [ 181 by 4 = FM,
y(mN
(9)
An algorithm was provided in [18] for the calculation of the biorthogonal function. If we write Eq. (4) in detail, we have go = F(N)[T,x, f1
=w)~r,-,~,
+ i-lx1 +
G..
+
rM_l~Y_l],
+ TOT1+ ... +rM-2+11,
(4)
where
gM-1 =F(N)[r,$+r2X1
4 = [a(O, O),~(0, l), . . . ,a(O, N - l), u(l,O), . . . ,
2$f-llT,
+ ..*+roxM_l]. (10)
u((M- 1, N - l)]’ = Lao, t1, ...
+ N - 1)
(5)
We define the Gabor transform at block time t in order to derive a block time-recursive relation
C. Lu et al. 1 Signal Processing 57 (1997) 195-203
198
between the updated Gabor transform and the previous computed Gabor transform. This is consistent with the new block of data sliding in and the old block of data sliding out, the block size is N, and the transform size is L = MN. Thus,
Denote the FZT of Gaussian function g as G( i, 1). An algorithm for the computation of Gabor coefficients of critically sampled signals of size L is given as follows [lo]: Step 1. Compute the FTZ X(i, 1) of signal x. Case I. G never vanishes in IMM,N. Step 2. Compute
al(t) = fxW~M-I$
+ ~0$+1 + ...
(14)
+ rM-$*+M- 11,
(11) whereZM,,=((i,1):0di
+
roxt+,-II,
and at block time t + 1, the Gabor transform can be expressed as ao(r + 1) = al(t) + F(NVM-r(Xt+M
- x,)9
ar(r + 1) = a&) + F(JWM-,(x*+,
- X,)>
C
i=O
I=0
O
CM-l@
+
1)
=
$I@)
+
F(W~o($+,
-
-
xtx
-I
i, I)
=
G(i’‘)’
(15)
G,441, 6, I) E k 6, I) E IM,N.
0,
$1.
Eq. (12) is called the block time-recursive Gabor transform (BTRGT). The Gabor coefficients A(t + 1) at block time t + 1 are given by the pre;iously computed A(t) at block time t with an error correction term - F(N)ri(Xt+M - X,), i = 0, 1, 2, , . . , A4 - 1. The same formulation can be derived through the finite Zak transform (FZT).
3. Block time-recursive via the FZT
O
XG, 1)
p( F(w~l(xt+M
P(i, l)wLMmlwNni,
Case II. Let 5 be the zero set of G in IM,N. If the signal x has its FZT X vanishing on i, then Step 2. Compute
(12) aM-z(r + I) = aM-l(t) +
M-l
a(m, n) = C
Step 3. Compute P( i, 1): N-l
the 2D M x N inverse FFT of
M-l
a(m, n) = C
C
i=O
I=0
O
P(i, l)~i”~w;~‘,
O
The block time-recursive Gabor (BTRGT) can be derived as follows: The FZT of x at time t,
Gabor transform
(17)
transform
M-l
X(i, 1, t) =
Define
c x(i + mN + tN)w$‘. m=o
(18)
The FZT of x at time t + 1: M-l
X(i, 1) =
1
x(i + mN)wg,
X(i, 1,t + 1) = f
m=O
i,lEZ,
wM=e
x(i + mN + tN)w$‘“-”
m=l - jZn/M
.
(13)
X(i, 1) is called the jinite Zak transform (FZT) of x.
=
X(i, 1, t)wk”I + [x(i + MN + tN) - x(i + tN)] w;‘.
(19)
199
C. Lu et al. / Signal Processing 57 (1997) 195-203
Compute P(i, 1,t + 1) =
X(i, 1,t + 1) G(i, 0
X(L 1,r) =mWM
-1
+ [x(i + MN + tN) - x(i + tN)]
wiMI
(20)
2
G(i, 1)
I
then the Gabor transform of x at t + 1, a(m, n, t + 1)
I
Err(m, n, t)
a(m +
1, n, 1) + Err(m, n, t)
I
c
o(m, n, t + 1)
rxti + MN
T.2 F.2
+zi=o
A-‘
+ tm - x(i + tmi
G(i, ;)
I=0
Fig. 1. Implementation flow diagram of Gabor coefficient computation on infinite data.
~
‘-
(21)
Step 2. Compute N-point FFTs on E,(m, i, t): N-l
Set
Err(m, n, t) = c M-l
e&n,
i) =
c
2!!!-
lzo
G(i,
E,(m, i, t)wNin,
i=O
-I -h
w 1)
M
(22)
’
Odm
(25)
Step 3. Add the error term Err(m, n, t) to Gabor coefficients at time t:
Then a(m, It, t + 1)
a(m, n, t + 1) = a(m + 1, n, t) + Err(m, n, t), N-l
=
a(m
+
1, n, t)
+
c
[x(i
+
AffN
+
O
tw
(26)
i=O -
x(i + tN)] e,(m, i)wivi”
(23)
When we write Eq.(23) in detail, it has the same format as Eq. (12), the block time-recursive complex Gabor transform. Procedure to compute the BTRGT at block time t+ 1: 1. Pre-compute e&m, i) and store it in memory. 2. Compute the first set Gabor coejicients at t = 0: a(m, n, 0). 3. Time-recursively compute: Step I. Compute
The implementation flow chart of the time-recursive algorithm is given in Fig. 1.
4. Block time-recursive real Gabor transform For an integer L = MN > 0, the generalized real Gabor expansion of the real finite discrete-time signal x(k) is defined as
x cos((2k + l)nn/2N)
E,(m, i, t) = [x(i + MN + tN) - x(i + tN)]e,(m, i), -!-
O
andO
(24)
‘3
1
b(m, O)g(k - mN) . (27)
200
C. Lu et al. / Signal Processing 57 (1997) 195-203
The matrix formulation for the computation expansion coefficients b(m, n) is given by [ 183
of
If we write Eq. (28) in detail, we have + Jrlx, _ + ... + JfM_1~,+1],
LJ,,= C(N)[f,x,, B = C&,
(28)
where Jj = [b(O, O),b(O, l), . . . ,b(0,iv - l), 41, O),
b, = C(N)JIJfM_lxo
+ foxI
+ Lf-23-11, ‘.
>
b((M - 1, N - l)lT
bypI
= Cbo, b1, .‘. >~hf-IIT,
= C(N)J[J&
+ f2g
-t ...
(29) + b$-
x=
+ ...
(35)
11.
[x,~x,,...,xN-l,xN,....X2N-l. XZN, '.. ,XL-IIT
(30)
= C~O,~I, .'.,$f-llT.
The matrix C is block-diagonal with C(N) the N-point discrete cosine transform (DCT) matrices on its main diagonal (N being an even integer). Thus,
WV C(N)J
C=
L
C(N)J
and the biorthogonal Jf,
f, JpMel
PO
12
(31)
We now define the real Gabor transform at block time t in order to derive a block time-recursive relation between the updated Gabor transform and the previous computed Gabor transform. This is consistent with the new block of data sliding in and the old block of data sliding out, the block size is N, and the transform size is L = MN. Thus, I b,(r) = C(N)CQ, + Jr,xt+l + ...
+ Jhf-&+M-11, !v~l(t) = C(N)JIJf’M-lxc
_!
+ &,+l
+ ...
+ LZ~t+M-Il,
function matrix p is given as ... ...
JpM+ fMm2
&l(t)
= C(N)J[J&,
+ &xt+l
+ ...
^
(32)
+ roxt+&f-1 I,
(36)
and at block time t + 1, the real Gabor transform can be given as with each f, is an N x N diagonal matrix
r?(mN) ?(mN
%, =
+
bo(r + 1) = WVC(C(N)J)-‘~I@)
+ Jhe
1)
I(X~+M - ?,)I,
b,(t+ 1) = C(N)J[C+(N)&(t) f(mN + N - 1)
+ L2(3+h4
- g1,
(33)
and J is the N x N row exchange matrix
ro ...
0
b~-z(t + 1)=
1-l
C(N)C(C(W-
lb-
+ Jh(xt+, (34)
10) - x,)1,
bMml(t + 1) = C(N)JIC-l(N)b,(t) +
rotxt+hl - x,)1.
(37)
C. Lu et al. / Signal Processing 57 (1997) 195-203
It can be proved that l-1
C(N)JCl(N)
=
0
0 -1
...
0
...
0
1
.
Lo
0
...
(38)
- 1 -I
Set s = C(N)JC’(N);
(39)
then we have b,(t - 1) = ~~1,(t) + C(N)JLf-1(x,+,
- x,)3
b,(t + 1) = @b,(t) + C(N)JL,(:,+,
- x,), Fig. 2. Parallel lattice structure of block time-recursive Gabor transform.
&$.f-,(t + 1) =
&f-,(t) + wvJ~&,+,
t)M-l@ + 1) = &o(t) + C(W~o(~*+,
- xt), - xc). (40)
Eq. (40) is called the block time-recursive real Gabor transform. The Gabor coefficients lj(t + 1) at block time t + 1 are given by the previously computed B(t) at block time t with an error correction term C(N)JTi(x,+,$$ - x,), i = 1,2, . . . , M - 1. Eq. (40) has a similar format to the block complex Gabor transform in time-recursive Eq. (12).
5. Parallel lattice structure of block time-recursive Gabor transforms Both the block time-recursive complex (Eq. (12)) and real (Eq. (40)) Gabor transforms have the same lattice structure. Their implementation by a parallel lattice array consisting of M lattice modules is illustrated in Fig. 2. It should be noted that data processing at each stage is by an array size of N. The latest sliding-in data block x,+~ of size N subtracts the sliding-out data block xt, and then the difference is multiplied by the biorthogonal coefficients Ti at each module. An N-point FFT (or DCT for real) is calculated at each module followed by an addition of shifted previously computed Gabor coefficients. An array of size N needs to be transferred to the nearest-neighbor module in circular order.
To obtain a block of N transformed data, and update the other M - 1 blocks, each module of the lattice array carries N-point multiplications, N-point additions, and an N-point FFT (or DCT). The N-point FFT requires an order of N log N additions and an order of N log N multiplications for the Cooley-Tukey FFT algorithm, or an order of N multiplications for the Good-Thomas prime factor algorithm combined with the Winograd algorithm for each relatively prime factor. If we assume the Cooley-Tukey FFT algorithm for the N-point FFT computation, then the computational complexity is of the order of MN + MNlog N.
(41)
For the block time-recursive real Gabor transform, the parallel lattice structure in Fig. 2 will include a rotation of J and a trivial multiplication of 3 at each lattice module; but they do not contribute any computational complexity. The N-point FFT will be replaced by an N-point DCT, which has similar computational complexity to an N-point real FFT.
6. Parallel computer implementation Assume that a distributed-memory parallel system has M processors. A parallel implementation
C. Lu et al. 1 Signal Processing 57 (1997) 195-203
202 Table 1 Timing results on a 4 node IBM SP2 Size L
M
N
Time (ms)
16 384 65 536 262 144 589 824 1048 576 1327 104
128 256 512 768 1024 1152
128 256 512 168 1024 1152
3 14 101 229 441 834
and local; interprocessor communication is only required, therefore, at the last error correction stage of size N. The algorithms are suitable for processing sequential input data, and it is very attractive for parallel implementation. Timing results show that real-time signal processing is possible.
References Dl
of the block time-recursive algorithm can be described by the following procedures: - Broadcast to each processor the latest input block of size N. - Compute the difference between the sliding-in and sliding-out data blocks. - Perform point-wise multiplications by the precalculated e&r, i). - Compute N-point FFTs or DCTs. - Transfer data to the next-neighbor processor. - Perform point-wise additions with the previously calculated Gabor coefficients. Timing results of the time-recursive algorithm on a 4-node IBM SP2 configuration are given in Table 1, where L is the Gabor transform size, N is the sliding in/out block size for each update, and L = MN. The timing results given in Table 1 are the CPU time for each recursion. The simulation programs use double precision complex data and complex FFTs. Optimized FFT routines from the scientific libaries ESSL (SP-2) were used to compute the 1D FFTs. For real data each stage of the computation will be real or Hermitian, and computational complexity should be reduced approximately by 50%.
7. Conclusions In this paper, the block time-recursive Gabor transforms (BTRGTs) for both complex and real infinite signals are defined with window size L = MN. They are unified by a parallel lattice structure. Fast algorithms for their computation are presented with 1-D FFTs or DCTs at the center stage. Every stage of the computation is modular
M.G. Amin, “A new approach to recursive Fourier transform”, Proc. IEEE, Vol. 75, No. 11, November 1987, pp. 1537-1538. 121M. An, A. Brodzik, G. Kechriotis, C. Lu and R. Tolimieri, “Weyl-Heisenberg systems and the finite Zak transform”, Signal Processing, Submitted. r31 L. Auslander, I. Gertner and R. Tolimieri, “The finite Zak transforms and the finite Fourier transforms”, in: IMA on Radar and Sonar, Part II, Vol. 39, Springer, New York, 1991, pp. 21-36. M L. Auslander, I. Gertner and R. Tolimieri, “The discrete Zak transform application to time-frequency analysis and synthesis of non-stationary signals”, IEEE Trans. Signal Process., Vol. 39, No. 4, April 1991, pp. 825-835. WI M.J. Bastiaans, “Gabor signal expansion and degree of freedom of a signal”, Opt. Acta, Vol. 29, April 1982, pp. 8255835. C61 R.R. Bitmead, “On recursive discrete Fourier transform”, IEEE Trans. Acoust. Speech Signal Process., Vol. 30, No. 2, April 1982, pp. 319-322. c71 W. Chen, N. Kehtarnavaz and T.W. Spencer, “An efficient recursive algorithm for time-varying Fourier transform”, IEEE Trans. Signal Process., Vol. 41, No. 7, July 1993, pp. 2488-2490. VI D. Gabor, “Theory of communications”, J. IEE, Vol. 93, No. III, November 1946, pp. 429-457. c91 K.J.R. Liu and C.T. Chiu, “Unified parallel lattice structures for time-recursive discrete cosine/sine/Hartley transforms”, IEEE Trans. Signal Process., Vol. 41, No. 3, March 1993, pp. 1357-1377. UOI C. Lu, G. Kechriotis, M. An and R. Tolimieri, “The computation of Weyl-Heisenberg coefficients for critically sampled and oversampled signals”, Proc. ICSPAT’94, Dallas, TX, USA, 18-21 October 1994, pp. 824-829. Cl11 C. Lu, G. Kechriotis, M. An and R. Tolimieri, “A timerecursive algorithm for the computation of Gabor coefficients”, Proc. ICSPAT’95, Boston, MA, USA, 24-26 October 1995, pp. 642-646. 121 Y. Lu and J.M. Morris, “On novel critically-sampled Gabor transforms”, Proc. 29th Ann. Conf. on Inform. Science and Systems, The Johns Hopkins University, Baltimore, MD, March 1995, p. 95. C131 J.M. Morris and H. Xie, “Fast algorithms for generalized discrete Gabor expansion”, Signal Processing, Vol. 39, No. 3, September 1994, pp. 317-331.
C. Lu et al. 1 Signal Processing 57 (1997) 195-203 [14] S.H. Nawab and T.F. Quatieri, “Short-time Fourier transforms”, in: J.S. Lim and A.V. Oppenheim, eds., Advanced Topics in Signal Processing, Chapter 6, Prentice-Hall,
Englewood Cliffs, NJ, 1988. [lS] R.S. Orr, “The order of computation for finite discrete Gabor transforms”, IEEE Trans. Signal Process., Vol. 41, No. 1, January 1993, pp. 122-130. [16] S. Qian and J.M. Morris, “Wigner Distribution decomposition and cross-term deleted representation”, Signal Processing, Vol. 27, No. 2, May 1992, pp. 125-144.
203
[17] K.R. Rao and P. Yip, Discrete Cosine Transform: Algorithms, Advantages, Applications, Academic Press, New York, 1990. [18] D.F. Stewart and L.C. Potter, “Computationally attractive real Gabor transforms”, IEEE Trans. Signal Process., Vol. 43, No. 1, January 1995, pp. 77-83. [19] M. Zibulski and Y.Y. Zeevi, “Oversampling in the Gabor scheme”, IEEE Trans. Signal Process., Vol. 41, No. 8, August 1993, pp. 2679-2687.