IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 41. NO. 3, MARCH 1993 ~ 1453
1
.o
0.8
0.6
0.4
0.2
0
<
)Estimate - Originall %Original Algorithm / Direct Algorithm d4
%
AR P = Algorithm 30 cap Length Record Length0.2
0.4
0.6
Fig. 4. Errors in spectral estimates for damped-sine-wave data set pendent of gap length over the range tested and were the best for gap lengths greater than 2 0 % .
REFERENCES
[l] S . L. Marple, Jr., Digital Specrral Analysis with Applica- tions. Englewood Cliffs, NJ: Prentice-Hall, 1987.
121 S . Kay, Modern Spectral Estimation: Theory and Application. En- glewood Cliffs, NJ: Prentice-Hall, 1988.
[3] A. K. Jain and S . Ranganath, “Extrapolation algorithms for discrete signals with application in spectral estimation,” IEEE Trans. Acoust. , Speech, Signal Processing, vol. 29, pp. 830-845, Aug. 1981. [4] A. V . Steinhardt, R. K . Goodrich. and R. A. Roberts, “Spectral es-
timation via minimum energy correlation extension,” IEEE Trans. Acoust., Speech, Signal Processing, vol. 33, pp. 1509-1515, Dec. 1985.
[5] N. Rozario and A. Papoulis, “Spectral estimation from nonconse- cutive data,” IEEE Trans. Inform. Theory, vol. 3 3 , pp. 889-894, Nov. 1987.
[6] S. D. Cabrera and T. W . Parks, “Extrapolation and spectral esti- mation with iterative weighted norm modification,” IEEE Trans. Sig- n a l f r o c e s s i n g , vol. 39, pp. 842-851, Apr. 1991.
[7] R. W . Schafer and L. R. Rabiner, “A digital signal processing ap- proach to interpolation,” Proc. IEEE, vol. 61, pp. 692-702, June
1973.
[8] A. Papoulis, “A new algorithm in spectral analysis and band-limited extrapolation,” IEEE Trans. Circuits Syst., vol. 22, pp. 735-742, Sept. 1975.
[9] J. A. Cadzow, “An extrapolation procedure for band-limited sig- nals,” IEEE Trans. Acoust., Speech, Signal Processing. vol. 27, pp. 4-12, Feb. 1979.
[lo] R . J . Marks, 11, “Restoring lost samples from an oversampled band- limited signal,” IEEE Trans. Acousr., Speech, Signal Processing,
vol. 31, pp. 752-755, June 1983.
[ l l ] W . Xu and C. Chamzas, “On the extrapolation of band-limited func- tions with energy constraints,’’ IEEE Trans. Acoust. , Speech, Signal Processing, vol. 31, pp. 1222-1234, Oct. 1983.
[12] J. L. C. Sanz and T. S . Huang, “A unified approach to noniterative linear signal restoration,” IEEE Trans. Acoust. , Speech, Signal Pro- cessing, vol. 32, pp. 403-409, Apr. 1984.
[13] R. J. Marks, 11, and D. Radbel, “Error of linear estimation of lost samples in an oversampled band-limited signal, ” lEEE Trans.
Acoust., Speech, Signal Processing. vol. 32, pp. 648-654, June 1984. [14] P. J . Thompson and M. T. Manry, “A frequency-domain weighting function approach to extrapolation,” IEEE Trans. Acousi., Speech, Signal Processing, vol. 38, pp. 1395-1402, Aug. 1990.
[15] S . T. Fathima and B . Yegnanarayana, “A maximum entropy ap- proach to interpolation,” Signal Processing, vol. 21, pp. 17-24,
1990.
1161 R. W. Schafer, R. M. Mersereau, and M. A. Richards, “Constrained iterative restoration algorithms,” Proc. IEEE, vol. 69, pp. 432-450, Apr. 1981.
[I71 B . J. Sullivan and B . Liu, “On the use of singular-value decompo- sition and decimation in discrete-time band-limited signal extrapola- ’ tion,” IEEE Trans. Acousr., Speech, Signal Processing, vol. 32, pp.
1201-1212, Dec. 1984.
[I81 S . Singh and S . N . Tandon, “An iterative restoration technique,”
Signal Processing, vol. 11, pp. 1-1 1, 1986.
Fast Orthogonalization Algorithm and Parallel
Architecture for AR Spectral Estimation Based
on Forward-Backward Linear Prediction
K. J. Ray Liu and S. F . HsiehAbstract-The truncated QR methods have been shown to be com- parable to the SVD-based methods for the sinusoidal frequency esti- mation based on the forward-backward linear prediction (FBLP) model. However, without exploiting the special structure of the FBLP matrix, the QR decomposition (QRD) of the FBLP matrix has the com- putational complexity on the order of 2(6m - n ) n 2 / 3
+
O ( n 2 ) for a 2m x n FBLP matrix. Here we propose a fast algorithm to perform the QRD of the FBLP matrix by exploiting its special Toeplitz-Hankel form. The computational complexity is then reduced to 10n2+
4mn+
O ( n ) . The fast algorithm can also be easily implemented onto a linear systolic array. The number of time steps required is further reduced to 2m+
5n - 4 by using the parallel implementation.I. I N T R O D U C T I O N
High-resolution spectral estimation i s an important subject I n
many applications of modem signal processing. The fundamental problem in applying various high-resolution spectral estimation al- gorithms is the computational complexity. In the pioneering paper of Tufts and Kumaresan [ l ] , a SVD-based method for solving the forward-backward linear prediction (FBLP) least squares (LS) problem was used to resolve the frequencies of closely spaced si- nusoids from a limited amount of data samples. By imposing an excessive order in the FBLP model and then truncating small sin- gular values to zero, this truncated SVD method yields a low sig- nal-to-noise ratio (SNR) threshold and greatly suppresses spurious frequencies. However, the massive computation required by SVD makes it unsuitable for real-time superresolution applications.
Recently, the truncated QR methods [4] have been shown to be comparable to the SVD-based methods in various situations. It is
very effective for the sinusoidal frequency estimation based on the FBLP model. However, without considering the special structure
of the FBLP matrix, the QR decomposition (QRD) o f the FBLP matrix still has the computational complexity on the order of 0 ( n 3 ) .
Manuscript received April 21, 1991; revised May 12, 1992.
K . J. R. Liu is with the Department of Electrical Engineering, Systems Research Center, University of Maryland, College Park, MD 20742.
S . F. Hsieh is with the Department of Communication Engineering, Na- tional Chiao Tung University, Hsinchu, Taiwan 30039, Republic of China.
IEEE Log Number 920601 1 . 1053-587X/93$03.00 0 1993 IEEE
-
14.54
Seeking fast algorithms for specially structured matrices has cap- tured lots of attention recently, especially the Toeplitz-structured matrices are used in many signal processing applications [ 2 ] , [ 3 ] , [ 7 ] - [ l o ] . However, exploiting the special structure of the FBLP
matrix for fast algorithm implementation has not yet been consid- ered so far. Here we propose a fast algorithm to perform the QRD of the FBLP matrix. The computational cost of the truncated .QR
methods can be further reduced from O ( n 3 ) to O ( n 2 ) which makes it more attractive than the SVD-based methods. Without exploiting the special structure of the FBLP matrix, the straightforward Q R D of the FBLP matrix has the computational complexity on the order of 2(6m - n ) n 2 / 3
+
O ( n 2 ) for a 2m X n FBLP matrix. The pro- posed fast algorithm reduces it to the order of 10n2 + 4mn+
O ( n ) .We will also show that the proposed fast algorithm is amendable to parallel processing. A fully pipelined linear systolic array based on the multiphase operations is used to implement the fast algo- rithm parallelly. The required time steps is further reduced to 2m
+
5 n - 4 .The idea of FBLP was originated by Burg 1.51 for the lattice pre- dictors. To improve the performance, Tufts and Kumaresan [ l ] de-
veloped a modified FBLP method which is very effective for esti- mating closely spaced frequencies. The FBLP method is to minimize the sum of the FBLP errors energy,
N
min
C
[Jef(i)l2+
l e b ( i ) 1 2 ~ (1)i = M + I
where e, and eb are the forward, and backward residuals, respec-
tively. The FBLP method is to solve the LS problem [ 1 3 ] :
I
IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 41. NO. 3, MARCH 1993
min ) ) A V - b J J 2 w with and
u(M
+
1)u(M)
. . .
T = [ u ( N'
-
1) u ( N - 2 ). . .
u ( N )u(N
- 1). . .
where J is an exchange matrix. The matrix of the form as given in
( 3 ) is called the Toeplitz-Hankel matrix. In fact, the augmented
matrix of the FBLP problem is of the Toeplitz-Hankel form with a special property, i.e., H = TJ. This special property can be used for developing a fast algorithm that will be considered in the fol- lowing sections.
11. EXPLOITING T H E TOEPLITZ-HANKEL S T R U C T U R E In using the truncated QR method for the high-resolution AR
spectral estimation, the key computational issue is to solve the FBLP LS problem based on the QR decomposition (QRD). Without
considering the special structure, a conventional Q R D requires
= 4 ( N - M ) M '
+
O ( M 2 ) multiplications to obtain the upper tri- angular matrix R . This is on the order of O ( N M 2 ) since usually N>>
M . Thus, a reasonable approach is to find a fast algorithm for the FBLP LS problem by exploiting its special Toeplitz-Hankel structure. This problem has not been considered so far, though the LS problem with Toeplitz structure has been studied extensively121, (31. [61-t101.
The Toeplitz part of the Toeplitz-Hankel matrix can be parti- tioned as where L u ( N - 1 ) u ( N - 2 ) .
. .
u ( N1
M ) ] x T = [U(M),. . .
, u ( 2 ) . u ( l ) ] Y T = [ u ( M+
2 ) ,.
*.
,
u ( N - l ) , u ( N ) ] [ ~ ( l ) , u ( 2 ) , * *.
, u ( N - M - l ) ] U T = U T = [ u ( N ) , u ( N - l ) ,. . .
, u ( N - M+
l ) ]and the Hankel part of the Toeplitz-Hankel matrix can be parti- tioned as H X B T
u ( M +
1 ) [ I ( N - M ) v B T ] = [ H y H = T J = where H =TJ
vBr = [ u ( N - M+ l ) ,
. . .
, u ( N - I), u ( N ) ] vTJ u ( N - M+
1 )u(N
- M ) (4) and X B T = [ u ( l ) , u ( 2 ) ,. . .
, u ( M ) ] = x T J .Here B denotes the backward arrangement of a vector. can be partitioned as follows:
Now, from the above partitions, the Toeplitz-Hankel matrix K
( 7 )
and
+ u'(N -
M )+
u'(M+
1 )+
y'y u'H+
u ( N - M ) v B T+
u(M
+
l ) x T+
y r THTu
+
v B u ( N - M )+
xu(M+
1 )+
T3
H T H+
FrT
+
vBvBT+
x x T__ ._
I
IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 41, NO. 3, MARCH 1993 1455
Also, the matrix K can be partitioned as R:R? = R: R~
+
vB vB1R:R~ = R:R? - U V ~
(9) R ~ R , = R ~ R , - r,rT (19)
I where R I , Rh, and R, , i = 1, 2 , 3, 4 are all M X M upper triangular matrices, and all of the vectors involved are M-dimensional. In and with this partition, we have
(10)
x T x B T
+
H T H+
F7F
+
V V ~ x B u ( M+
1)+
H r y+
F7u+
VU(N
- M )u ( M
+
1 ) ~ ~ ’+
y 7 H+
uTT+
u ( N - M)v‘ u’(M+
1)+
y 7 y+
U ’ U+
u L ( N - M )K ~ K =
Let the QRD of the matrix
K
be K = QR, where R E IS an upper triangular matrix and it can also be par-@ M + l ) X ( M + I ) ’
titioned as follows:
where Rh E
a M x M
is the principal bottom submatrix of R , RI E(RM
is the principal top submatrix of R, andrT = tr1.2, r1.5
.
* . 7 T I . M +I 1
rT = [ T I . M + I , T ? . M + I .. . .
3 T M . M +I1
order to start the above recursions, the first row of R, (or R ) must be available. In general, there is no shortcut for obtaining this row and it can be done by a sequence of Givens rotations on the matrix
K to zero out the first column of K , except its leading element on
the diagonal. By denoting # a s a “don’t care” element or vector, the fast algorithm is summarized in Table I.
As we can see, for the initialization (obtaining the first row of
R ) , the computational cost is = 4 ( N - M ) M multiplications (since only half of the rotation needed to be done). Following this, the recursions in the main iterations are then started. As there are five rotation-like upidowndatings, the computational cost is
Note that both Rb and RI are upper triangular matrices. Since the
matrix Q is orthogonal, we have 5 x ( 4 M + 4 ( M - 1 ) +
. . .
+ 4 . I )K ~ K = R ~ R (12)
and
(for multiplication). Therefore, the total computational complexity is = 10M’
+
4 ( N - M ) M (for multiplication) for a 2 ( N - M ) XM Toeplitz-Hankel matrix. As mentioned before, without consid-
ering the special structure, by using the conventional Q R D , the computational complexity is of = 4 N M 2
+
O ( M 2 ) . Obviously, the proposed fast algorithm has an improvement of an order of mag- nitude. In general, for the Q R D of a 2m X n Toeplitz-Hankel ma- trix, the fast algorithm needs Ion2+
4mn+
O ( n ) multiplications,while a conventional implementation needs 2(6m - n ) n 2 / 3 +
O(n’), where m
>>
n.It the least squares weight vector is of interested, a backward substitution can then be used for computing the weight vector. For the truncated QR method, a truncation of the noise subspace is nec- essary before computing the weight vector [4].
(13)
d.1
r1.1.T d R , rTr2+
r L + I . M + I Define-
A
K =[-+-I
(I4) then we haveK T K
H T H +t T 7 .
(15)From the lower right submatrices of (8) and (13), we obtain
RLRb
+
r l r r = KTrt+
x x T+
vB v”.Also, from the upper left submatrices of (IO) and (13), we have (16)
RTR,
= K T K+
+
x B x B T . (17)Substituting (17) to (16), we obtain the relation between Rb and RI
as given by
RIRb = RYR,
+
x x T - x B r B T+
vBvBT - vuT - r , r : . (18) 111. T H E FAST ALGORITHMIt is clear how to perform the updating and downdating of the Cholesky factors [2]. As we can see, in [18], there are two rank-I updatings and three rank-1 downdatings. Let us split (18) into a sequence of five upldowndating equations given by
R:R, = R:R,
+
XXIR ~ R , = R ~ R , - x B x B 1
IV. PAKALLFL IMPLEMENTATION
The fast algorithm obtained in the previous section not only re- duces the computational complexity, but is also amenable for par- allel implementation. From the fact that only the first row of the upper triangular matrix R has to be obtained first, a linear array of M
+
1 processing cells, as shown in Fig. I , can be used to rotate the matrix K such that the first column can be zeroed out and when the initialization phase is finished, the first row of the matrix R is kept in the linear array. Fig. 2 shows the initialization to obtain the first row of R. The operations of the processing cells are given in Table 11. The data matrix is arranged in a skewed manner for the systolic array implementation. The idea is similar to the triangular array for the Q R D proposed by Gentleman and Kung [ I 11. The difference is that their scheme is a general one without considering any special structure of the data matrix. Accordingly, a full trian- gular array is needed.I
1456 IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 41, NO. 3, MARCH 1993
TABLE I
SUMMARY OF T H E FAST ALGORITHM (Initialization)
X7'(ll = [U(M), u(M - l ) , .
. .
, u ( l ) ]xBr"" = [ u ( l ) , u(2), . . . , u(M)]
U""' = [ u ( N ) , u ( N - I ) , . . ., u(N - M
+
I)]u ~ ~ " ' ' = [ u ( N - M
+
I ) , u(N - M+
2), . . . , u ( N ) ]kli =
kt F o r i = I t o 2 ( N - M ) - 1 , End For; [rl.l. rrl = rr;.#I7
= k L , N - M i ) r;'"' = rT, (Main Iterations) Fori = 1 t o M - 1, (Phase 1 )[
O,;:l,l] = I (Phase 2) (Phase 3) (Phase 4) (Phase 5)U:'+ I = r f r excluding the last one, End For.
X X I
Y Fig. 1. The linear systolic array and its processing cells
K
f
Fig. 2 . The initialization.
Due to the consideration of the special Toeplitz-Hankel struc- ture, once the first row of the matrix R is available, the subsequent rows of R can be generated one by one by the main iterations given in the fast algorithm. T o start the main iterations, r',' is needed. Fortunately, it is the first M elements of the first row of R that are stored in the first M processing cells. The main iterations are now started with inputs x r"", xRT"", x""", uT"", and r ; ' " ' , and the outputs , and r:"', respectively, as illustrated in Fig. 3. The outputs have one less dimension than their inputs do. The vector rp' can now be obtained on the linear array.
As given in the main iterations, there are five different phases. The operations of the processing cells for different phases are given in Table 11. Based on the multiphase concept proposed in [12], the outputs are fedback to the input ports for another iteration of dif- ferent phases. Note that the outputs are obtained from PE2 to PEM.
The feedback is, however, directed to the processors from PE1 t o
PE(M - I ) . Since r r take the first M - 1 elements of rir, it oc- cupies the first M - 1 processing cells. The second iteration is started once the fedback data are available. It is fully pipelined without any intermediate data arrangement and interrupt. The it- erations are then continued until all the rows of R are obtained. The overall data arrangement is given in Fig. 4. Since only a left shift is performed in the feedback loop, it is obvious that a linear array as shown in Fig. 5 can be used without the need of feedback.
The number of time steps required for this linear array imple- mentation is now being further reduced to 2 ( N
-
M )+
( 5 ( M -1)
+
1) = 2N+
3M - 4 (or 2m+
Sn - 4 for a 2m X n Toeplitz- Hankel matrix) which is linearly proportional to either M or N (mor n ) .
If the LS weight vector is of interest, another phase for the back- ward substitution can be started easily since all the data are now available in the linear array. The details of the operations of the backward substitution using a linear array can be found in [ I 11. are To' x B T ' t t Xb'T'll
,
V . CONCLUSIONS
In this correspondence, we propose a fast algorithm for the Q R D of a Toeplitz-Hankel matrix. The computational complexity for the QRD of a 2m X n Toeplitz-Hankel matrix is 10n2
+
4mn+
O ( n )multiplications, which has an order of magnitude improvement over conventional algorithms. This algorithm can also be implemented onto a fully pipelined multiphase linear systolic array. The number of time steps required is further reduced to 2m
+ 5 - 4 for the
parallel implementation. An interesting point for the QRD of the specially structured matrices such as Toeplitz and Toeplitz-Hankel forms is that there is no need to store all the generated rows of the
IEEE TRANSACTIONS ON SIGNAL PROCESSING. VOL. 41. NO. 3 . MARCH 1993 1457
TABLE I1
THE O P E R A T l O l S O F THF PROCESSING C E L L S IN DIFFERFNT PHASES
Initialization Phases 1 and 3 Phases 2 , 4 , and 5 PE 1 d = m d = . F s d = ./,= c' = r / d . s = . r / d r = (I r = d r = d c = r / d , s = . r / d c = d / r , s = x / r PEi. 1 C- i 5 M
+
1 r = ('r+
sx r = cr+
sx v = - s r+
ex J = - s r+
ex r = ( r - sx)/c v = -.sr+
cx1
c
Fig. 3 . The first iteration with five different phases
-M
1458
I
IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 41, NO. 3, MARCH 1993
Fig. 5. The linear array for parallel implementations without feedback upper triangular matrix R . As long as the first row of R is known, all the subsequent row can be generated recursively, and this is also the basic principle of the proposed fast algorithm.
REFERENCES I l l I21 131 141 151 161 171 181 191 I101 1111 [I21
D. W. Tufts and R. Kumaresan, “Estimation of frequencies of mul- tiple sinusoids: Making linear prediction perform like maximum likelihood,” Proc. IEEE, vol. 70, no. 9, pp. 975-989, Sept. 1982. A. W. Bojanczyk. R. P. Brent, and F. R. de Hoog, “ Q R factorization of Toeplitz matrices,” Numer. M a t h . , vol. 49, pp. 81-94, 1986. D. R. Sweet, “Fast Toeplitz orthogonalization,” Numer. Math., vol. 43, pp. 1-21, 1984.
S . F. Hsieh, K. J . R . Liu, and K. Yao, “Estimation of multiple sin- usoidal frequencies using truncated least squares methods,” f E E E
Trans. Signal Processing, vol. 41, no. 2, Feb. 1993.
J. P. Burg, “Maximum entropy spectral analysis.” presented at the 37th Annu. Int. Meeting, Soc. Explor. Geophys.. Oklahoma City, OK, 1967.
C. P. Rialan and L. L. Scharf, “Fast algorithms for computing Q R
and Cholesky factors of Toeplitz operators,” fEEE Trans. Acou.st., Speech, Signal Processing, vol. 36, pp. 1740-1748, Nov. 1988. T. Kailath, S . - Y . Kung, and M. Mod, “Displacement ranks of matrices and linear equations,” J . Math. Anal. A p p l . , vol. 68, pp. S.-Y. Kung and Y . H. Hu, “A highly concurrent algorithm and pipe- lined architecture for solving Toeplitz systems,” IEEE Trans. Acoust., Speech, Signal Processing, vol. 3 1 , Feb. 1983.
A. W. Bojanczyk, R. P. Brent, and F. R. de Hoog. “Linearly con- nected arrays for Toeplitz least squares problems,” J . Parallel Dis- tributed Computing, vol. 9, pp. 261-270, 1990.
J . Chun, T. Kailath. and H. Lev-Ari, “Fast parallel algorithms for
Q R and triangular factorization,” SIAM J . Sci. Stat. Cornput.. vol.
8. pp. 899-9 12, Nov. 1987.
W. M. Gentleman and H. T . Kung, “Matrix triangularization by sys- tolic array,’’ Proc. SPIE, vol. 298, pp. 298, 1981.
K. J. R . Liu and K . Yao, “Multiphase systolic algorithms for spectral decomposition,” IEEE Trans. Signul Processing, vol. 40, no. 1. pp.
190-201, Jan. 1992. 395-407, 1979.
1131 K. J . R . Liu and S. F. Hsieh, “Fast Toeplitz-Hankel orthogonali- zation for high-resolution spectral estimation.” i n Proc. IEEE ICASSP, May 1991, pp. 2249-2252.
The Relationship Between Instantaneous Frequency
and Time-Frequency Representations
Brian C . Lovell, Robert C . Williamson, and Boualem BoashashAbstract-We give the relationship between instantaneous frequency estimation via the derivative of the phase of the analytic signal and the first moment of general time-frequency representations from Cohen’s
Manuscript received March, 2 , 1990; revised January 25, 1992. This work was supported by the Australian Research Council.
B. C. Lovell and B. Boashash are with the Department of Electrical En- gineering, University of Queensland, St. Lucia, Brisbane, Australia 4067. R. C. Williamson is with the Department of Systems Engineering and Interdisciplinary Program in Engineering, Australian National University, Canberra, Australia 2601.
IEEE Log Number 9206002.
. ,
1053-587X/93$03.00 @ 1993 IEEE
class in both the continuous and discrete-time domains. Many re- searchers have applied the standard linear definition of first moment to discrete-time time-frequency representations although this leads to biased instantaneous frequency estimators with high variance; we show that periodic (circular) definitions of moments must be used to account for the periodization of the frequency variable due to sampling.
I. INTRODUCTION
Several authors [3], [9] have investigated the possibility of using the first moments of time-frequency representations with respect to the frequency variable as estimators of instantaneous frequency. This correspondence derives the relationships between instanta- neous frequency and the first moments of the general class of time- frequency representations for both continuous and discrete-time signals.
11. CONTINUOUS-TIME ESTIMATION
Consider a frequency modulated sinusoidal signal of the form
x ( t ) = a, ( t ) cos
4 ( t ) ,
where a, represents the envelope function and 4 is the cumulative phase of the signal. W e define the instan- taneous frequency of this signal by1
X((t)
=2a dt
If x is sufficiently narrow band, a good estimate of the cumulative phase reduced modulo 27r may be obtained from the phase of the
analytic signal defined as follows:
Definition 1 : Analytic Signal: The analytic signal z associated with the real signal x is defined by
z
= A [x], where A [x] = x+
j H ( x ] is the operator which forms the analytic signal and H [ ] isthe Hilbert transform defined by
We use the derivative of the phase of the analytic signal to define the following instantaneous frequency estimator.
Definition 2 : Analytic Derivative Estimator: Let z = A [ x ] , where x is a real signal. Then the instantaneous frequency of x at time t is estimated by
where (( )),= denotes reduction modulo 27r and ((dldt)),, denotes the appropriate differentiation of a quantity which is reduced mod-
ulo 21r a s shown above.
The spectrogram (or magnitude-squared short-time Fourier transform) and time-frequency distributions such as the Wigner- Ville, Bom-Jordan-Cohen, Margenau-Hill-Rihaczek, and Choi- Williams exponential distributions can all be examined within the framework of Cohen’s general class of time-frequency represen- tations [ 2 ] .
Definition 3: Cohen ’s Class of Time-Frequency Representations for Analytic Signals: Each member of this class of bilinear rep-
