Illustration of the effects of sampling on higher-order spectra

SIGNAL PROCESSING ELSEVIER

Signal Processing 36 (1994) 375-390

Illustration of the effects of sampling on higher-order spectra Joel

Le Roux*, Christophe Coroyer, Delphine Rossille

13S, University of Nice, CURS . bat . 4, 250 rue A . Einstein, Sophia Antipolis, 06560 Valbonne, France

Received 3 August 1992 ; revised 5 June 1993 and 10 November 1993

Abstract This paper illustrates some of the consequences of sampling on higher-order spectrum (HOS) analysis, and especially the important difference between deterministic and random stationary signals . In the case of deterministic signals, it provides the conditions enabling the construction of the sampled signal HOS by periodization of the analog signal HOS, in the same way as for random stationary signals . In this last case, it insists on the precautions to be taken in applying digital system identification techniques to actual analog signals . It also proposes a method for performing the simulation of stationary analog random signal sampling as a mean of checking the validity of identification algorithms .

Zusammenfassung Dieser Artikel erlautert einige der Konsequenzen der Abtastung auf die Analyse von Spektren hoherer Ordnung and besonders den wichtigen Unterschied zwischen deterministischen Signalen and stationaren Zufallssignalen . Im Fall von deterministischen Signalen werden die Vorraussetzungen dafur geschaffen, um die Entwicklung des Spektrums hoherer Ordnung von abgetasteten Signalen durch die Periodisierung des Spektrums hoherer Ordnung von analogen Signalen, auf die gleiche Weise zu erlauben, wie fur stationare Zufallssignale . Im diesem letzten Fall wird auf die Einhaltung von VorsichtsmaBnahmen bei der Anwendung von Methoden zur Identifikation digitaler Signale auf in Wirklichkeit analoge Signale bestanden . Es wird auch ein Verfahren zur Simulation der Abtastung von stationaren analogen Zufallssignalen als eine Methode vorgeschlagen, um die Gultigkeit von Identifikationsalgorithmen zu uberprufen .

Resume Cet article illustre quelques-unes des consequences de l'echantillonnage sur l'analyse spectrale d'ordre superieur (HOS), particulierement la difference importante entre signaux deterministes et signaux aleatoires stationnaires . Dans le cas deterministe, nous donnons les conditions permettant la construction du HOS du signal echantillonne par periodization du HOS du signal analogique, de la meme maniere que dans le cas aleatoire stationnaire . Dans ce dernier cas, nous insistons sur les precautions devant titre prises lots de l'application des techniques d'identification de systemes digitaux a des signaux analogiques reels. Nous proposons egalement one methode de simulation de 1'echantillonnage de signaux analogiques aleatoires en tant que moyen de verification de la validity des algorithmes d'identification .

"Corresponding author . 0165-1684/94/$7 .00 '1994 Elsevier Science B .V . All rights reserved SSD1 0165-1684(93)E0123-3

376

.1. Le Rmcr ei al . i Signal Processing 36 (1994) 375-390

Key words. Bispectrum; Deterministic signal; Random stationary signals; Effects of sampling ; Higher-order spectrum ;

Pseudo-random sequences ; Stationarity ; Third-order spectrum ; Third-order correlation ; Third-order pseudo-random sequences; Triple correlation

1 . Introduction

Most results concerning second-order spectrum analysis of sampled signals share the same expressions in the case of deterministic and random stationary signals . This is generally not valid for higher-order spectrum (HOS) analysis of sampled signals . For example, Brillinger and Rosenblatt [1] have shown that the third-order spectrum (TOS) of the sampled signal deduced from a band-limited analog random stationary signal vanishes in the triangles ABC and EFG of Fig . 1 ; it is obtained by periodization of the TOS of the analog signal . (by `analog' signals, we mean functions defined for a continuous variable .) This is generally not true for deterministic or random non-stationary signals ; the TOS of these signals may be non-zero in these triangles . Hinich uses this property in order to detect `transient' signals [5] . It appears that higher-order spectra of sampled deterministic signals often hold more information than those computed after sampling random stationary analog signals . A specific reconstruction algorithm may be deduced from the properties of these deterministic signals HOS . For example, in the case of analog deterministic or random stationary signals, Matsuoka and Ulrych obtain the best least-squares estimate of a signal Fourier transform phase from the TOS phase through the solution of a set of linear equations [8, 12] . In the case of deterministic sampled signals, this solution is expressed simply as an average [6] . This last type of algorithm could not be developed when using only third-order spectral data measured on the analog signal, that is to say data in ACDEGH . The second section illustrates this difference between deterministic and random sampled signals in one important case, when higher-order spectra are factorizable (see [11, 13]). The main result on the effects of sampling is that the TOS of a sampled deterministic signal can be deduced by periodization of the TOS of the initial analog signal only if

the range of this signal is less than a third of the sampling frequency . This is detailed in Section 2 .1 . Some examples are given in Section 2.2 . The result is extended to higher-order spectra in Section 2 .3 . If we consider random signals, there is an important difference between two classes of `stationary' sampled signals : those that are obtained by sampling a band-limited random stationary signal and those that are not, although the sequence of samples, given for discrete values of time, is stationary . We think that it is important to insist on this difference for the following reason : many recent publications on HOS system modelling such as those mentioned in tutorial papers [3, 4, 9] generally consider the analysis of sampled signals and not the problem of sampling itself. The validation of the algorithms is obtained through digital simulations and it seems that they are very seldom applied (after sampling), to actual analog signals . It appears that several algorithms proposed in these papers, especially those on MA modelling, are developed for the second class of signals and cannot be applied (after sampling) to band-limited random stationary signals . These algorithms often suppose the representation of the TOS of the sampled signal in a factorized form . For example, the validity of a finite-order MA model supposes that the TOS is non-zero almost everywhere . This is generally not compatible with the periodization of the TOS of band limited stationary analog signals as in [1] . The third section illustrates this fundamental difference between the two classes of random signals in the context of system identification . In Section 3 .1, we recall this problem of system identification based on HOS factorization . We also give the conditions allowing the simultaneous factorizability of the TOS of analog and sampled signals, for which there is no difference between deterministic and random stationary signal TOS . In Section 3 .2, we give an example showing that a signal obtained by (digital) filtering a sequence of random independent samples cannot be the sampling of a third-order



.1. Le Roux of al . / Signal Processing 36 (1994) 375 390

stationary band-limited analog signal . However, in order to check a numerical algorithm before applying it to actual signals, one needs to generate a sampled pseudo-random sequence that is equivalent to the sampling of a third-order stationary analog signal . In Section 3 .3, we propose a method for constructing such a signal . Finally, Section 4 discusses several conditions allowing the correct application of identification algorithms, In the following we recall some basic definitions and give the notations for deterministic and random signal TOS .

V

D

C

B

E

A

F

H

G

Fig . 1 . Third-order spectrum support for deterministic bandlimited and sampled signals . Dotted lines indicate the symmetry axes : the non-redundant region is OAK in the case of analog band-limited signals and DAL in the case of sampled signals .

1 .1 . Third-order correlation and specs m o deterministic signals The third-order correlation (TOC) of an analog deterministic real signal f,(t) is r,(x,y)=J f(t)f(t+x).fo(t+y)dt,

37 7

is the Fourier transform of the (N + I )th-order correlation

(1)

T

N

where T is the support of f;(t) . Its Fourier transform is the TOS that may be computed from the Fourier transform F J (u) of f (t) : Ba(u, v) = Fa(u)Fa(v)Fa( - u - v) .

(2)

We assume that Fa (u) is zero outside ] - rc, tc[ . Then the support of B,(u, v) is the hexagon ACDEGH of Fig. 1, (3)

Jul, Ivl, In + vl < R .

For integer values oft, f (t) is sampled in f(k) with (periodic) Fourier transform F(u) and TOS B(u,v)=F(u)F(v)F(-u-vmod2it) .

(4)

B(u, v) is the Fourier transform of the (discrete)

TOC of f(k), rim, n)=Ef(k)f(k+m)f(k+n) .

(5)

k

r .v+I(mt, . . .,MN)=Y_ k

1 .2 . TOC and spectrum signals

f(k) 6fik+m^l .

(7)

n-1

of stationary random

In the case of stationary random analog signals, the TOC is ra(x,y)=E{f(t).fa(t+ .x)f(t+Y)} .

(8)

When we consider random signals, the averages E{ } are computed on independent realizations of the random process. The important problem of higher-order ergodicity is not considered in the present communication . The TOS Ba (u, v) is the Fourier transform of ra (x, y). When f(t) is sampled in f(k), Brillinger and Rosenblatt [1] have shown that B(u,v), TOS of f(k), Fourier transform of the TOC

Similarly, the (N + 1)th-order spectrum of f(k), rim, n) = E { f(k)f(k + m) f(k + it)}, BNII(nt . . . ,u,V) N

= F

- ~, u ^ mod2it

M1'

IJ F(u ^ ),

(6)

(9)

vanishes in ABC and EFG of Fig . I and is deduced from B 5 (u, v) by a simple square periodization of BDFH . (This result is generally not true in the case



378

J. Le Roux el al.

of deterministic signals .) When able in

B a (u,

Signal Processing 36 (1994) 3 75 -390

v) is factoriz-

Ba (u, v) = G5(u)G0(v)G,( - u - v),

We intend to give the conditions for which B(u, v) is equal to the periodization (following a square geometry) of Ba (u, v), that is,

(10) Bt(u, v)

fa (t)

may be considered as the output of a linear filter with frequency response Ga (u) driven by a non-Gaussian zero-mean third-order stationary white noise . G a (u) may be non-minimum phase and an important application of higher-order spectral analysis is the reconstruction of G i (u) from measures of Ba (u, v). However, the factorizability of B a (u, v) is generally not equivalent to the factorizability of B(u, v) . Similar definitions and properties are directly extended to the case of higher-order spectra and correlations .

2 . Effects of sampling on the TOS of deterministic signals

_

Y,

Yr

Fa(u - 2Pn)Fa(v - 2qn)

p - -"^ q-

x F a (- u -

v + 2pit + 2qn) .

B'(u, v) is the Fourier transform of the signal obtained by applying a square geometry sampling to ra (x, y). (A detailed presentation of 2-D sampling is given in [2] .) B(u, v) and B'(u, v) being periodic of period 2n in u and v, it is sufficient to analyze their restriction to the square support [ - n, n[ x [ - it, it[. The equality holds if no couple (u, v) belongs simultaneously to this square and to an oblique strip where F a ( - u - c - 2g 3 n) is nonzero for a non-zero value of q 3 . Suppose that F 4 (u)is non-zero only for Jul
2 .1 . Relation between the third-order spectra of analog and sampled deterministic signals

(11)

B(u, v) is the TOS deduced from F(u) . Replacing F(u) by the periodization of F 0 (u) in (4), we get B(u,v)=

Fa(u L 4L =-

-

2grx) J

r

® Fa(v - 2g2'n)

Y

Support of de Fouler transform of Ne sampled triple rorrelstimt ofthocnnunuauslignal

J 11111

x

Fa(-u-v L qJ --

U

-

2gsr<)J .

(12)

(14)

Then the support of B(u, v) is the repeated hatched pattern in Fig . 2 (the hexagon and the two

The basic theorem on signal sampling states that the reconstruction of a(t) from f(k) is possible if the support of Fa (u) is in a range less than half the sampling frequency (Jul < n) . In this case. there is no spectral aliasing and F(u) is given by the periodization of Fa (u) . Then, the support of B a (u, v) is the hexagon (3), ACDEGH in Fig. 1, but the periodization of Ba (u, v) is not necessarily equal to B(u, v) . This is valid only if the support of F a (u) is Jul < 2t/3 .

(13)

Support of the brspccwm of a sanded signal (Pounce transform of the triple correlation of the sampled signal

Fig . 2 . Overlapping effect due to sampling .



J. Le Roux et al . / Signal Processing 36 (1994) 3 75 390

triangles) where the three factors of (12) are nonzero simultaneously while the support of B'(u, v) is given by the repetition of the hexagon only . The couples (u. v) belonging to [ - n, n [ x [ - it, n [ may strip also belong to an oblique lu + v - 2g 3 rz1 < n only if q 3 = + 1 . If there is no couple (u . v) of [ - n, 7c[ x [ - n, n [ so that 2n-h < lu+vl,

(15)

there is no overlapping of the three (horizontal, vertical and oblique) strips and at least one of the three factors will vanish in ABC and EFG . If (15) does not hold, there may exist an interaction between a frequency u in (0, h), a frequency v in (0, h) and a frequency ( - u - v mod 2n) in (x, 2n - h) yielding a non-zero signal in the two triangular supports, ABC and EFG . According to (14), the largest value of In + vl is 2h and B'(u, v) = B(u, v)

379

ation, the Fourier transform of f(k) is non-zero in (h, 2n - h) ; consequently, its TOS is non-zero in ABC .

2.2 .2 . Third-order spectra of band-limited signals The TOS of the sampled signal 2 sin(nk/2)/nk is

the periodization of the TOS of the analog signal 2 sin(nt/2)/nt following a square pattern; the TOS of 4 sin(3xk/4)/ink is not obtained by periodization of the TOS of 4 sin(3nt/4)/ant . In the first case, (11) holds, while it does not in the second . We have computed the TOC of the two sampled signals 2sin(xk/2)/xk and 4sin(3nk/4)/3xk; the TOS, Fourier transforms of these TOC, show this effect of sampling (Fig . 4).

2 .3 . Extension to higher-order spectra

(16)

only if the bandwith of the signal is one-third of the sampling frequency, h < 21r/3 .

(17)

Then, r(m, n) is obtained by sampling ra (x, y) . If (17) does not hold, B(u, v) and B'(u, v) coincide inside the hexagon (3) but not in the two triangles - n< I u l, I v I< n, I u+ v(> n and r(m, n) is not obtained by sampling r a (x, y).

The commutativity of the operations of sampling and computing the TOC is possible only if the square of width 2h and centered at the origin Jul, Ivl < h < n does not intersect the lines In + vl = 2n - h. Similarly, in the case of the Nthorder spectrum, if Fa (u) has a finite support Jul < h < n, the commutativity holds only if the hypercube lull, . . . , luxl

< h

(18)

does not intersect the hyperplanes Ju l

2 .2 . Examples

+ . . . + u,l = 2n - h,

(19)

that is, only if 2 .2 .1 . Third-order spectra

of band

pass signals

There exist analog signals (for example cos 2nt/ 3 + 0) with null TOC while the TOC of the sampled signal is non-zero (here cos 2n(m + n)/ 3 + 0) . The TOS of this last function, computed after application of a weighting window, is shown in Fig . 3 . This type of phenomenon may appear in particular when the spectral support of the signal is between a fourth and a half of the sampling frequency : Consider a band-limited real signal f (t) with spectral support (h, n) with 0 < h < x . Its TOS is zero if h > x/2 . When it is sampled at rate 1, there is an overlapping effect in ABC of Fig. 1 . By periodiz-

h < 2x/(N + 1).

(20)

3 . Third-order spectra of random sampled signals In this section, we intend to illustrate the important difference induced by sampling between the HOS of random stationary signals corresponding to the hypotheses of Brillinger and Rosenblatt [I] and the HOS of non-stationary analog signals which appear `stationary' after sampling . A signal obtained by digital filtering a sequence of independent and identically distributed (i .i .d) samples



380

.1. Le Roux et al. / Signal Processing 36 (1994) 375-390

COs

I

I

T7'

COS

1

T-

t

0 .8 0 .6 0 .4 0 .2

-10

0

10

20

30

3 .a. Analyzed signal (after application of a weighting window)

3 .b . Bispectrum modulus Fig. 3 .

Third-order spectrum

of

a sine wave whose frequency

belongs to this last class and cannot be obtained in sampling a band-limited stationary random signal . In fact many studies on HOS, especially in identification, are confirmed through simulations on this kind of signals . This observation yields two consequences : - the conditions for applying identification algorithms to actual signals must not be disregarded. This point is illustrated in Sections 3 .1 and 3.2 . - A method for generating a signal corresponding to the stationarity hypotheses of Brillinger and

is

one-third

of

the sampling frequency .

Rosenblatt is needed. Such a method is proposed in Section 3 .3 .

3 .1 . TOS foetorizability after sampling An important class of potential applications of higher-order spectrum analysis is the identification of non-minimum-phase, linear, time-invariant systems with transfer function G a (u) from the measurement of higher-order statistics of its output f,(t) . The underlying model in most HOS identification



J . Le Roux el al . I/ Signal Processing 36 (1994) 375-390

38 1

4 .a.1 . Analyzed signal (after application of a weighting window)

4 .a .2 . Triple correlation

4 .a.3 . Bispectrum modulus

Fig. 4 . Third-order spectrum of sine functions . (a) Third-order spectrum of a sampled sine function whose largest frequency is the fourth of the sampling frequency : there is no overlapping effect . (a .l) Analyzed signal (after application of a weighting window). (a .2) Third-order correlation. (a.3) Third-order spectrum modulus .



3 82

J. Le Roux et al. i Signal Processing 36 (1994) 375-390

4 .b .l . Analyzed signal (after application of a weighting window)

4 .b .2 . Triple correlation

4 .b .3 . Bispectram modulus Fig. 4 . (b) Third . order spectrum of a sampled sine function whose largest frequency is 3118 of the sampling frequency : the overlapping effect appears . (b .l) Analyzed signal (after application of a weighting window). (b .2) Third-order correlation . (b .3) Third-order spectrum modulus .



J . Le Roux et al .( Signal Processing 36 (1994) 375-390

techniques consists in stationary non-Gaussian white noise w(t) feeding the analyzed system and producing fa(t) . Then the TOS of a ll) is

f

Ba (u,v)=G a (u)G a (v)G a (-u-v)E{w 3 } .

noise (independent samples) w(k) with constant spectral density in ( - n, n) . The digital output signal f(k) is

(21) f(k)

E{ } are averages on several independent realizations and not an average in time on one experiment. In order to perform the identification, one may use digital techniques, and needs to compute the third-order statistics of a stationary sampled signal . If f (t) is stationary and band-limited in ( - h, h) with h < n, and sampled in f(k), then B(u, v), TOS of f(k), vanishes in ABC and EFG and is the periodization of B a (u, v) [1] . A development similar to that of Section 2 shows that, when Ba (u, v) is factorizable as in (21), B(u, v) will be factorizable in B(u, v) = G(u)G(v)G( - u - v).

3 83

(22)

where

_ Y

g(a)w(k - a),

(24)

where g(a) are the samples of the inverse Fourier transform of G(u) . In this case the sampling effect appears just as in the case of deterministic signals ; the TOC of f(k) is r(m,n)=E{ Y la=-w

g(a)w(k-a)

xw(k+m-/f)

Y

g(fl)

q(y)w(k + n - 7) (25)

or r(m,n)=

L g(a) L- „ g(f) Y g(!) B=

G(u)

_ Y

Ga(u - 2pn),

(23)

xE[w(k-a)w(k+m-/i)w(k+n- )} .

P-

only if the sampling frequency is larger than 3h . Many digital reconstruction algorithms are based on the factorizability (22) and special care should be taken when applying to f(k), as the periodization of B a (u, v) is generally not factorizable . However, the sampled signal may be `stationary' (considering only the sampling instants) while the band-limited signal to be sampled is not stationary . In those cases, the TOS is not necessarily zero in the triangle ABC of Fig . 1 and B(u, v) may be factorizable. As we show in the following example, this is what occurs in simulations where .f(k) is generated by feeding a digital filter with a sequence of samples .

3 .2 . Filtering third-order stationary digital i .i.d. sequences

We show an example where the TOS is factorizable as in (22) (see [13]), i .e . when the analyzed signal f(k) may be considered as the output of a linear system with impulse response g(k) driven by a non-Gaussian stationary zero-mean digital

(26) As w(k) is stationary, white, zero-mean and nonGaussian (we suppose that E{w 3 } # 0), the math-

ematical expectation is non-zero only when

- a=m - lj=n - y

(27)

So, r(m, n) = E{w 3 }

Y,

y(a)g(m + a)g(n + a) . (28)

r(m, n) is proportional to the TOC of g(a) . Hence, if

the overlapping effect due to sampling and described in the case of deterministic signals appears for g(ot), it will also appear in the TOC and TOS of the sampled random signal f(k) . This shows that a pseudo-random sequence cannot be interpreted as the result of sampling a third-order stationary band-limited analog signal. However, (28) shows that in this case, the application of the algorithms developed for factorizable deterministic signals is licit . EXAMPLE. Consider a non-Gaussian zero-mean noise w(t) made of independent samples taking the value 2 with probability 1/3 and - 1 with



3 84

I Le Roux et al. / Signal Processing 36 (1994) 375-390 2 15 1 0 .5 0 -0 .5

0

.30

.20

-10

0

10

20

30

40

-10

0

10

20

30

40

a. Zero mean white noise 3 2

0

-3 -40

-30

-20

b .Filtered noise Fig. 5. Generation and filtering of a non-Gaussian digital `white' noise. (a) Zero mean white noise (input). (b) Noise filtered by I - z

probability 2/3 for integer values oft and 0 for noninteger values of t (Fig . 5) . w(t) is filtered by an analog filter with frequency response Ga(u)= 1 -exp(-ju) for Jul S it, for lul > n, G,(u) = 0

(29)

yielding a band-limited random signal J;(t). Then the TOS of f,(t) is non-zero only in the hexagon ACDEGH : U v u+ll B,(u, v) = 8j sin 2 sin sin 2 2 for Jul, lvl and lu + vl < a, (30) B,(u, v) = 0 for Jul, Ivl or lu+vl > t.

When f,(t) is sampled for integer values of t, it is equivalent to a signal f(k), result of filtering w(k) by the filter 1 - z - ' (Fig. 5) . The TOC r(m, n) of f(k) is computed (after sampling) ; its samples are given in Table L The Fourier transform of r(rn, n) gives the TOS of f(k) :

B(u, v) = 8j sin u sin

V

u sin

2 for all u and v . (31)

It cannot be obtained by periodization of B,(u, v) . This has been confirmed by simulation . The TOS modulus of the impulse response of I - z - ' is given in Fig . 6; the values of B(u, v) in the two triangles



J . Le Roux et al. I Signal Processing 36 (1994) 375-390

385

Table I Third-order correlation of a digital white noise filtered by 1 - z ' (all other values are equal to 0) 0

1

Fig . 8 . Third-order correlation of a computer generated pseudo-random sequence as described in Section 3 .2 . It presents significant peaks in (-31, -3) and in the five positions deduced from it by symmetry.

REMARK. The poor quality of the estimated TOS

Fig. 6 . Third-order spectrum of the impulse response of the filter I - - - ' (modulus) .

Fig. 7 . Estimated third-order spectrum of J(m), non-Gaussian digital white-noise filtered by I - z - ' (modulus) .

ABC and EFG cannot be obtained by periodization from the hexagon ACDEGH . This effect remains when one considers the TOC and TOS of the filtered random signal as described above (Fig . 7) . This TOS is the Fourier transform of the average of 100000 TOC; note that the estimation is still very noisy .

is mainly due to the imperfections of the pseudorandom noise generator . These imperfections appear clearly in the TOC of w(k) computed on an IBM RS 6000 . We generate a pseudo-random number w ra „ d (k) between 0 and 1 ; if w ra „d (k) 2/3 then w(k) = 2/3, else w(k) = - 1 /3 . The corresponding TOC is shown in Fig . 8 . It presents significant peaks in ( - 31, - 3) and in the five positions deduced from it by symmetry . The characteristics of pseudo-random generators may probably be found through the analysis of higher-order statistics . These generators are to be used with caution in HOS simulations .

3 .3 . Simulation and sampling of third-order stationary pseudo-random sequences

In order to check the validity of identification algorithms before applying them to actual stationary analog signals, it is very useful to build a pseudo-random sequence that can be considered as the sampling of a third-order stationary random analog and band-limited signal . In this section, we propose the construction of such a signal . We generate a random non-Gaussian zero-mean sequence of independent samples s(pAt) taking for example the values 2 with probability 1/3 and





386

.L Le Roux et al . / Signal Procevsing 36 (1994) 37 .5-390

- I with probability 2/3 for integer values of p and 0 elsewhere as in Section 3 .2 (Fig . 5) and apply on it a low pass filter with impulse response sin(nt)/nt in order to allow its sampling at period 1 . In the next subsection, we show that this analog signal is stationary. In Section 3 .3 .2, we use this signal to show an illustration of the results of Brillinger and Rosenblatt on the effects of sampling .

A shift of the variables t, t' and t" yields R a (u, v, w) = E{s 3 } H,(u)H,(v)H,(w)

exp[j(u + v + w)p] At, p

(36)

=- A

where H,(u)=I

3 .3 .1 . Stationarity

if ui cm

and

Ha (u)=0

In order to check the stationarity of the bandlimited random analog signal z(t)= v=

- z'

stn n(t - p At) s(PAt) n(t-p At)

Hence, (32)

we follow the approaches of Brillinger and Rosenblatt [1] and Hinich [5] and compute its 'thirdorder covariance' : E {z(t)z(t')z(t") }

v=-x

R a (u,

Mt ,

~'

p'At)

s(p"At) d

sin Tilt" - p" At) X (t „ n

(33)

As the samples s(p At) are independent, (33) reduces to 11

E{z(t)z(t')z(t")}=E{s3}

Y

sin tilt - p At) n(t-pAt)

l H.

(38)

wl < n

and

(u + v + w)At = 2qn,

(39)

where q is integer (see Fig . 9) . If At < 2/3, (38), (39) holds only for q = 0 . Then, R a (u, v, w) is non-zero only for u + v + w = 0 and n
x

v, w) = 0,

unless Jul,

a(pA,)sin 7c (t -pAt) n(t - p A t )

sin n (t' - p' At) s(p'At) XP'

if lul>n . (37)

v) = R,(u, v, - u - v)

sin n(t' - p At) sin n(t" - p At) n(t'-pAt) n(t"-pAt)

= Ha (u)H a (v)Ha ( - u - v)E {s 3 },

(40)

(34) non-zero in ACDEGH only . If At > 2/3, (39) has solutions for q # 0, Ra (u, v, iv) has non-zero values outside u + v + w = 0 and z(t) is not stationary .

The Fourier transform of E{z(t)z(t')z(t")} is R,(u, v . w)

J ~Jx

=E i s 3 } fYwJ

sinn(t-pAt) n(t - p At)

3 .3 .2. Effect of sampling

sin n(t' - p At) sin n(t" - p At) x n(t , - p At) 7c(t" - p At) x exp[ - j(ut + ut' + wt"] dtdt' dt" .

(35)

We verify that the results of Brillinger and Rosenblatt on the sampling of stationary signals [1] apply to the signal z(t) generated as in Section



J . Le Roux el al. / Signal Processing 36 (1994) 375-390

sin n(k'

387

- p At) sin n(k" - p At)

X

n(k' - p At)

n(k" - p

At)

x exp [ - j (uk + vk' + wk")] .

(42)

The Fourier transform of the sampled signal sin n(k - p At)/n(k - p At) is the periodization of the function equal to exp(jupAt) in -n < u < it . Consequently, all the terms of R(u, v, w) arc periodic and R(u, v, w) is the periodization in u, v and w of the function equal to R5(u,v,w)=E{s3}

I

exp[j(u + r + w)p At, (43)

2n-27%

A

2,cjAt - 2.

s

Fig. 9 . Third-order spectrum (Ra (u, u, w)) support of the analog band-Limited signal z(t) generated in Section 3 .3 .1 . This support consists in an hexagon and two triangles . After sampling, the support is the 3-D periodization of this figure . When z(t) is stationary (At < 2/3) . the support reduces to the hexagon in the plane u + v + w = 0 and R,(u, v, w) is zero in the two triangles u + r + w = ± Zit/At . Then the support of B(u, r) is the projection of the three-dimensional support of R ,(u, v. w) .

3 .3 .1 : Consider z(k) obtained by sampling z(t) at rate 1 . Its third-order covariance is E{z(k)z(k')z(k")} sin n(k-pAt)sinn(k'-pAt)

=E1s3}

r=-m n(k-pAt)

a(k'-pAt)

sin n(k" - p At) x n(k"

pAt)

Its Fourier transform is R(u, v, w) sin n(k = E1s3} y = -ck

-- `n k'--~

k"=

a

- in At)

n(k - p At)

in the cube - n < u, v, w < ic. If At < 2/3, R a (u, v, w) is non-zero only in u + v + w = 0 . In that case, z(k) is a stationary discrete sequence and its TOS, B(u, v) is the projection of R(u, v, w) on the plane w = 0 (Fig . 9) . It is the periodization of the function constant in the hexagon ACDEGH and equal to zero in the triangles ABC and EFG : B(u, v) = B a (u, for -

v) = R 5 (u,

n < u, r < n .

v, - u - r) (44)

was computed for At =1/2 . The corresponding B(u, v) is given in Fig . 10 . Note that when B(u, v) is estimated on one record of the signal, it is not close to zero in ABC, as it is factorizable by construction . B(u, v) decreases significantly in ABC only for an average on several records (here 100000 records of 128 samples each) . EXAMPLE . z(k)

This method for generating the sampling of a third-order stationary signal is an illustration of the result of Brillinger and Rosenblatt on the effect of sampling [1] . However, when At = I, z(k) is also a third-order stationary discrete sequence, although it is the sampling of a cyclo-stationary (and not stationary) band-limited analog signal . Then the correct interpretation of the TOS is the one given in Section 3 .2 .

388

.1 . Le Roux et al . / Signal Processing 36 (1994, 375 390

H

Fig . 10. Third-order spectrum modulus of the sampling of a third-order stationary analog signal y(n) generated as in (32) for At = 1/2 (average on 100000 records of 128 samples each) . Two representations are given simultaneously : a 3-D perspective view and a projection view in the plane (u, v) represented in contours . The amplitude is around 18-22 dB in ACDEGH and - 5 to - 15 dB in ABC and EFG.

sampled system output is factorizable ; this is generally not compatible with the assumption on the initial analog signal stationarity . So, when applying digital modelling techniques based on the factorizability of the TOS (assuming that the sampled signal is stationary), it is necessary to check whether the analog signal is stationary or not- The detection of significant values in ABC is an indication of non-stationarity [5] ; however, there are non-stationary signals whose bispectra, computed as averages on several experiments, vanish in ABC . (2) If the analog signal is stationary, there are several alternatives . Among them the following ones may be considered : - One can sample at more than three times the largest frequency of the band-limited signal; then both the analog and the sampled signal are factorizable simultaneously and one may apply algorithms developed in the case of deterministic signals HOS ; this may yield numerical difficulties as the TOS vanishes in 66% of the frequency domain . - When sampling the signal at twice the largest frequency, B(u, v) vanishes in ABC . The identification problem (22) is replaced by the factorization of B(u, v) in B(u, v) = G(u)G(v)G( - u - r)H(u, v),

4 . Remarks concerning the application of digital HOS techniques to sampled random signals (1) Several identification algorithms based on HOS suppose that the higher-order spectrum of the

where H(u, v) is the periodization of the function equal to one in ACDEGH and equal to zero in ABC and EFG_ Identification techniques have been developed in the frequency domain [8] . This problem needs to be

Table 2 Summary of the different kinds of signals in system identification based on HOS Type of signal Deterministic

Random

Sampled

Sampling of analog cyclostationary or digital i .i.d . sequence (Section 3 .2)

(Section 2)

FactorizableinB(u,r)=F(u)F(c)F(-u-u)

(45)

Sampling of analog stationary (Section 3 331 Btu, rJ=F(u)F(v)F(-u-v)H(u,cl

non-factorizable except when F„(u) = 0 for Jug 3 2n/3



J. Le Roux er al. 7 Signal Processing 36 (1994) 375-390

studied in the time domain . Then it can be stated as follows : Find g(k) such that r(m, n)

_ Y

p. LI q

g(k)g(k + p)g(k + q)

k

x h (m - p, n -

q),

I (46)

H(u, v):

for m # 0 :

h(0,0)=3/4,

h(m, 0) = h(0, m) = h(m, m) I-(-I)2 >< z m z

for mn(m - n)

540:

density of often known, it is generally simpler to use

IB(u,v)1'=1F(u)h1F(v)h1F(u+v)1 2.

where r(m, n) is the measured TOC and h(m, n) is the inverse Fourier transform of for m=n=0 :

389

h(m, n) _

( - U"n - ( - I)"`m + ( - I)'°-"(m - n) 2n 2 nm(m - n) (47) However, when one assumes the factorizability of the sampled system impulse response TOS as in (22), it is possible to reconstruct the data in ABC and EFG before applying identification algorithms developed in the frequency domain ; the value of the phase of B(u, v) in ABC and EFG can be deduced from the values of B(u, v) in ACDEGH through an extension of a formula proposed by Marron, Sanchez and Sullivan for solving the problem of phase unwrapping [7] . If 13(u, v) is factorizable,

( 49)

(3) If the input signal applied to the analyzed system is controlled (for instance in communication applications), one may use a non-stationary signal made of random independent pulses spaced at the sampling period and apply directly the algorithms developed for deterministic signals on the sampled output . Note that this third case is often the only one considered in HOS reconstruction algorithms . So, many of the proposed algorithms can be applied correctly only if one controls the input of the analyzed system or at least the phase of the input sampling clock . (4) In order to simulate a signal corresponding to the sampling of a band-limited third-order stationary input, one may use the random signal proposed in (32) for At < 2/3, for example At = 1/2 :

e(k)

= Y

s(p/2) s n(k(k

p/2) ,

(50)

where s(pAt) is a non-Gaussian sequence of i.i .d samples . So, the input will be white at the second-order and the TOS of the analyzed signal will (theoretically) vanish in ABC and EFG, after averaging on several experiments . Note that the quality of the pseudo-random generator producing s(p At) may be questionable and require improvements.

phase B(u, v) n u)B(v, = phase B(u, - n - u- v) B(u+v,it-u-v)

(48)

Through this formula, a value of phase B(u . v) in ABC is deduced from two values on the edge of the hexagon, phase B(u, n - u) on AC and phase B(u + v, n - u - v) on EG and one inside ACDEGH, phase B(v, n - u - v), that is phase B,(v, n - u - v) . The estimation of phase Btu, n - u) raises theoretical and practical difficulties . Similar formulas can be developed in the case of the modulus . However, as the signal spectral

5 . Conclusion In this paper, we have attempted to illustrate through examples the important difference between the TOS of deterministic and random stationary signals after sampling that was considered theoretically by Brillinger, Rosenblatt and Hinich . We have shown that this difference disappears if the sampling rate is three times larger than the bandwidth of the analog signal . We have also considered the problem of random signal simulation in HOS analysis, and proposed



3 90

.1 . Le Roux er al . / Signal Processing 36 (1994) 375- 390

a method for generating a digital signal with constant second-order spectrum that is equivalent to the sampling of a third-order stationary band-limited signal . Finally, we have insisted on the care needed when applying digital identification algorithms to actual random signals when the sampling frequency is less than three times the signal bandwidth . References [1] D. Brillinger and M . Rosenblatt, "Computation and interpretation of kth-order spectra" . in: B . Harris, ed . . Spectral Analysis of Time Series, Wiley. New York, 1967, pp . 189-232 . [2] D .E . Dudgeon and R .M . Mcrsereau, Multidimensional Signal Processing, Prentice-Hall, Englewood Cliffs, NJ, 1984 . [3] G .B . Giannakis . "Signal reconstruction from multiple correlations : Frequency and time domain approaches", J . Opt . Soc. Amer . A, Vol . 6, No. 5, May 1989, pp. 682-697 . [4] G .B . Giannakis and J.M . Mendel, "Identification of nonminimum phase systems using higher order statistics", IEEE Trans . Acoust . Speech Signal Process . . Vol . 37, No. 3, March 1989, pp . 360-377. [5] M .J. Hinich, "Detecting a transient signal by bispectral analysis", IEEE Trans. Acoust. Speech Signal Process ., Vol . 38, No. 7, July 1990, pp. 1277 1283 . [6] J . Le Roux, "Reconstruction of a sampled signal Fourier transform from its bispectrum", Internal, Signal Processing

Workshop on Higher Order Statistics . Chamrousse (France), 10-12 July 1991, pp . 233-236 ; "Bispectral analysis and reconstruction in the frequency domain of mono and bidimensional deterministic sampled signals" . Multidimensional Systems and Signal Processing, No . 4, 1993, pp . 39-66 . [7] J.C . Marron, P .P . Sanchez and R .C . Sullivan, "Unwrapping algorithm for least-squares phase recovery from the modulo 2n bispectrum phase", J. Opt. Soc. Amer ., Vol . 7, No . 1 . January 1990, pp. 14-20. [8] T. Matsuoka and T .J . Ulrych, "Phase estimation using the bispectrum", Proc . IEEE, Vol . 72, No . 10, October 1984, pp. 1403-1411 . [9] J.M . Mendel, "Tutorial on higher-order statistics (spectra) in signal processing and system theory : Theoretical results and applications" . Proc . IEEE, Vol . 79, No. 3. March 1991, pp . 278 305 . [10] C.L . Nikias and M . Raghuveer, "Bispectrum estimation : A digital signal processing framework", Proc . IEEE Vol . 75, 1987, pp . 869-891 . [11] R . Pan and C .L . Nikias, "The complex ccpstrum of higherorder moments", IEEE Trans . Acoust . Speech Signal Process ., Vol . 36, 1988, pp. 186-205. [12] M . Rangoussi and G . Giannakis, "FIR modeling using log-bispectra : Weighted least squares algorithms and performance analysis", IEEE Trans. Circuits Systems, March 1991, pp. 281-296 . [13] A .M . Tekalp and A .T. Erdem, "Higher order spectrum factorization in one and two dimensions with applications in signal modelling and nonminimum phase system identification". IEEE Trans Acoust . Speech Signal Process ., Vol . 37, No. 10, October 1989 . pp. 1537-1549 .