- Email: [email protected]
Parallel lattice structures of block time-recursive discrete Gabor transform and its inverse transform
ARTICLE IN PRESS
Signal Processing 88 (2008) 407–414 www.elsevier.com/locate/sigpro
Fast communication
Parallel lattice structures of block time-recursive discrete Gabor transform and its inverse transform Liang Taoa,, Hon Keung Kwanb a
b
School of Computer Science and Technology, Anhui University, Hefei, Anhui 230039, P.R. China Department of Electrical and Computer Engineering, University of Windsor, 401 Sunset Avenue, Windsor, Ontario, Canada N9B 3P4 Received 15 March 2007; received in revised form 17 July 2007; accepted 7 August 2007 Available online 15 August 2007
Abstract The traditional discrete Gabor transform (DGT) is briefly reviewed in this paper. Block time-recursive algorithms for the efficient and fast computation of the DGT coefficients and for the fast reconstruction of the original signal from the coefficients are then developed in both the critical sampling case and the over-sampling case. Unified parallel lattice structures for the implementation of the algorithms are studied. Computational complexity analysis and comparison show that the proposed algorithms provide a more efficient and faster approach as compared to the existing DGT algorithms for the computation of DGTs. r 2007 Published by Elsevier B.V. Keywords: Discrete Gabor transform; Parallel lattice structure; Block time-recursive algorithms; Discrete Fourier transform (DFT)
1. Introduction The Gabor transform was first proposed by Gabor [1] to represent a signal in both time and frequency domains. Although the Gabor transform has been recognized as being very useful in diverse areas such as speech and image processing; radar, sonar and seismic data processing and interpretation, its real time applications are limited due to difficulties associated with the high complexity of the computation of the discrete Gabor transform (DGT) coefficients and the reconstruction of the original signal from the transform coefficients. Therefore, many approaches have been proposed Corresponding author.
E-mail addresses: [email protected] (L. Tao), [email protected] (H.K. Kwan). 0165-1684/$ - see front matter r 2007 Published by Elsevier B.V. doi:10.1016/j.sigpro.2007.08.005
to solve the problem based on finding an analytic solution [2–7]. In 1997, the parallel lattice structure of block time-recursive DGTs was developed by Lu et al. [8], which, so far, has been the fastest algorithm among the existing DGT algorithms using the analytic solution method. However, the algorithm is limited to the critical sampling case and the computation of the DGT coefficients. The oversampling case and the reconstruction of the original signal from the transform coefficients (i.e., the inverse DGT) are not included in [8]. In this paper, we shall show that the method used in [8] can be generalized to the over-sampling case, and the unified parallel lattice structures of block timerecursive DGT can be developed, not only for the computation of the DGT coefficients but also for the reconstruction of the original signal from the coefficients.
ARTICLE IN PRESS L. Tao, H.K. Kwan / Signal Processing 88 (2008) 407–414
408
The outline of the paper is organized as follows. Section 2 briefly introduces the traditional DGT. The block time-recursive DGT algorithms for the fast computation of the DGT coefficients and for the fast reconstruction of the original signal from the coefficients are developed, respectively, in Sections 3 and 4. In Section 5, the two algorithms will be, respectively, implemented by a unified parallel lattice structure, and the computational complexity analysis and comparison of the proposed algorithms are provided. Numerical simulations are given in Section 6. Finally, we make our conclusions in Section 7.
2. Review of DGT Let x(k) denote a real finite and periodic discretetime signal with a period L, the traditional discrete Gabor expansion [3,4] is defined by
xðkÞ ¼
M1 X N1 X
¯ expðj2pnk=NÞ, aðm; nÞgðk mNÞ
m¼0 n¼0
(1) pffiffiffiffiffiffiffi where j ¼ 1 and the coefficients a(m, n) can be obtained by
aðm; nÞ ¼
L1 X
following biorthogonality condition [3]: L1 X
¯ gðk þ mNÞ expðj2pnk=NÞgðkÞ ¼ ðL=MNÞdm dn
k¼0
¯ 1; 0pnpN ¯ 1, 0pmpM
ð3Þ
where dk denotes the Kronecker delta. Note that the DGT exists only if the biorthogonality condition is satisfied. 3. Block time-recursive DGT algorithm for the computation of the DGT coefficients In the over-sampling case or the critical sampling ¯ ¼ MN ¯ case, LpMN and L ¼ M N should be satisfied, thus, the over-sampling rate b ¼ MN=L ¼ ¯ ¼ N=N. ¯ Obviously, b ¼ 1 and b41 are M=M corresponding to the critical sampling and the oversampling, respectively. Properly choose the sampling ¯ or N and N, ¯ to make the overparameters M and M, sampling rate b a positive integer. Suppose that ¯ 1, and i ¼ 0; 1; . . . ; b 1, j ¼ 0; 1; . . . ; M xj ¼ ½xðjNÞ; xðjN þ 1Þ; . . . ; xðjN þ N 1ÞT ,
(4)
x ¼ ½xT0 ; xT1 ; . . . ; xTM1 T , ¯
(5)
aij ¼ ½aði þ jb; 0Þ; aði þ jb; 1Þ; ; aði þ jb; N 1ÞT , (6)
¯ expðj2pnk=NÞ, xðkÞgðk mNÞ
(2)
k¼0
where Eq. (2) is termed the DGT [3,4] and Eq. (1) can also be called the inverse DGT. In (1) and (2), ¯ ¯ M and N are the numbers of L ¼ NM ¼ N M, ¯ sampling points in time and frequency domains. M ¯ and N are the frequency and time sampling intervals, respectively. The condition MNXL must be satisfied for a stable reconstruction. The critical ¯M ¯ ¼ NM ¼ L (The numsampling occurs when N ber of coefficients a(m, n) is equal to the number of time samples x(k)), and the over-sampling occurs when MN4L. There may be a loss of information in an under-sampling condition (MNoL). Because the critical sampling leads to bad localization of the analysis window g(k), the over-sampling DGT is often used in signal analysis and processing such as noise reduction for MNR FID signals [9]. Note that the synthesis window g(k) and the analysis window g(k) are all real and periodic in L, and satisfy the
Ai ¼ ½ðai0 ÞT ; ðai1 ÞT ; . . . ; ðaiM1 ÞT T . ¯
(7)
Now let m ¼ i þ jb, (2) can be rewritten as aði þ jb; nÞ ¼
L1 X
¯ jNÞ xðkÞgðk iN
k¼0
expðj2pnk=NÞ
ð8Þ
and the matrix formation of (8) is given by Ai ¼ FCi X
(9)
¯ N-point where F is a block-diagonal matrix with M fast discrete Fourier transform (DFT) matrices on its diagonal 3 2 FðNÞ 7 6 FðNÞ 7 6 7 6 (10) F¼6 7 .. 7 6 . 5 4 FðNÞ
ARTICLE IN PRESS L. Tao, H.K. Kwan / Signal Processing 88 (2008) 407–414
and Ci is a block circulant matrix 3 2 i Ci1 CiM1 C0 ¯ 7 6 i Ci0 CiM2 7 6 CM1 ¯ 7 6 ¯ i , C ¼6 . .. 7 .. .. 6 .. . 7 . . 5 4 Ci1 Ci2 . . . Ci0
and at block time t+1, the DGT coefficients can be expressed as
(11)
with each block Cij being an N N diagonal matrix: 2 6 6 6 Cij ¼ 6 6 4
3
¯ gðjN iNÞ ¯ þ 1Þ gðjN iN ..
. ¯ þ N 1Þ gðjN iN
7 7 7 7. 7 5
(12) The diagonal elements from each block in the first row of Ci comprise one period of the analysis window gðkÞ. Eq. (9) can be written in detail as follows: ai0 ¼ FðNÞ½Ci0 x0 þ Ci1 x1 þ þ CiM1 xM1 ; ¯ ¯ x0 þ Ci0 x1 þ þ CiM2 xM1 ; ai1 ¼ FðNÞ½CiM1 ¯ ¯ ¯ .. . aiM1 ¯
¼
FðNÞ½Ci1 x0
þ
Ci2 x1
þ þ
Ci0 xM1 : ¯
Suppose the block xt is inputted into a block ¯ 1 in series (see delay line from t ¼ 0 to t ¼ M Fig. 1). The time interval of the series inputs of blocks is defined as a block-time interval. We define the DGT at block time t (t ¼ 0; 1; 2; . . .) in order to derive a block time-recursive relation between the updated DGT coefficients and the previous computed DGT coefficients. 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 ¯ ¼ MN. ¯ transform size is L ¼ M N Thus, x¯ tþM1 ai0 ðtÞ ¼ FðNÞ½Ci0 x¯ t þ Ci1 x¯ tþ1 þ þ CiM1 ; ¯ ¯ x¯ t þ Ci0 x¯ tþ1 þ þ CiM2 x¯ tþM1 ai1 ðtÞ ¼ FðNÞ½CiM1 ; ¯ ¯ ¯ .. . aiM1 ðtÞ ¼ FðNÞ½Ci1 x¯ t þ Ci2 x¯ tþ1 þ þ Ci0 x¯ tþM1 ¯ ¯
xt+M-1
···
xt+1
Fig. 1. A block delay line.
xt
ðx¯ tþM¯ x¯ t Þ; ai1 ðt þ 1Þ ¼ ai2 ðtÞ þ FðNÞCiM2 ¯ .. . aiM2 ðt ¯
þ 1Þ ¼
(15)
þ FðNÞCi1 ðx¯ tþM¯ x¯ t Þ; FðNÞCi0 ðx¯ tþM¯ x¯ t Þ:
aiM1 ðtÞ ¯
aiM1 ðt þ 1Þ ¼ ai0 ðtÞ þ ¯
Eq. (15) is called the block time-recursive DGT for the computation of the DGT coefficients. The DGT coefficients Ai ðt þ 1Þ at block time t+1 are given by Ai ðtÞ at block time t with an error correction FðNÞCij ðx¯ tþM¯ x¯ t Þ, i ¼ 0; 1; . . . ; b 1, j ¼ 0; 1; . . . ; ¯ 1. The error correction becomes zero when M ¯ due to the periodicity of x(k), and then the t ¼ M, time-recursive process in (15) gets into stable stage. Note that x¯ t ¼ 0 in the data delay line when ¯ 1. tpM Eq. (15) can also be rewritten as follows:
ðx¯ tþM¯ x¯ t Þ; a~ i1 ðt þ 1Þ ¼ a~ i2 ðtÞ þ CiM2 ¯ .. . a~ iM2 ðt þ 1Þ ¼ a~ iM1 ðtÞ þ Ci1 ðx¯ tþM¯ x¯ t Þ; ¯ ¯
(16)
ðt þ 1Þ ¼ a~ i0 ðtÞ þ Ci0 ðx¯ tþM¯ x¯ t Þ a~ iM1 ¯ ¯ compute and when t ¼ M, aij ¼ FðNÞ~aij
¯ 1. i ¼ 0; 1; . . . ; b 1; j ¼ 0; 1; . . . ; M
(17) Obviously, using (16) and (17) instead of (15) will ¯ 1Þ ðNpoint fast DFTÞ computation. save MðM Note that Ref. [8] only used the formula similar to (15) for the computation of the transform coefficients in the critical sampling case. To summarize, a procedure of the above algorithm is given as follows. Algorithm 1.
(14) xt+M
ðx¯ tþM¯ x¯ t Þ; ai0 ðt þ 1Þ ¼ ai1 ðtÞ þ FðNÞCiM1 ¯
ðx¯ tþM¯ x¯ t Þ; a~ i0 ðt þ 1Þ ¼ a~ i1 ðtÞ þ CiM1 ¯
(13)
xt
409
(a) Let t ( 0, clear all the block units in the block delay line and clear all of a~ ij ði ¼ 0; 1; . . . ; b 1; ¯ 1Þ to be zero. j ¼ 0; 1; . . . ; M (b) Input the signal block xt into the block delay line to let x¯ tþM¯ ( xt .
ARTICLE IN PRESS L. Tao, H.K. Kwan / Signal Processing 88 (2008) 407–414
410
¯ 1, com(c) For i ¼ 0; 1; . . . ; b 1; j ¼ 0; 1; . . . ; M pute
a ti
i at+M
x¯ tþM¯ a~ ij ðt þ 1Þ ¼ a~ ijþ1 ðtÞ þ CiMj1 ¯
i at+M-1
···
i at+1
a ti
Fig. 2. A block delay line.
¯ 1, bring all of Note that a~ iM¯ ¼ a~ i0 . If t ¼ M ¯ 1Þ to Step ¼ 0; 1; . . . ; b 1; j ¼ 0; 1; . . . ; M (d), otherwise let t ( t þ 1 and go to Step (b). ¯ 1, com(d) For i ¼ 0; 1; . . . ; b 1; j ¼ 0; 1; . . . ; M pute a~ ij ði
aij ¼ FðNÞ~aij At present, we obtain the DGT coefficients A ¼ ½Ai b1 i¼0 by Eq. (7).
The diagonal elements from each block in the first row of G i comprise one period of the synthesis window g(k). Eq. (19) can be written in detail as follows:
x0 ¼
b1 X i ¯ i i ¯ ½G i0 FðNÞa 0 þ G 1 FðNÞa1 i¼0 i ¯ FðNÞa þ þ G iM1 , ¯ ¯ M1
4. Block time-recursive DGT algorithm for the reconstruction of the original signal
x1 ¼
b1 X i ¯ i i ¯ FðNÞa ½G iM1 ¯ 0 þ G 0 FðNÞa1 i¼0
¯ ¼ N; bM ¯ ¼ M, Let m ¼ i þ jb, and note that bN Eq. (1) can be rewritten as xðkÞ ¼
¯ b1 M1 N 1 X X X
xM1 ¼ ¯
2pnk exp j N ¼
¯ jNÞ gðk iN
i¼0 j¼0
2pnk exp j . N
i ¯ þ þ G i0 FðNÞa . ¯ M1 N 1 X
aði þ jb; nÞ
n¼0
ð18Þ
(19)
i¼0
¯ where F¯ ¼ F 1 with each diagonal block FðNÞbeing ¯ the inverse matrix of F(N) (i.e., FðNÞ is an N-point fast inverse DFT), G i is also a block circulant matrix: 3 2 i G0 G i1 G iM1 ¯ 7 6 i G i0 G iM2 7 6 G M1 ¯ 7 6 ¯ i (20) G ¼ 6. 7, .. .. .. 7 6 .. . . . 5 4 G i1 G i2 . . . G i0 with each block 2
6 6 6 i Gj ¼ 6 6 4
b1 X i ¯ i i ¯ ½G i1 FðNÞa 0 þ G 2 FðNÞa1 i¼0
The matrix formation of (18) is given by b1 X ¯ i, X¼ G i FA
G ij
.. .
¯ jNÞ aði þ jb; nÞgðk iN
i¼0 j¼0 n¼0
¯ b1 M1 X X
i ¯ FðNÞa þ þ G iM2 , ¯ ¯ M1
As in the block time-recursive algorithm for the computation of the DGT coefficients, we define the DGT at block time t in order to derive a block timerecursive relationship between the updated reconstructed signal and the previous computed reconstructed signal. This is consistent with the new block of the DGT coefficients sliding in and the old block of the DGT coefficients sliding out (see Fig. 2), the block size is N, and the transform size is L ¼ ¯ ¼ MN. ¯ MN Thus,
x0 ðtÞ ¼
¯ FðNÞ¯ aitþM1 þ þ G iM1 , ¯ ¯ x1 ðtÞ ¼
¯ þ 1Þ gðjN iN ..
. ¯ þ N 1Þ gðjN iN
(21)
b1 X ¯ ¯ FðNÞ¯ ait þ G 0 FðNÞ¯ aitþ1 ½G M1 ¯ i¼0
¯ FðNÞ¯ aitþM1 þ þ G M2 , ¯ ¯ .. .
3
7 7 7 7. 7 5
b1 X ¯ ¯ ait þ G i1 FðNÞ¯ aitþ1 ½G i0 FðNÞ¯ i¼0
being an N N diagonal matrix,
¯ gðjN iNÞ
ð22Þ
xM1 ðtÞ ¼ ¯
b1 X ¯ ¯ ait þ G i2 FðNÞ¯ aitþ1 ½G i1 FðNÞ¯ i¼0
¯ aitþM1 þ þ G i0 FðNÞ¯ ¯
ð23Þ
ARTICLE IN PRESS L. Tao, H.K. Kwan / Signal Processing 88 (2008) 407–414
and at block time t+1, the reconstructed signal can be expressed as x0 ðt þ 1Þ ¼ x1 ðtÞ þ
b1 P i¼0
x1 ðt þ 1Þ ¼ x2 ðtÞ þ
b1 P i¼0
¯ FðNÞð¯ aitþM¯ a¯ it Þ; G iM1 ¯
(b) For i ¼ 0; 1; . . . ; b 1, input the DGT coefficient block ait into the ith block delay line to let a¯ itþM¯ ( ait . ¯ 1, compute (c) For j ¼ 0; 1; . . . ; M xj ðt þ 1Þ ¼ xjþ1 ðtÞ þ
¯ FðNÞð¯ aitþM¯ a¯ it Þ; G iM2 ¯
xM1 ðt þ 1Þ ¼ x0 ðtÞ þ ¯
b1 P i¼0
b1 P
¯ aitþM¯ a¯ it Þ; G i1 FðNÞð¯
¯ aitþM¯ a¯ it Þ: G i0 FðNÞð¯
i¼0
Eq. (24) is called the block time-recursive DGT for the reconstruction of the original signal. The reconstructed signal Xðt þ 1Þ at block time t+1 is given by XðtÞ at block time t with an error correction ¯ aitþM¯ a¯ it Þ; G ij FðNÞð¯
¯ FðNÞ¯ aitþM¯ . G iMj1 ¯
i¼0
(24)
b1 X
b1 X
¯ 1, we obtain Note that xM¯ ¼ x0 . If t ¼ M the reconstructed signal X by Eq. (5), otherwise, let t ( t þ 1 and go to Step (b).
.. . xM2 ðt þ 1Þ ¼ xM1 ðtÞ þ ¯ ¯
411
¯ 1. j ¼ 0; 1; . . . ; M
i¼0
To summarize, a procedure of the above algorithm is given as follows. Algorithm 2. (a) Let t ( 0, clear all the block units in the block delay line and clear all of xj ðj ¼ 0; 1; . . . ; ¯ 1; Þ, to be zero. M
5. Unified parallel lattice structures of block timerecursive DGT and its inverse transform 5.1. Unified parallel lattice structures The block time-recursive Eqs. (16) and (17) for the computation of the DGT coefficients can be implemented by b unified parallel lattice arrays as shown in Fig. 3. Each of the parallel lattice arrays ¯ unified lattice modules. It should be consists of M noted that data processing at each stage is by an array size N. The latest sliding-in data block x¯ tþM¯ of size N subtracts the sliding-out data block x¯ t , and then the difference is multiplied by the analysis window block Cij at each module followed by an addition of shifted and previously computed data. A block data of size N needs to be transferred to the nearest-neighbor module in circular order. Finally, with one period of signal sliding in, an N-point fast DFT is computed at the end of each module.
xt i M-1
xt+M
a 0i
i M-2
+
F (N)
a 1i
1i
+
F (N)
i a M-2
F (N)
i a M-1
+
-1 xt+1
F (N) Z-1
xt+M-1 . . .
+
Z-1
Z-1 xt
0i
+ Z-1
Fig. 3. Unified parallel lattice structure of block time-recursive DGT for the computation of the DGT coefficients (i ¼ 0; 1; . . . ; b 1).
ARTICLE IN PRESS L. Tao, H.K. Kwan / Signal Processing 88 (2008) 407–414
412 a ti
i=0, 1, · · ·, −1 i
i
G M−1
at+M
+
+
x0
i at+M-1
+
F (N)
i=0, 1, · · ·, −1 Z-1 i G1
+
+
+
+
xM-2
-1 i at+1
i=0, 1, · · ·, −1 Z-1
a ti
i
G0
xM-1 Z-1
Fig. 4. Unified parallel lattice structure of block time-recursive DGT for the reconstruction of the original signal from the DGT coefficients.
The block time-recursive Eq. (24) for the reconstruction of the original signal can also be implemented by b unified parallel lattice arrays. Each ¯ unified lattice modules, which of them consists of M is illustrated in Fig. 4. Data processing at each stage is by an array size N. The latest sliding-in data block a¯ itþM¯ of size N subtracts the sliding-out data block a¯ it , and then an N-point fast inverse DFT of the difference is applied. The resulting data is multiplied by the synthesis window block G ij at each module followed by b1 additions of the resulting block data from the other b1 parallel lattice arrays, and one addition of shifted and previously computed reconstructed signal. A block data of size N needs to be transferred to the nearest-neighbor (output) module in circular order. Obviously, the oversampling case shown in Figs. 3 and 4 will become the critical sampling case when ¯ and N ¼ N ¯ (i.e., b ¼ 1). M¼M 5.2. Computation time and complexity analysis Since the total computational complexity in a parallel algorithm is equally shared by each unified parallel lattice module, the computation time of the parallel algorithm depends on the computational complexity of each unified parallel lattice module, while the computation time of a series algorithm depends on its total computational complexity. In Fig. 3, each parallel module carries N-point multiplications and N-point additions at each block time, and an N-point fast DFT (with computational complexity 0:5N log2 N) at the last block time.
Table 1 Comparison of computational complexity related to total computational time References
Computational complexity
Applicability
[4] [6] [7] [8] Proposed algorithms
ML þ 0:5ML log2 N ML+0.5L log2 N 0.5L log2(4L)+L log2 M L þ 0:5L log2 N L þ 0:5N log2 N L þ 0:5L log2 N
CS, CS, CS, CS, CS, CS,
OS, Coef Coef, Reco Coef, Reco Coef OS, Coef OS, Reco
Therefore, to obtain a block of N transformed data, the computational complexity is of the order of L þ 0:5N log2 N.
(25)
In Fig. 4, at each block time, every parallel module of each parallel array carries an N-point fast inverse DFT, N-point multiplications followed by b Npoint additions. Therefore, to obtain a block of N reconstructed data, the computational complexity is of the order of L þ 0:5L log2 N.
(26)
Table 1 gives a comparison between the proposed algorithms and the main existing DGT algorithms. In the table, the symbols CS, OS, Coef and Reco indicate four cases of algorithm applicability which are respectively the critical sampling case, the oversampling case, the computation of the DGT coefficients, and the reconstruction of the original signal from the coefficients (i.e., the inverse DGT).
ARTICLE IN PRESS L. Tao, H.K. Kwan / Signal Processing 88 (2008) 407–414
We can see from the table that the computational complexity in (25) or (26) related to the total computation time of the proposed parallel algorithms is far less than that of the existing DGT algorithms. Thus, the proposed algorithms are more efficient and faster in computing the DGT coefficients and reconstructing the original signal from the coefficients.
413
DGT coefficients |a(m,n)|
Fourier spectrum Fx (f) 8 kHz
127
4 kHz
0
0 kHz 0
127
0
70
A speech signal x (k) 1
6. Computer simulations
0
In this section, we present some numerical simulation results, which are computed by the proposed algorithms. In the oversampling case ¯ ¼M ¯ ¼ 16, b ¼ 8), the top plot (N ¼ M ¼ 128, N of Fig. 5 depicts a Gaussian synthesis window, qffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffi gðkÞ ¼ 2=40 expfp½ðk 1023:5Þ=402 g, with length L ¼ 2048, corresponding to its analysis window g (k) in the bottom plot of Fig. 5, which is computed by Eq. (3). g(k) and g(k) are used to compute the DGT coefficients of a speech signal x(k) and to reconstruct the original speech signal by the DGT coefficients. The speech signal x(k) is shown in Fig. 6 with length L ¼ 2048, i.e., 256 ms of the word ‘‘yes’’ spoken by a man. x(k) is partitioned by Eq. (5) into 16 nonoverlapping block vectors of ¯ 1Þ. Correspondlength 128 ðxj ; j ¼ 0; 1; . . . ; M ingly, the DGT coefficient matrix a(m, n) of the speech signal is also partitioned by Eq. (7) into ¯ nonoverlapping block vectors 8 16 (i.e., b M) of length 128 ðaij ; i ¼ 0; 1; . . . ; b 1; j ¼ 0; 1; . . . ; ¯ 1Þ. By using the proposed algorithm for the M computation of the DGT coefficients, Fig. 6 displays the computed DGT coefficients (magnitude) of x(k). In order to compare with the DGT coefficients, Fig. 6 also shows one period of the
g (k)
0.2 0.1 0 0
512
1024 k
1536
2048
0
512
1024
1536
2048
(k)
0.04 0.02 0 k Fig. 5. Gaussian window g(k) and its analysis window g (k). ¯ ¼ 16 (oversampling). L ¼ 2048, N ¼ 128, N
-1 0 ms
256 ms
Fig. 6. DGT coefficients jaðm; nÞj and Fourier spectrum Fx(f) of a speech signal x(k).
Table 2 Comparison of total number of multiplications related to total computational time References
Total number of multiplications (applicability)
[4] [6] [7] [8] Proposed algorithms
786432 (CS, Coef), 1179648 (OS, Coef) 266240 (CS, Coef, Reco) 27648 (CS, Coef, Reco) 6144 (CS, Coef) 2080 (CS, Coef), 2496(OS, Coef) 6144 (CS, Reco), 9216(OS, Reco)
Fourier spectrum Fx(f) of x(k), which reveals four frequency tones in the range from 0 to 4 kHz (half the sampling frequency 8 kHz). From the spectrum alone, we can not tell how those frequencies change over time; however, from the DGT coefficients, not only can we see how frequencies (n) change over time (m), but we can also see the intensity of the frequencies as shown by the relative gray levels of the image. The proposed algorithm for the reconstruction of the original signal is also simulated. The signal reconstruction error (MSE) is around 10–15, which is virtually error-free reconstruction. A numeric comparison on the total number of multiplications related to the total computational time between the proposed algorithms and the main existing DGT algorithms is given in Table 2, which further indicates the advantages of the proposed algorithms. In Table 2, the parameters L ¼ 2048, M ¼ N ¼ 128 for oversampling and L ¼ 2048, M ¼ 128, N ¼ 16 for critical sampling are used. For example, ‘‘6144 (CS, Coef)’’ for Ref. [8] means that, for computing the DGT coefficients in the
ARTICLE IN PRESS 414
L. Tao, H.K. Kwan / Signal Processing 88 (2008) 407–414
critical sampling case (L ¼ 2048, M ¼ 128, N ¼ 16), the total number of multiplications of the algorithm in Ref. [8] is 6144, which is calculated using the formula of Ref. [8, Table 1].
Talent Development Foundation under Grant no. 2005Z029, the Anhui University Program for Creative Teams, and the MOE’s Key Laboratory of Intelligent Computing and Signal Processing in Anhui University.
7. Conclusions In this paper, we developed two block timerecursive DGT algorithms in both the critical sampling case and the oversampling case, not only for the computation of the DGT coefficients but also for the reconstruction of the original signal from the coefficients. The two algorithms were implemented respectively by a unified parallel lattice structure. Computational complexity analysis and comparison showed that the proposed algorithms are more efficient and faster than the existing DGT algorithms in computing the DGT coefficients and reconstructing the original signal from the coefficients. Moreover, the proposed algorithms can also be generalized to the DHT (discrete Hartley transform)-based real-valued DGT, which was presented in our previous work [10]. Acknowledgments The authors would like to thank the anonymous reviewers for their constructive comments. The first author would like to acknowledge the supports of the National Natural Science Foundation of China under Grant no. 60572128, the Anhui Provincial
References [1] D. Gabor, Theory of communication, J. Inst. Electr. Eng. 93 (3) (1946) 429–457. [2] M. Bastiaans, Gabor’s expansion of a signal into Gaussian elementary signals, Opt. Eng. 20 (4) (1981) 594–598. [3] J. Wexler, S. Raz, Discrete Gabor expansions, Signal Processing 21 (3) (1990) 207–220. [4] S. Qian, D. Chen, Discrete Gabor transforms, IEEE Trans. Signal Process. 21 (7) (1993) 2429–2438. [5] S. Qian, D. Chen, Joint time–frequency analysis, IEEE Signal Process. Mag. 6 (2) (1999) 52–67. [6] L. Wang, C. Chen, W. Lin, An efficient algorithm to compute the complete set of discrete Gabor coefficients, IEEE Trans. Image Process. 3 (1) (1994) 87–92. [7] S. Qiu, F. Zhou, P.E. Crandall, Discrete Gabor transforms with complexity O(Nz log N), Signal Processing 77 (2) (1999) 159–170. [8] C. Lu, S. Joshi, J.M. Morris, Parallel lattice structure of block time-recursive generalized Gabor transforms, Signal Processing 57 (2) (1997) 195–203. [9] B. Friedlander, A. Zeira, Oversampled Gabor representation for transient signals, IEEE Trans. Signal Process. 43 (9) (1995) 2088–2094. [10] L. Tao, H.K. Kwan, Block time-recursive real-valued discrete Gabor transform implemented by unified parallel lattice structures, IEICE Trans. Inform. Systems E88-D (7) (July 2005) 1472–1478.
Signal Processing 88 (2008) 407–414 www.elsevier.com/locate/sigpro
Fast communication
Parallel lattice structures of block time-recursive discrete Gabor transform and its inverse transform Liang Taoa,, Hon Keung Kwanb a
b
School of Computer Science and Technology, Anhui University, Hefei, Anhui 230039, P.R. China Department of Electrical and Computer Engineering, University of Windsor, 401 Sunset Avenue, Windsor, Ontario, Canada N9B 3P4 Received 15 March 2007; received in revised form 17 July 2007; accepted 7 August 2007 Available online 15 August 2007
Abstract The traditional discrete Gabor transform (DGT) is briefly reviewed in this paper. Block time-recursive algorithms for the efficient and fast computation of the DGT coefficients and for the fast reconstruction of the original signal from the coefficients are then developed in both the critical sampling case and the over-sampling case. Unified parallel lattice structures for the implementation of the algorithms are studied. Computational complexity analysis and comparison show that the proposed algorithms provide a more efficient and faster approach as compared to the existing DGT algorithms for the computation of DGTs. r 2007 Published by Elsevier B.V. Keywords: Discrete Gabor transform; Parallel lattice structure; Block time-recursive algorithms; Discrete Fourier transform (DFT)
1. Introduction The Gabor transform was first proposed by Gabor [1] to represent a signal in both time and frequency domains. Although the Gabor transform has been recognized as being very useful in diverse areas such as speech and image processing; radar, sonar and seismic data processing and interpretation, its real time applications are limited due to difficulties associated with the high complexity of the computation of the discrete Gabor transform (DGT) coefficients and the reconstruction of the original signal from the transform coefficients. Therefore, many approaches have been proposed Corresponding author.
E-mail addresses: [email protected] (L. Tao), [email protected] (H.K. Kwan). 0165-1684/$ - see front matter r 2007 Published by Elsevier B.V. doi:10.1016/j.sigpro.2007.08.005
to solve the problem based on finding an analytic solution [2–7]. In 1997, the parallel lattice structure of block time-recursive DGTs was developed by Lu et al. [8], which, so far, has been the fastest algorithm among the existing DGT algorithms using the analytic solution method. However, the algorithm is limited to the critical sampling case and the computation of the DGT coefficients. The oversampling case and the reconstruction of the original signal from the transform coefficients (i.e., the inverse DGT) are not included in [8]. In this paper, we shall show that the method used in [8] can be generalized to the over-sampling case, and the unified parallel lattice structures of block timerecursive DGT can be developed, not only for the computation of the DGT coefficients but also for the reconstruction of the original signal from the coefficients.
ARTICLE IN PRESS L. Tao, H.K. Kwan / Signal Processing 88 (2008) 407–414
408
The outline of the paper is organized as follows. Section 2 briefly introduces the traditional DGT. The block time-recursive DGT algorithms for the fast computation of the DGT coefficients and for the fast reconstruction of the original signal from the coefficients are developed, respectively, in Sections 3 and 4. In Section 5, the two algorithms will be, respectively, implemented by a unified parallel lattice structure, and the computational complexity analysis and comparison of the proposed algorithms are provided. Numerical simulations are given in Section 6. Finally, we make our conclusions in Section 7.
2. Review of DGT Let x(k) denote a real finite and periodic discretetime signal with a period L, the traditional discrete Gabor expansion [3,4] is defined by
xðkÞ ¼
M1 X N1 X
¯ expðj2pnk=NÞ, aðm; nÞgðk mNÞ
m¼0 n¼0
(1) pffiffiffiffiffiffiffi where j ¼ 1 and the coefficients a(m, n) can be obtained by
aðm; nÞ ¼
L1 X
following biorthogonality condition [3]: L1 X
¯ gðk þ mNÞ expðj2pnk=NÞgðkÞ ¼ ðL=MNÞdm dn
k¼0
¯ 1; 0pnpN ¯ 1, 0pmpM
ð3Þ
where dk denotes the Kronecker delta. Note that the DGT exists only if the biorthogonality condition is satisfied. 3. Block time-recursive DGT algorithm for the computation of the DGT coefficients In the over-sampling case or the critical sampling ¯ ¼ MN ¯ case, LpMN and L ¼ M N should be satisfied, thus, the over-sampling rate b ¼ MN=L ¼ ¯ ¼ N=N. ¯ Obviously, b ¼ 1 and b41 are M=M corresponding to the critical sampling and the oversampling, respectively. Properly choose the sampling ¯ or N and N, ¯ to make the overparameters M and M, sampling rate b a positive integer. Suppose that ¯ 1, and i ¼ 0; 1; . . . ; b 1, j ¼ 0; 1; . . . ; M xj ¼ ½xðjNÞ; xðjN þ 1Þ; . . . ; xðjN þ N 1ÞT ,
(4)
x ¼ ½xT0 ; xT1 ; . . . ; xTM1 T , ¯
(5)
aij ¼ ½aði þ jb; 0Þ; aði þ jb; 1Þ; ; aði þ jb; N 1ÞT , (6)
¯ expðj2pnk=NÞ, xðkÞgðk mNÞ
(2)
k¼0
where Eq. (2) is termed the DGT [3,4] and Eq. (1) can also be called the inverse DGT. In (1) and (2), ¯ ¯ M and N are the numbers of L ¼ NM ¼ N M, ¯ sampling points in time and frequency domains. M ¯ and N are the frequency and time sampling intervals, respectively. The condition MNXL must be satisfied for a stable reconstruction. The critical ¯M ¯ ¼ NM ¼ L (The numsampling occurs when N ber of coefficients a(m, n) is equal to the number of time samples x(k)), and the over-sampling occurs when MN4L. There may be a loss of information in an under-sampling condition (MNoL). Because the critical sampling leads to bad localization of the analysis window g(k), the over-sampling DGT is often used in signal analysis and processing such as noise reduction for MNR FID signals [9]. Note that the synthesis window g(k) and the analysis window g(k) are all real and periodic in L, and satisfy the
Ai ¼ ½ðai0 ÞT ; ðai1 ÞT ; . . . ; ðaiM1 ÞT T . ¯
(7)
Now let m ¼ i þ jb, (2) can be rewritten as aði þ jb; nÞ ¼
L1 X
¯ jNÞ xðkÞgðk iN
k¼0
expðj2pnk=NÞ
ð8Þ
and the matrix formation of (8) is given by Ai ¼ FCi X
(9)
¯ N-point where F is a block-diagonal matrix with M fast discrete Fourier transform (DFT) matrices on its diagonal 3 2 FðNÞ 7 6 FðNÞ 7 6 7 6 (10) F¼6 7 .. 7 6 . 5 4 FðNÞ
ARTICLE IN PRESS L. Tao, H.K. Kwan / Signal Processing 88 (2008) 407–414
and Ci is a block circulant matrix 3 2 i Ci1 CiM1 C0 ¯ 7 6 i Ci0 CiM2 7 6 CM1 ¯ 7 6 ¯ i , C ¼6 . .. 7 .. .. 6 .. . 7 . . 5 4 Ci1 Ci2 . . . Ci0
and at block time t+1, the DGT coefficients can be expressed as
(11)
with each block Cij being an N N diagonal matrix: 2 6 6 6 Cij ¼ 6 6 4
3
¯ gðjN iNÞ ¯ þ 1Þ gðjN iN ..
. ¯ þ N 1Þ gðjN iN
7 7 7 7. 7 5
(12) The diagonal elements from each block in the first row of Ci comprise one period of the analysis window gðkÞ. Eq. (9) can be written in detail as follows: ai0 ¼ FðNÞ½Ci0 x0 þ Ci1 x1 þ þ CiM1 xM1 ; ¯ ¯ x0 þ Ci0 x1 þ þ CiM2 xM1 ; ai1 ¼ FðNÞ½CiM1 ¯ ¯ ¯ .. . aiM1 ¯
¼
FðNÞ½Ci1 x0
þ
Ci2 x1
þ þ
Ci0 xM1 : ¯
Suppose the block xt is inputted into a block ¯ 1 in series (see delay line from t ¼ 0 to t ¼ M Fig. 1). The time interval of the series inputs of blocks is defined as a block-time interval. We define the DGT at block time t (t ¼ 0; 1; 2; . . .) in order to derive a block time-recursive relation between the updated DGT coefficients and the previous computed DGT coefficients. 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 ¯ ¼ MN. ¯ transform size is L ¼ M N Thus, x¯ tþM1 ai0 ðtÞ ¼ FðNÞ½Ci0 x¯ t þ Ci1 x¯ tþ1 þ þ CiM1 ; ¯ ¯ x¯ t þ Ci0 x¯ tþ1 þ þ CiM2 x¯ tþM1 ai1 ðtÞ ¼ FðNÞ½CiM1 ; ¯ ¯ ¯ .. . aiM1 ðtÞ ¼ FðNÞ½Ci1 x¯ t þ Ci2 x¯ tþ1 þ þ Ci0 x¯ tþM1 ¯ ¯
xt+M-1
···
xt+1
Fig. 1. A block delay line.
xt
ðx¯ tþM¯ x¯ t Þ; ai1 ðt þ 1Þ ¼ ai2 ðtÞ þ FðNÞCiM2 ¯ .. . aiM2 ðt ¯
þ 1Þ ¼
(15)
þ FðNÞCi1 ðx¯ tþM¯ x¯ t Þ; FðNÞCi0 ðx¯ tþM¯ x¯ t Þ:
aiM1 ðtÞ ¯
aiM1 ðt þ 1Þ ¼ ai0 ðtÞ þ ¯
Eq. (15) is called the block time-recursive DGT for the computation of the DGT coefficients. The DGT coefficients Ai ðt þ 1Þ at block time t+1 are given by Ai ðtÞ at block time t with an error correction FðNÞCij ðx¯ tþM¯ x¯ t Þ, i ¼ 0; 1; . . . ; b 1, j ¼ 0; 1; . . . ; ¯ 1. The error correction becomes zero when M ¯ due to the periodicity of x(k), and then the t ¼ M, time-recursive process in (15) gets into stable stage. Note that x¯ t ¼ 0 in the data delay line when ¯ 1. tpM Eq. (15) can also be rewritten as follows:
ðx¯ tþM¯ x¯ t Þ; a~ i1 ðt þ 1Þ ¼ a~ i2 ðtÞ þ CiM2 ¯ .. . a~ iM2 ðt þ 1Þ ¼ a~ iM1 ðtÞ þ Ci1 ðx¯ tþM¯ x¯ t Þ; ¯ ¯
(16)
ðt þ 1Þ ¼ a~ i0 ðtÞ þ Ci0 ðx¯ tþM¯ x¯ t Þ a~ iM1 ¯ ¯ compute and when t ¼ M, aij ¼ FðNÞ~aij
¯ 1. i ¼ 0; 1; . . . ; b 1; j ¼ 0; 1; . . . ; M
(17) Obviously, using (16) and (17) instead of (15) will ¯ 1Þ ðNpoint fast DFTÞ computation. save MðM Note that Ref. [8] only used the formula similar to (15) for the computation of the transform coefficients in the critical sampling case. To summarize, a procedure of the above algorithm is given as follows. Algorithm 1.
(14) xt+M
ðx¯ tþM¯ x¯ t Þ; ai0 ðt þ 1Þ ¼ ai1 ðtÞ þ FðNÞCiM1 ¯
ðx¯ tþM¯ x¯ t Þ; a~ i0 ðt þ 1Þ ¼ a~ i1 ðtÞ þ CiM1 ¯
(13)
xt
409
(a) Let t ( 0, clear all the block units in the block delay line and clear all of a~ ij ði ¼ 0; 1; . . . ; b 1; ¯ 1Þ to be zero. j ¼ 0; 1; . . . ; M (b) Input the signal block xt into the block delay line to let x¯ tþM¯ ( xt .
ARTICLE IN PRESS L. Tao, H.K. Kwan / Signal Processing 88 (2008) 407–414
410
¯ 1, com(c) For i ¼ 0; 1; . . . ; b 1; j ¼ 0; 1; . . . ; M pute
a ti
i at+M
x¯ tþM¯ a~ ij ðt þ 1Þ ¼ a~ ijþ1 ðtÞ þ CiMj1 ¯
i at+M-1
···
i at+1
a ti
Fig. 2. A block delay line.
¯ 1, bring all of Note that a~ iM¯ ¼ a~ i0 . If t ¼ M ¯ 1Þ to Step ¼ 0; 1; . . . ; b 1; j ¼ 0; 1; . . . ; M (d), otherwise let t ( t þ 1 and go to Step (b). ¯ 1, com(d) For i ¼ 0; 1; . . . ; b 1; j ¼ 0; 1; . . . ; M pute a~ ij ði
aij ¼ FðNÞ~aij At present, we obtain the DGT coefficients A ¼ ½Ai b1 i¼0 by Eq. (7).
The diagonal elements from each block in the first row of G i comprise one period of the synthesis window g(k). Eq. (19) can be written in detail as follows:
x0 ¼
b1 X i ¯ i i ¯ ½G i0 FðNÞa 0 þ G 1 FðNÞa1 i¼0 i ¯ FðNÞa þ þ G iM1 , ¯ ¯ M1
4. Block time-recursive DGT algorithm for the reconstruction of the original signal
x1 ¼
b1 X i ¯ i i ¯ FðNÞa ½G iM1 ¯ 0 þ G 0 FðNÞa1 i¼0
¯ ¼ N; bM ¯ ¼ M, Let m ¼ i þ jb, and note that bN Eq. (1) can be rewritten as xðkÞ ¼
¯ b1 M1 N 1 X X X
xM1 ¼ ¯
2pnk exp j N ¼
¯ jNÞ gðk iN
i¼0 j¼0
2pnk exp j . N
i ¯ þ þ G i0 FðNÞa . ¯ M1 N 1 X
aði þ jb; nÞ
n¼0
ð18Þ
(19)
i¼0
¯ where F¯ ¼ F 1 with each diagonal block FðNÞbeing ¯ the inverse matrix of F(N) (i.e., FðNÞ is an N-point fast inverse DFT), G i is also a block circulant matrix: 3 2 i G0 G i1 G iM1 ¯ 7 6 i G i0 G iM2 7 6 G M1 ¯ 7 6 ¯ i (20) G ¼ 6. 7, .. .. .. 7 6 .. . . . 5 4 G i1 G i2 . . . G i0 with each block 2
6 6 6 i Gj ¼ 6 6 4
b1 X i ¯ i i ¯ ½G i1 FðNÞa 0 þ G 2 FðNÞa1 i¼0
The matrix formation of (18) is given by b1 X ¯ i, X¼ G i FA
G ij
.. .
¯ jNÞ aði þ jb; nÞgðk iN
i¼0 j¼0 n¼0
¯ b1 M1 X X
i ¯ FðNÞa þ þ G iM2 , ¯ ¯ M1
As in the block time-recursive algorithm for the computation of the DGT coefficients, we define the DGT at block time t in order to derive a block timerecursive relationship between the updated reconstructed signal and the previous computed reconstructed signal. This is consistent with the new block of the DGT coefficients sliding in and the old block of the DGT coefficients sliding out (see Fig. 2), the block size is N, and the transform size is L ¼ ¯ ¼ MN. ¯ MN Thus,
x0 ðtÞ ¼
¯ FðNÞ¯ aitþM1 þ þ G iM1 , ¯ ¯ x1 ðtÞ ¼
¯ þ 1Þ gðjN iN ..
. ¯ þ N 1Þ gðjN iN
(21)
b1 X ¯ ¯ FðNÞ¯ ait þ G 0 FðNÞ¯ aitþ1 ½G M1 ¯ i¼0
¯ FðNÞ¯ aitþM1 þ þ G M2 , ¯ ¯ .. .
3
7 7 7 7. 7 5
b1 X ¯ ¯ ait þ G i1 FðNÞ¯ aitþ1 ½G i0 FðNÞ¯ i¼0
being an N N diagonal matrix,
¯ gðjN iNÞ
ð22Þ
xM1 ðtÞ ¼ ¯
b1 X ¯ ¯ ait þ G i2 FðNÞ¯ aitþ1 ½G i1 FðNÞ¯ i¼0
¯ aitþM1 þ þ G i0 FðNÞ¯ ¯
ð23Þ
ARTICLE IN PRESS L. Tao, H.K. Kwan / Signal Processing 88 (2008) 407–414
and at block time t+1, the reconstructed signal can be expressed as x0 ðt þ 1Þ ¼ x1 ðtÞ þ
b1 P i¼0
x1 ðt þ 1Þ ¼ x2 ðtÞ þ
b1 P i¼0
¯ FðNÞð¯ aitþM¯ a¯ it Þ; G iM1 ¯
(b) For i ¼ 0; 1; . . . ; b 1, input the DGT coefficient block ait into the ith block delay line to let a¯ itþM¯ ( ait . ¯ 1, compute (c) For j ¼ 0; 1; . . . ; M xj ðt þ 1Þ ¼ xjþ1 ðtÞ þ
¯ FðNÞð¯ aitþM¯ a¯ it Þ; G iM2 ¯
xM1 ðt þ 1Þ ¼ x0 ðtÞ þ ¯
b1 P i¼0
b1 P
¯ aitþM¯ a¯ it Þ; G i1 FðNÞð¯
¯ aitþM¯ a¯ it Þ: G i0 FðNÞð¯
i¼0
Eq. (24) is called the block time-recursive DGT for the reconstruction of the original signal. The reconstructed signal Xðt þ 1Þ at block time t+1 is given by XðtÞ at block time t with an error correction ¯ aitþM¯ a¯ it Þ; G ij FðNÞð¯
¯ FðNÞ¯ aitþM¯ . G iMj1 ¯
i¼0
(24)
b1 X
b1 X
¯ 1, we obtain Note that xM¯ ¼ x0 . If t ¼ M the reconstructed signal X by Eq. (5), otherwise, let t ( t þ 1 and go to Step (b).
.. . xM2 ðt þ 1Þ ¼ xM1 ðtÞ þ ¯ ¯
411
¯ 1. j ¼ 0; 1; . . . ; M
i¼0
To summarize, a procedure of the above algorithm is given as follows. Algorithm 2. (a) Let t ( 0, clear all the block units in the block delay line and clear all of xj ðj ¼ 0; 1; . . . ; ¯ 1; Þ, to be zero. M
5. Unified parallel lattice structures of block timerecursive DGT and its inverse transform 5.1. Unified parallel lattice structures The block time-recursive Eqs. (16) and (17) for the computation of the DGT coefficients can be implemented by b unified parallel lattice arrays as shown in Fig. 3. Each of the parallel lattice arrays ¯ unified lattice modules. It should be consists of M noted that data processing at each stage is by an array size N. The latest sliding-in data block x¯ tþM¯ of size N subtracts the sliding-out data block x¯ t , and then the difference is multiplied by the analysis window block Cij at each module followed by an addition of shifted and previously computed data. A block data of size N needs to be transferred to the nearest-neighbor module in circular order. Finally, with one period of signal sliding in, an N-point fast DFT is computed at the end of each module.
xt i M-1
xt+M
a 0i
i M-2
+
F (N)
a 1i
1i
+
F (N)
i a M-2
F (N)
i a M-1
+
-1 xt+1
F (N) Z-1
xt+M-1 . . .
+
Z-1
Z-1 xt
0i
+ Z-1
Fig. 3. Unified parallel lattice structure of block time-recursive DGT for the computation of the DGT coefficients (i ¼ 0; 1; . . . ; b 1).
ARTICLE IN PRESS L. Tao, H.K. Kwan / Signal Processing 88 (2008) 407–414
412 a ti
i=0, 1, · · ·, −1 i
i
G M−1
at+M
+
+
x0
i at+M-1
+
F (N)
i=0, 1, · · ·, −1 Z-1 i G1
+
+
+
+
xM-2
-1 i at+1
i=0, 1, · · ·, −1 Z-1
a ti
i
G0
xM-1 Z-1
Fig. 4. Unified parallel lattice structure of block time-recursive DGT for the reconstruction of the original signal from the DGT coefficients.
The block time-recursive Eq. (24) for the reconstruction of the original signal can also be implemented by b unified parallel lattice arrays. Each ¯ unified lattice modules, which of them consists of M is illustrated in Fig. 4. Data processing at each stage is by an array size N. The latest sliding-in data block a¯ itþM¯ of size N subtracts the sliding-out data block a¯ it , and then an N-point fast inverse DFT of the difference is applied. The resulting data is multiplied by the synthesis window block G ij at each module followed by b1 additions of the resulting block data from the other b1 parallel lattice arrays, and one addition of shifted and previously computed reconstructed signal. A block data of size N needs to be transferred to the nearest-neighbor (output) module in circular order. Obviously, the oversampling case shown in Figs. 3 and 4 will become the critical sampling case when ¯ and N ¼ N ¯ (i.e., b ¼ 1). M¼M 5.2. Computation time and complexity analysis Since the total computational complexity in a parallel algorithm is equally shared by each unified parallel lattice module, the computation time of the parallel algorithm depends on the computational complexity of each unified parallel lattice module, while the computation time of a series algorithm depends on its total computational complexity. In Fig. 3, each parallel module carries N-point multiplications and N-point additions at each block time, and an N-point fast DFT (with computational complexity 0:5N log2 N) at the last block time.
Table 1 Comparison of computational complexity related to total computational time References
Computational complexity
Applicability
[4] [6] [7] [8] Proposed algorithms
ML þ 0:5ML log2 N ML+0.5L log2 N 0.5L log2(4L)+L log2 M L þ 0:5L log2 N L þ 0:5N log2 N L þ 0:5L log2 N
CS, CS, CS, CS, CS, CS,
OS, Coef Coef, Reco Coef, Reco Coef OS, Coef OS, Reco
Therefore, to obtain a block of N transformed data, the computational complexity is of the order of L þ 0:5N log2 N.
(25)
In Fig. 4, at each block time, every parallel module of each parallel array carries an N-point fast inverse DFT, N-point multiplications followed by b Npoint additions. Therefore, to obtain a block of N reconstructed data, the computational complexity is of the order of L þ 0:5L log2 N.
(26)
Table 1 gives a comparison between the proposed algorithms and the main existing DGT algorithms. In the table, the symbols CS, OS, Coef and Reco indicate four cases of algorithm applicability which are respectively the critical sampling case, the oversampling case, the computation of the DGT coefficients, and the reconstruction of the original signal from the coefficients (i.e., the inverse DGT).
ARTICLE IN PRESS L. Tao, H.K. Kwan / Signal Processing 88 (2008) 407–414
We can see from the table that the computational complexity in (25) or (26) related to the total computation time of the proposed parallel algorithms is far less than that of the existing DGT algorithms. Thus, the proposed algorithms are more efficient and faster in computing the DGT coefficients and reconstructing the original signal from the coefficients.
413
DGT coefficients |a(m,n)|
Fourier spectrum Fx (f) 8 kHz
127
4 kHz
0
0 kHz 0
127
0
70
A speech signal x (k) 1
6. Computer simulations
0
In this section, we present some numerical simulation results, which are computed by the proposed algorithms. In the oversampling case ¯ ¼M ¯ ¼ 16, b ¼ 8), the top plot (N ¼ M ¼ 128, N of Fig. 5 depicts a Gaussian synthesis window, qffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffi gðkÞ ¼ 2=40 expfp½ðk 1023:5Þ=402 g, with length L ¼ 2048, corresponding to its analysis window g (k) in the bottom plot of Fig. 5, which is computed by Eq. (3). g(k) and g(k) are used to compute the DGT coefficients of a speech signal x(k) and to reconstruct the original speech signal by the DGT coefficients. The speech signal x(k) is shown in Fig. 6 with length L ¼ 2048, i.e., 256 ms of the word ‘‘yes’’ spoken by a man. x(k) is partitioned by Eq. (5) into 16 nonoverlapping block vectors of ¯ 1Þ. Correspondlength 128 ðxj ; j ¼ 0; 1; . . . ; M ingly, the DGT coefficient matrix a(m, n) of the speech signal is also partitioned by Eq. (7) into ¯ nonoverlapping block vectors 8 16 (i.e., b M) of length 128 ðaij ; i ¼ 0; 1; . . . ; b 1; j ¼ 0; 1; . . . ; ¯ 1Þ. By using the proposed algorithm for the M computation of the DGT coefficients, Fig. 6 displays the computed DGT coefficients (magnitude) of x(k). In order to compare with the DGT coefficients, Fig. 6 also shows one period of the
g (k)
0.2 0.1 0 0
512
1024 k
1536
2048
0
512
1024
1536
2048
(k)
0.04 0.02 0 k Fig. 5. Gaussian window g(k) and its analysis window g (k). ¯ ¼ 16 (oversampling). L ¼ 2048, N ¼ 128, N
-1 0 ms
256 ms
Fig. 6. DGT coefficients jaðm; nÞj and Fourier spectrum Fx(f) of a speech signal x(k).
Table 2 Comparison of total number of multiplications related to total computational time References
Total number of multiplications (applicability)
[4] [6] [7] [8] Proposed algorithms
786432 (CS, Coef), 1179648 (OS, Coef) 266240 (CS, Coef, Reco) 27648 (CS, Coef, Reco) 6144 (CS, Coef) 2080 (CS, Coef), 2496(OS, Coef) 6144 (CS, Reco), 9216(OS, Reco)
Fourier spectrum Fx(f) of x(k), which reveals four frequency tones in the range from 0 to 4 kHz (half the sampling frequency 8 kHz). From the spectrum alone, we can not tell how those frequencies change over time; however, from the DGT coefficients, not only can we see how frequencies (n) change over time (m), but we can also see the intensity of the frequencies as shown by the relative gray levels of the image. The proposed algorithm for the reconstruction of the original signal is also simulated. The signal reconstruction error (MSE) is around 10–15, which is virtually error-free reconstruction. A numeric comparison on the total number of multiplications related to the total computational time between the proposed algorithms and the main existing DGT algorithms is given in Table 2, which further indicates the advantages of the proposed algorithms. In Table 2, the parameters L ¼ 2048, M ¼ N ¼ 128 for oversampling and L ¼ 2048, M ¼ 128, N ¼ 16 for critical sampling are used. For example, ‘‘6144 (CS, Coef)’’ for Ref. [8] means that, for computing the DGT coefficients in the
ARTICLE IN PRESS 414
L. Tao, H.K. Kwan / Signal Processing 88 (2008) 407–414
critical sampling case (L ¼ 2048, M ¼ 128, N ¼ 16), the total number of multiplications of the algorithm in Ref. [8] is 6144, which is calculated using the formula of Ref. [8, Table 1].
Talent Development Foundation under Grant no. 2005Z029, the Anhui University Program for Creative Teams, and the MOE’s Key Laboratory of Intelligent Computing and Signal Processing in Anhui University.
7. Conclusions In this paper, we developed two block timerecursive DGT algorithms in both the critical sampling case and the oversampling case, not only for the computation of the DGT coefficients but also for the reconstruction of the original signal from the coefficients. The two algorithms were implemented respectively by a unified parallel lattice structure. Computational complexity analysis and comparison showed that the proposed algorithms are more efficient and faster than the existing DGT algorithms in computing the DGT coefficients and reconstructing the original signal from the coefficients. Moreover, the proposed algorithms can also be generalized to the DHT (discrete Hartley transform)-based real-valued DGT, which was presented in our previous work [10]. Acknowledgments The authors would like to thank the anonymous reviewers for their constructive comments. The first author would like to acknowledge the supports of the National Natural Science Foundation of China under Grant no. 60572128, the Anhui Provincial
References [1] D. Gabor, Theory of communication, J. Inst. Electr. Eng. 93 (3) (1946) 429–457. [2] M. Bastiaans, Gabor’s expansion of a signal into Gaussian elementary signals, Opt. Eng. 20 (4) (1981) 594–598. [3] J. Wexler, S. Raz, Discrete Gabor expansions, Signal Processing 21 (3) (1990) 207–220. [4] S. Qian, D. Chen, Discrete Gabor transforms, IEEE Trans. Signal Process. 21 (7) (1993) 2429–2438. [5] S. Qian, D. Chen, Joint time–frequency analysis, IEEE Signal Process. Mag. 6 (2) (1999) 52–67. [6] L. Wang, C. Chen, W. Lin, An efficient algorithm to compute the complete set of discrete Gabor coefficients, IEEE Trans. Image Process. 3 (1) (1994) 87–92. [7] S. Qiu, F. Zhou, P.E. Crandall, Discrete Gabor transforms with complexity O(Nz log N), Signal Processing 77 (2) (1999) 159–170. [8] C. Lu, S. Joshi, J.M. Morris, Parallel lattice structure of block time-recursive generalized Gabor transforms, Signal Processing 57 (2) (1997) 195–203. [9] B. Friedlander, A. Zeira, Oversampled Gabor representation for transient signals, IEEE Trans. Signal Process. 43 (9) (1995) 2088–2094. [10] L. Tao, H.K. Kwan, Block time-recursive real-valued discrete Gabor transform implemented by unified parallel lattice structures, IEICE Trans. Inform. Systems E88-D (7) (July 2005) 1472–1478.