Synthesis of spectral densities using finite automata

1993

IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 41, NO. 5, MAY 1993

ACKNOWLEDGMENT

I. INTRODUCTION

The authors greatly appreciate the helpful suggestions from the anonymous referees. The first author thanks A. Shatz for his help in providing some of the experimental results.

REFERENCES [l] L. R. Rabiner, “A tutorial on hidden Markov models and selected applications in speech recognition,” Proc. IEEE, vol. 77, pp. 257286, Feb. 1989. [2] B. S . Atal, “Effectiveness of linear prediction characteristics of the speech wave for automatic speaker identification and verification,” J. Acoust. Soc. Amer., vol. 55, no. 6, pp. 1304-1312, June 1974. [3] B. S. Atal, “Automatic recognition of speakers from their voices,” Proc. IEEE, vol. 64, pp. 460-475, Apr. 1976. [4] L. R. Rabiner and J. G.Wilpon, “Some performance benchmarks for isolated world speech recognition systems,” Comput. Speech and Language, vol. 2, pp. 343-357, 1987. [5] L. R. Rabiner, B.-H. Juang, S. E. Levinson, and M. M. Sondhi, “Some properties of continuous hidden Markov model representations,” AT&T Tech. J., vol. 64, no. 6, pp. 1251-1270, July-Aug. 1985. [6] L. R. Rabiner, B.-H. Juang, S. E. Levinson, and M. M. Sondhi, ‘ ‘Rgcognitionof isolated digits using hidden Markov models and continuous mixture densities,” AT&TTech. J . , vol. 64, no. 6, pp. 12111234, July-Aug. 1985. [7] T. W. Anderson, An Introduction to the Statistical Analysis of Data. New York: Wiley, 1971. [8] S. Kay and J. Makoul, “On the statistics of the estimated reflection coefficients of an autoregressive process,” .IEEE Trans. Acoust. , Speech, Signal Processing, vol. ASSP-31, no. 6, pp, 1447-1455, Dec. 1983. [9] H. B. Mann and A. Wald, “On the statistical treatment of stochastic difference equations,” Econometrica, vol. 11, pp. 173-220, July 1943. [lo] D. R. Brillinger, Time Series Data Analysis and Theory. Holden Day, 1981. 1111 Y. Chen, “Cepstral domain stress compensation for robust speech recognition,” in Proc. IEEE Int. Con$ Acoust., Speech, Signal Processing, Apr. 1987, p~).717-720. [I21 D. B. Paul, “A speaker-stress resistant HMM isolated word recognizer,” in Proc. IEEE Int. Con$‘ Acoust., Speech, Signal Processing, Apr. 1987, pp. 713-716. 1131 R. P. Lippman, E. A. Martin, and D. B. Paul, “Multistyle training for robust isolated word recognition,” in Proc. IEEE Int. Con$ Acoust., Speech, Signal Processing, Apr. 1987, pp. 705-708. [I41 J. Makoul, “Linear prediction: A tutorial revkw,” Proc. IEEE, vol. 63, pp. 561-580, Apr. 1975. [15] R. M. Gray, “Toeplitz and circulant matrices: 11,” Tech. Rep. 6504-1, Inform. Syst. Lab., Stanford Univ., Apr. 1977. ‘

Synthesis of Spectral Densities Using Finite Automata

The interest for the synthesis of rational spectral densities by means of finite sequential machines (SM’s) is twofold: first, this synthesis may be applied to design encoders often used for spectral shaping purposes in digital data transmission or recording; second, it may be used for the numerical simulation of processes with assigned spectral density, as an alternative method to the classical approach of passing a white noise signal through a suitable digital filter. Let us recall that any rational spectral density may be represented in the form [l]

where * denotes conjugate and the zeros Ck and poles ek satisfy the conditions: I Fk 1 5 1, I ek 1 < 1. Moreover, if the process is real, to every zero and pole there corresponds its complex conjugate. The synthesis through a linear filter requires a factorization R(z) = o2h(z-’)h*(z*); then h ( z - ’ ) provides the transfer function of the filter and a’ gives the variance of the filter input. The general problem of synthesizing an assigned rational spectral density by means of an SM fed by a suitable input formed by independent and identically distributed (i.i.d.) symbols was dealt with by Mullis and Roberts [2], who proved that such a synthesis is possible for any rational spectral density. Unfortunately, the proof is not constructive and does not convey any suggestion about the choice of the SM and of the input probabilities. A complete solution to the problem has been given by Mullis and Steiglitz [3] for a particular case, namely, for those spectral densities which qn obtained by summing up elementary spectral densities haviqg only a single pole inside the unit disk. In their solution, each elefnentary spectral density requires a separate SM and the machine inputs must be independent. In the following we consider the same class of spectral densities as in [3], but use a different approach. Our solution, in which a fundamental role is played by circulant matrices as well, leads to a simpler structure involving a single sequential machine. 11. REVIEWOF THE SPECTRAL ANALYSIS OF THE OUTPUT O F AN

-+

s f + I = g(s,,

Carlo M. Monti, Gianfranco L. Pierobon, and Umberto Viaro

Abstract-A method for designing a finite automaton whose output exhibits a given rational power spectral density belonging to a particular class, is presented. The method exploits the properties of circulant matrices, which allow us to build a stochastic matrix with a set of assigned eigenvalues.

Manuscript received May 7, 1992; revised October 3, 1992. The associate editor coordinating the review of this correspondence and approving it for publication was Prof. Georgios B. Giannakis. The authors are with the Department of Electronics and Informatics, University of Padova, 35131 Padova, Italy. IEEE Log Number 9207551.

SM

As it is known, an SM (in particular, reference will be made to a Moore machine) may be specified [4]as a quintuple 3n = {a, a, S, g, h ) where 63 is the input set, Q. is the output set, S = {al, a,, . . ’- , a,} is the state set, g: S x 63 --t S is the state transition Q. is the output function; explicitly, function, and h: S

b,),

a,

=

h(sJ

(2)

where s,, b,, and a, denote the state, input, and output processes. If the input process b, is composed of i.i.d. symbols, the state process 9, is a homogeneous Markov chain [5] whose transition probability matrix ll can be determined from the probability mass function of the input and the state transition function of the automaton. By assuming that the Markov chain is ergodic, the state probability vectorp = [p(l),p ( 2 ) , . . . ,p ( l ) ] collecting the probabilitjesp(i) = Pr[s, = a,] is obtained from I

p =~II,

C p(i) = ,=I

1.

We recall that the transition probability matrix

1053-587X/93$03.00 0 1993 IEEE

(3)

II of an ergodic

IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL.

1994

Then, if W / & is the unitary matrix whose (k given by w(k)/&, we have

Markov chain is an irreducible stochastic matrix [6] which, by the Frobenius theorem [7], satisfies the following conditions: i) II has nonnegative entries; ii) the rows of II sum to 1; iii) II has no eigenvalue with magnitude greater than 1; iv) II has the simple eigenvalue X = 1; v) if II has more eigenvalues with magnitude equal to 1, say A,, A I , * . . , Ah - they are simple and are the hth roots of unity, i.e., A,,, = exp (i2?rm/h),m = 0, 1, . . , h - 1. The integer h gives the periodicity of the Markov chain, which is cyclic if h > 1 , and is acyclic if h = 1. If h > 1, then by a suitable joint permutation of rows and columns, II can be put into the following cyclic form where the zeros on the main diagonal are square matrices:

A

,,

.

. ..

. ..

0

0

0

IIhl

0

0

.

0

. ,h

O l

-

...

.

. . * XihZ,

0

where I , , Z2, . . , z h are identity matrices whose orders coincide with those of the zero matrices in homologous positions in (4). It has been shown [8] that the spectral density of the output process a , decomposes into two parts: a continuous component R,(z), and a discrete periodic component Rd(fl. The continuous part of the spectral density can be written in the form: - II)-’A A*(z-’Z - IIT)-’@d

R,(z)

(6)

i) A is the vectorA [h(u,),h(u2),* . . , h(u,)lTspecifying the output function, ii) Q = D - IITDII, in which D is the diagonal matrix D 4 diag { ~ ( l )~ ,( 2 1 ., * , PW}.

A*J:P~PJ,A

111. RELEVANTPROPERTIES

=

IPJA l2

OF CIRCULANT

(7)

MATRICES

Denoting by cJ,j = 0, 1, . . . , n - 1, the elements of the first row of a circulant matrix C of order n [3], it is easy to verify that the orthogonal vectors

,,,(k) = [ I ,

er2*k/n,

...

ef2r(n-l)k/n T

1

(8)

with k = 0, 1, . . . , n - 1, are right eigenvectors of C corresponding to the eigenvalues . .I J

=

0, 1,

. . . ,n -

1.

(9)

0, 1,

. . ,n

- 1.

(11)

Equation (6) shows that the poles of R,(z) only depend on II, whereas its zeros essentially depend on A. This suggests splitting the synthesis procedure into two steps: first, a matrix II with the desired eigenvalues (and possibly other eigenvalues) is looked for, and then the output vector A is determined so as to ensure the cancellation of the additional eigenvalues introduced in the first step and, at the same time, the realization of the desired zeros. Conceming the first step, as the poles of R,(z) correspond to the eigenvalues of II and their inverses, the problem is that of determining a stochastic matrix II whose eigenvalues include the poles ek of magnitude less than 1 of the desired spectral density. Since the poles p k are either real or in complex conjugate pairs, they may be arranged so as to form a Hermitian sequence. However, the Fourier transform of this sequence may exhibit negative values; in this case, one can resort to Corollary 1 and take the smallest value of L ensuring that the set { 1, &, . * , of;- ,} corresponds to a stochastic matrix C. To recover the desired eigenvalues, the following considerations are useful.

-

II-

el, = j = o c . e i 2 r k j / n , k

=

IV. SYNTHESIS PROCEDURE

The discrete component Rd(f) exhibits spectral lines, located at the frequencies fm k m / h , m = 0, 1, * . . , h - 1, and having amplitudes =

(10)

,

where z = * denotes conjugate transpose, Z is the identity matrix of order I, and

R:”

C = WAW-’

Furthermore, since the sequence {c, }, j = 0, 1, * * * , n - 1, is real, the eigenvalues e,, form a symmetric Hermitian sequence. It is also easy to see that, if the inverse of C exists, it is itself a circulant matrix. As we shall see shortly, our aim is to realize a stochastic circulant matrix which exhibits preassigned eigenvalues. Therefore, it is necessary that (1 1) supplies nonnegative real entries c,. Now, if the set of these eigenvalues includes, in addition to 1, pairs of complex conjugate eigenvalues (with magnitude less than 1) and if they are ordered so as to form a Hermitian sequence, then (1 1) indeed provides real entries. To ensure that these are nonnegative too, further conditions must be satisfied. A sufficient condition is provided by the following theorem whose proof [9] is omitted for brevity. Theorem 1: Given a set of n eigenvalues including unity and forming a Hermitian sequence, there exists a circulant stochastic matrix which exhibits these eigenvalues if the sum of their mag0 nitudes, excluding the eigenvalue 1, does not exceed 1. Of course, a preassigned set of eigenvalues may not satisfy the above conditions directly. However, as we shall see in the next section, by means of a suitable “contraction” procedure, we shall be able to obtain from the original set a new set meeting the conditions under consideration, and then to build a stochastic matrix with the required eigenvalues (perhaps along with some unwanted ones which will be disposed of later on). A direct consequence of Theorem 1 is the following. Corollary I : Given a set of n eigenvalues { 1 , e , , . . , e, - } forming a Hermitian sequence, there exists an integer L such that the set (1, e:, . . * , e ; - , } may be realized by a circulant sto0 chastic matrix.

(4)

* *

C ek e - J 2 * k J / n , j

k=O

IIh-1

For later use, it is convenient to introduce, form = 0, 1, - 1, the matrices

W-ICW,

n-I

nc =

0

... ...

=

+ 1)th column is

where A = diag {eo,e l , * . . , e n - ,}. From (9), it follows that the first row of C is the discrete Fourier transform of the ordered sequence (eo,e l , * . . , p, - I ) of its eigenvalues, except for the multiplicative constant n, i.e.,

0

n = l ..

41, NO. 5 , MAY 1993

~

1995

IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 41. NO. 5, MAY 1993

Making use of the same decomposition for the remaining two factors, we get

Let us consider the cyclic matrix

"=L!

:. :.

0 0

0

c

0

0

;]

. * .

1

R,(z) = - [za

nL

T + z2a; +

(12)

where, due to the lack of space, only the generic kth term of the diagonal matrix has been shown. Now, denoting by b, = [bmo,bml, . . . , b,,,, , I T the discrete Fourier transform of a,, we have W-'a, = b,,, so that (16) may be rewritten as

. [Z-Ibl, whose eigenvalues are the eigenvalues of C repeated L times. It follows that the eigenvalues of ll are the Lth roots of them. In particular, these eigenvalues comprise the required set { 1, el, . . . ,

e,,-Il. The problem is now to find a vectord such that: i) the amplitudes of the spectral lines (7) are set to zero, ii) the additional poles are cancelled, and iii) the desired zeros are realized. By assuming that ll has the structure (12) and partitioning the vector p into L subvectors of dimension n , i.e., p = [ p Ip 2 . . . p L ] , (3) is satisfied by p I = p 2 = . . . = p L , p L = p L C , so that every component of p is equal to 1 / n L and we simply have D = (1 / n L ) I . Conceming i), from (7) we obtain pJ,J = 0, m = 0 , 1, * * * , L - 1, and after partitioning the vector A like vector p , i.e., A = [a:.,' . . . we have

which implies that the sum of the elements of every subvector a, must be zero. We thus obtain L conditions on the elements of the output vector A. Conceming ii) and iii), we now derive the structure of the continuous part of the spectral density corresponding to the considered transition probability matrix (12). Using the above partition, after cumbersome algebraic manipulations, from (6) we get

nL

+ z52,*]

... ...

where the number of blocks in a row or column is L (the period of the Markov chain) and the matrices Io are identity matrices of order n. The Lth power of ll is

1 R,(z) = - [zaT

* *

+ z2a: + . . + $a,*] *

. [Io - t.cT]-'[Io - C T C ] [Io - z-Lc1-1 . [ z - l a , + 2-2u2 + . . .

+ z-La,].

(15)

The three factors appearing in the middle of (15) are circulant matrices, so that they can be diagonalized by means of the unitary matrix W/&. With reference to the first of them, we obtain

+ ~ - ~ b 2+r

'

..

+ z-LbLk]

(17)

where the term fork = 0 is missing since eo = 1. The residues at the poles to be cancelled are zero only if [bTk zbtk . . . + z L - ' b & ] is a factor of [ l - ( z ~ t )and ~ ] the same happens for the terms in z - ' ; therefore, the unknowns bmksatisfy the relations

+

+

b,,h = b , , e ; " - l ,

m

=

2 , 3,

...

, L, k = 1, 2,

...

,n - 1 (184

and (17) reduces to

By equating the residues at the poles values r,, we find

ek

of (19) to the required

so that

The unknowns b,nO,m = 1, 2, . . . , L, which do not appear in (18), are all equal to 0 due to the condition on the magnitude of the spectral lines previously considered. Observe that, in order for (18b) to be satisfied, the ratio r L / e k must be positive; in other words, the given spectral density can be realized by the suggested procedure only if the arguments of the residues match those of the corresponding poles. This condition implies that R , ( z ) is the sum of n terms, each of which is itself a spectral density (cf. (20)). The same synthesis limitation is inherent in the method presented by Mullis and Steiglitz in [ 3 ] where, however, the realization of each pair of reciprocal poles requires a separate machine and the machine inputs must be independent. To summarize, the proposed synthesis procedure for the spectral densities that may be expressed in the form (20) entails the following steps: i) order poles ek, 1 1 < 1, so as to form a Hermitian sequence and find the minimum integer L such that the eigenvalues (1, e:, . * * , ef;- I) satisfy the condition of Theorem 1;

IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 41,

1996

ii) compute the Fourier transform of the sequence (1, e!, * , ef;- ,), thus determining matrices Cand ll (cf. (1 1) and (12)); iii) compute the residues at the above poles and find the values Of bmk Using (1 8); iv) determine the components of the output vector A as partitioned in (14) by inverse Fourier transformation:

NO.

5 , MAY 1993

for determination of the angles and distances of the radially distributed contents of the autoterms in the ambiguity domain. The proposed kernel effectively reduces the cross-terms and noise for linear FM signals. The result is a tool for high-resolution time-frequency representation of nonstationary, primarily linear FM signals.

I. INTRODUCTION

.“-I

j = 0, 1,

. . . ,n

- 1.

(21)

V. EXAMPLE Let us consider the synthesis of a spectral density (20) with n = 4 , e l = (1 i ) / 3 , ez = 1/2, e3 = (1 - i ) / 3 and positive values of the ratios rk/& As it can easily be verified, the discrete Fourier transform (11) of [ l , e,, e*, e3] leads to a circulant matrix with some negative entries. However, according to Corollary 1 with L = 2, there exists a circulant stochastic matrix C with the eigenvalues 1, e : = 2i/9, e: = 1/4, e: = -2i/9. From ( l l ) , the entries of the first row of this circulant matrix turn out to be

+

1 -[45 144

11 45 431.

This enables us to find the transition probability matrix ll as in (12) with L = 2. In this way, the resulting SM has nL = 8 states. The output vector A can finally be obtained through (18) and (21). REFERENCES

J. L. Doob, Stochastic Processes. New York: Wiley, 1953. C. T. Mullis and R. A. Roberts, “On weak equivalence of linear systems and finite state systems,” SIAM J . Math. Anal., vol. 10, no. 3, pp. 498-511, 1979. C. T. Mullis and K. Steiglitz, “Circulant Markov chains as digital signal sources,” IEEE Trans. Audio Electroacoust., vol. AU-20, no. 4, pp. 246-248, 1972. Z. Kohavi, Switching and Finite Automara Theory. New York: McGraw-Hill, 1970. A. Paz, Introduction to Probabilistic Automata. New York: Academic, 1977. I. G. Kemeny and J. L. Snell, Finire Markov Chains. New York: Van Nostrand, 1960. P. Lancaster, Theory ofMatrices. New York: Academic, 1969. G. Bilardi, R. Padovani, and G. L. Pierobon, “Spectral analysis of functions of Markov chains with applications,” IEEE Trans. Commun., vol. COM-31, no. 7, pp. 853-861, 1983. C. M. Monti and G. L. Pierobon, “Circulant stochastic matrices with assigned eigenvalues,” Dep. Electron. Informatics, Int. Rep., Univ. Padova. 199 I .

Any bilinear time-frequency distribution (TFD) can be expressed as the two-dimensional Fourier transform of the product of the symmetrical ambiguity function of the signal (in text referred to as the ambiguity function, AF) and the kernel ~ ( 0 7) , [7]: P,(t,f) =

jymlmA , @ , -m

7)cp(0,

where the ambiguity function A, (0, A,@,

7)

1-

=

-m

z(t

+ 7/2)

7)

*

7)e12*ere-1zrfi de d7

(1)

is defined as [7]

z * ( t - 7/2)eJZdr dt.

(2)

A number of kernels were proposed [2], [7] in order to suppress the influence of the artifacts and noise. Some of them are fixed (signal independent) [5], [ 131 the others are adaptively designed for the signal under analysis (signal dependent) [I], [LO]. Only the signal independent kernel based TFD’s can preserve the desirable properties of the WVD [6]. However, since the kemel function is actually the transfer function of a two-dimensional filter (dimensions are 0 and 7 ) , for suppression of artifacts and noise, the kernel design (as the filter design) has to be signal dependent in order to obtain good performance for a broader class of signals. The problem addressed in this correspondence is a signal dependent kemel design for time-frequency signal analysis. In a variety of publications on the subject it is important to mention the work by Baraniuk and Jones [l] who treated the kernel design as an optimization problem. In this correspondence we propose a simpler approach based on the Radon transform of the modulus of the AF .

11. THEORETICAL CONSIDERATION FOR LINEARFM SIGNAL COMPONENTS

A . Background We model each signal segment as K

Z(t) =

C

k= 1

&(f)

+ n(f),

(0 5

f 5

T)

where K is the number of signal components of ~ ( t )n, ( t ) is the noise, and Tis the segment duration. We further model each component as a finite duration unit amplitude quadratic phase signal:

Design for Time-Frequency Using the Radon Transform Branko Ristic and Boualem Boashash Abstract-This correspondence presents a new kernel design technique for time-frequency signal analysis. The technique is based on the Radon transform of the modulus of the ambiguity function of the signal

Manuscript received December 5, 1991; revised June 29, 1992. The associate editor coordinating the review of this correspondence and approving it for publication was Prof. Miguel A. Lagunas. This work was supported by the Defence Science and Technology Organisation, Australia and the Australian Research Council. The authors are with the Signal Processing Research Centre, Queensland University of Technology, Brisbane, 4.4001 Australia. IEEE Log Number 9207544.

0 <

fok 5 f 5 tok

+ Tk 5

T.

(3)

The signal component k starts at tok and exists during the period Tk. For simplification, unit amplitudes are assumed. In case the signal components have amplitudes, only step 2 of the design procedure change (Section ”’). From the definition of the instantaneous frequency of a signal component, fk (t) = (1 /2?r) ( d + k ( t ) / d t )we get

f k ( f ) = alk

1053-587X/93$03.00 0 1993 IEEE

+ a2kf

0

5 tok 5 f 5 tok

+ Tk 5

T.

(4)