Higher-Order Spectral Analysis




Definitions and Properties of HOS


HOS Computation from Real Data Indirect Method • Direct Method


Linear Processes

Non-Parametric Methods • Parametric Methods

Athina P. Petropulu


Nonlinear Processes

Drexeо University


HOS in Biomedical Signal Processing


The past 20 years witnessed an expansion of power spectrum estimation techniques, which have proved essential in many applications, such as communications, sonar, radar, speech/image processing, geophys­ics, and biomedical signal processing [32], [27], [20]. In power spectrum estimation, the process under consideration is treated as a superposition of statistically uncorrelated harmonic components. The dis­tribution of power among these frequency components is the power spectrum. As such phase relations between frequency components are suppressed. The information in the power spectrum is essentially present in the autocorrelation sequence, which would suffice for the complete statistical description of a Gaussian process of known mean. However, there are applications where one would need to obtain information regarding deviations from the Gaussianity assumption and presence of nonlinearities. In these cases power spectrum is of little help, and a look beyond the power spectrum or autocorrelation domain is needed. Higher-Order Spectra (HOS) (of order greater than 2), which are defined in terms of higher-order cumulants of the data, do contain such information [35]. The third order spectrum is commonly referred to as bispectrum and the fourth-order as trispectrum. The power spectrum is also a member of the higher-order spectral class; it is the second-order spectrum.

HOS consist of higher-order moment spectra, which are defined for deterministic signals, and cumu­lant spectra, which are defined for random processes. In general there are three motivations behind the use of HOS in signal processing: (1) to suppress Gaussian noise of unknown mean and variance, (2) to reconstruct the phase as well as the magnitude response of signals or systems, and (3) to detect and characterize nonlinearities in the data.

The first motivation stems from the property of Gaussian processes to have zero higher-order spectra of order greater than 2). Due to this property, HOS are high signal-to-noise ratio domains in which one can perform detection, parameter estimation, or even signal reconstruction even if the time domain noise is spatially correlated. The same property of cumulant spectra can provide a means of detecting and characterizing deviations of the data from the Gaussian model.

The second motivation is based on the ability of cumulant spectra to preserve the Fourier phase of signals. In the modeling of time series, second-order statistics (autocorrelation) have been heavily used because they are the result of least-squares optimization criteria. However, an accurate phase reconstruc­tion in the autocorrelation domain can be achieved only if the signal is at minimum phase. Nonminimum phase signal reconstruction can be achieved only in the HOS domain, due to the HOS ability to preserve the phase.

Being nonlinear functions of the data, HOS are quite natural tools in the analysis of nonlinear systems operating under a random input. General relations for stationary random data passing through an arbitrary linear system exist and have been studied extensively. Such expressions, however, are not available for nonlinear systems, where each type of nonlinearity must be studied separately. Higher-order correlations between input and output can detect and characterize certain nonlinearities [56], and for this purpose several higher-order spectra based methods have been developed.

This chapter is organized as follows. In Section 57.2 definitions and properties of cumulants and higher-order spectra are introduced. In section 57.3 methods for estimation of HOS from finite length data are presented and the asymptotic statistics of the corresponding estimates are provided. In Section 57.4 methods for identification of a linear time-invariant system from HOS of the system output are outlined, while in Section 57.6, non-linear system estimation is considered. Section 57.6 summarizes the research activity related to HOS as applied to biomedical signal analysis and provides details on a particular application, namely the improvement of resolution of medical ultrasound images via HOS is discussed.

Definitions and Properties of HOS

In this chapter only random one-dimensional processes are considered. The definitions can be easily extended to the two-dimensional case [34].

The join moments of order r of the random variables x1, …, xn are given by [42]:

Higher-Order Spectral Analysis


Higher-Order Spectral Analysis


Where k1 + …+ kn = r, and O(-) is their joint characteristic function. The corresponding joint cumulants are defined as:

Higher-Order Spectral Analysis

□ □ 0 •

подпись: □ □ 0 •

Cumxk1, □, xn

подпись: cumxk1, □, xn(57.2)

For a stationary discrete time random process x(k), (k denotes discrete time), the moments of order n are given by:

Higher-Order Spectral Analysis


Where Ј{.} denotes expectation. The nth order cumulants are functions of the moments of order up to

N, i. e.,

1st order cumulants:

Higher-Order Spectral Analysis

3rd order cumulants:

CX[,^□□ mX[,^[[x[№№mX[№ mX[ — T№2[X[I (57.6)

4th order cumulants:

CX[,q, qj□□ m4X[,^,^[mX[QnX[ № mX[[<[ — q□

□ mX[[[ — q[

MX №X [-^, №-T № № [, № № mX [, q, [ mX [, q, [ (57.7)

[X[fe[№ m2X[[ mX[[ mX[ — q[ m2X[ □ T, [

□mX[ -q[6[x[

The general relationship between cumulants and moments can be found in [35].

Some important properties of moments and cumulants are summarized next.

[P1] If X(k) is a jointly Gaussian process, then cXn (t1, t2, …, Tn-1) = 0 for n > 2. In other words, all the information about a Gaussian process is contained in its first and second-order cumulants. This property can be used to suppress Gaussian noise, or as a measure for non-Gaussianity in time series.

[P2] If X(k) is symmetrically distributed, then cX3 (t1, t2) = 0. Third-order cumulants suppress not only Gaussian processes, but also all symmetrically distributed processes, such as uniform, Laplace, and Bernoulli-Gaussian.

[P3] Additivity holds for cumulants. If X(k) = s(k) + w(k), where s(k), w(k) are stationary and statistically independent random processes, then cXn (t1, t2, …, Tn-1) = c5n (t1, t2, …, Tn-1) + cwn (t1, t2, …, Tn-1). It is important to note that additivity does not hold for moments.

If w(k) is a Gaussian representing noise which corrupts the signal of interest, s(k), then by means of (P2) and (P3), we conclude that cXn (t1, t2, …, Tn-1) = c5n (t1, t2, …, Tn-1), for n > 2. In other words, in higher-order cumulant domains the signal of interest propagates noise free. Property (P3) can also provide a measure of the statistical dependence of two processes.

[P4] if X(k) has zero mean, then cXn (t1, …, Tn-1) = mXn (t1, …, Tn-1), for n < 3.

Higher-order spectra are defined in terms of either cumulants (e. g., cumulant spectra) or moments (e. g., moment spectra).

Assuming that the nth order cumulant sequence is absolutely summable, the nth order cumulant spectrum of X(k), CXn (ro1, ro2, ron-1), exists, and is defined to be the (n-1)-dimensional Fourier

Transform of the nth order cumulant sequence. In general, CXn (ro1, ro2, ron-1) is complex, i. e., it has

Magnitude and phase. In an analogous manner, moment spectrum is the multidimensional Fourier transform of the moment sequence.

If v(k) is a stationary non-Gaussian process with zero mean and nth order cumulant sequence

Where S(.) is the delta function, v(k) is said to be nth order white. Its nth order cumulant spectrum is then flat and equal to yxn.

Cumulant spectra are more useful in processing random signals than moment spectra since they posses properties that the moment spectra do not share, i. e., properties P1, P3, and flatness of spectra corre­sponding to higher-order white noise.

HOS Computation from Real Data

The definitions of cumulants presented in the previous section are based on expectation operations, and they assume infinite length data. In practice we always deal with data of finite length, therefore, the cumulants can only be estimated. Two methods for cumulants and spectra estimation are presented next for the third-order case.

Indirect Method:

Let x(k), k = 0,., N — 1 be the available data.

Segment the data into K records of M samples each. Let x'(k), k = 0, …, M — 1, represent the ith record.

Subtract the mean of each record.

Estimate the moments of each segment x'(k) as follows:

1 i2

M^[, T2)□ — □ x! Wx'[I[]

M i□ 4 (57.9)

L1 □ max[, □ □!, □ □2 [ i2 □ min(M □ 1□ □!, M □ 1□ □2, M □ l[ [ □ L, |d^| □ L, i □ 1, 2, □ , K.

Since each segment has zero mean, its third-order moments and cumulants are identical, i. e., cx3 (T!, t2) = mx3(T!, t2).

Compute the average cumulants as:


Cxx [, T2 □□ — □ m3 [T1, T2 □ (57.10)


Obtain the third-order spectrum (bispectrum) estimate as the two-dimensional Discrete Fourier Transform of size Mx M of the windowed c| (tj, t2), i. e.,

L L,□ 2 □ 2 □ □

C3x(k[)[] ^c3x[Tl, T2MkA0w[Tl, T2□ k1, k2 □ 0, □, M□ 1 (57.11)

Tj □□l T2 □□l

Where L < M-1, and w(t1, t2) is a two-dimensional window of bounded support, introduced to smooth out edge effects. The bandwidth of the final bispectrum estimate is A = 1/L.

A complete description of appropriate windows that can be used in Eq. 57.11 and their properties can be found in [35]. A good choice of cumulant window is:




подпись: l
N |sin -nql+Di □y tos f

Which is known as the minimum bispectrum bias supremum [36].

Direct Method

Let X(k), k = 1, … N be the available data.

Segment the data into K records of M samples each. Let X'(k), k = 0, M — 1, represent the ith record.

Subtract the mean of each record.

Compute the Discrete Fourier Transform X'(k) of each segment, based on M points, i. e.,


подпись: (57.14)X1 [[□□ x1 [Q°J^, k □ 0,1, □ , M□ 1, i □ 1,2, □ , K.

N □ 0

The third-order spectrum of each segment is obtained as

C3Xi[, k2[-1Xi[[[ [X* [ □ K[ i □ 1, □ , K. (57.15)

Due to the bispectrum symmetry properties, CX3i(kpk2) needs to be computed only in the triangular region 0 < k2 < kj, kj +k2 < M / 2.

In order to reduce the variance of the estimate, additional smoothing over a rectangular window of size (M3 x M3) can be performed around each frequency, assuming that the third order spectrum is smooth enough, i. e.,

M3/2Q M3/2Q

^3Xi [, K [□^ □ □ C3Xi [□ [, k2 □ n2 [] (57.16)

M3 n □ □M3/2 «2 □ □M3/2

Finally, the third-order spectrum is given as the average over all third-order spectra, i. e.,

C"X Qik[□1 [] CXi [,k2[ (57.17)

K o

The final bandwidth of this bispectrum estimate is A = M3/M, which is the spacing between frequency samples in the bispectrum domain.

For large N, and as long as

Both the direct and the indirect methods produce asymptotically unbiased and consistent bispectrum estimates, with real and imaginary part variances [50]:

Var|=Re[[3x (, k2) □ varj=jlm[[3x D, k

Rr2C2xD)(k2DC2x[k, □ k2D indirect (57-19)

1 cxDk DDk DDt n k r= 2 1 2 2 2 1 2

#2N C2 ^ 1(2)2 Lki n k2D U. qDt,)[>2Dp;D nk2D direct,

Where V is the energy of the bispectrum window.

From the above expressions, it becomes apparent that the bispectrum estimate variance can be reduced by increasing the number of records, or reducing the size of the region of support of the window in the cumulant domain (L), or increasing the size of the frequency smoothing window (M3), etc. The relation between the parameters M, K, L, M3 should be such that 57.18 is satisfied.

Linear Processes

Let x(k) be generated by exiting a linear, time-invariant (LTI), exponentially stable, mixed-phase system with frequency response H(ro) (or impulse response h(k)), with an nth order white, zero-mean, non — Gaussian stationary process v(k), i. e.,

X Dk D= vDf ph Dk D (57.20)

Where * denotes convolution.

The output Hth order cumulant equals [10]:

C;D,□,□□□:□ hD®nqD-hin33 (57.21)


Where yn is a scalar constant and equals the nth order spectrum of v(k). The output nth order cumulant spectrum is given by

C;D □ 2,□ , □ „□Dн>І1D■■H>}naDh*D n — +□ nit] (57.22)

For a linear non-Gaussian random process x(k), the nth order spectrum can be factorized as in Eq. 57.22 for every order n, while for a nonlinear process such a factorization might be valid for some orders only (it is always valid for n = 2).

It can be easily shown that the cumulant spectra of successive orders are related as follows:

CxD □ 2,□ ,0D= CX-1 D □ 2, □ , □ n-2. (57.23)

For example, the power spectrum of a non-Gaussian linear process can be reconstructed from the bispectrum up to a scalar factor, i. e.,

N v

While for n = 2, Eq. 57.22 provides the magnitude only of H(ro), as a function of CXn(m 1,ro2, …, ron-1), for n > 2 both the magnitude and phase of H(m) can be obtained. A significant amount of research has been devoted to estimating the system phase from the phase of the bispectrum or trispectrum [33], [35],

, and [47].


Y&Q x[ Q [ (57.25)

Where X(k) was defined in Eq. 57.20, and w(k) is zero-mean Gaussian noise, uncorrelated to X(k). Methods for reconstructing h(k) from y(k) can be divided into two main categories: non-parametric and para­metric.

Non-Parametric Methods

Non-parametric methods reconstruct the system by recovering its Fourier phase and magnitude. Some utilize the whole bispectrum (or polyspectrum) information [33], [7], [43], [53], [41], [30], and others use fixed, one-dimensional polyspectrum slices [31], [13]. In [45] and [49] methods for system recovery from arbitrarily selected HOS slices have been proposed. A system method based on preselected polyspec­trum slices, allows regions where polyspectrum estimates exhibit high variance, or regions where the ideal polyspectrum is expected to be zero, such as in the case of bandlimited systems to be avoided. Assuming that there existed a mechanism to select “good” slices, it can also lead to computational savings.

Let the bispectrum of y(k) in Eq. 57.25, i. e., CX3(m 1,ro2), be discretized on a square grid of size 2n/N x

N/N, and let {Cy3(m, l), m = 0, …, N-1}, {Cy3(m, l + r), m = 0,…, N-1} be two bispectral slices. In [49] it was shown that unique identification of the system impulse response can be performed base on two horizontal slices of the discretized third-order spectrum of the system output, as long as the distance between slices, i. e., r and N are co-prime numbers. Let

H □[log[[□ ,logH[□ 1[[□ 1[ 1[ (57.26)

Where H(k) is the Discrete Fourier Transform of h(n). Then h can be obtained as the solution of the system:

All □ b, (57.27)


BO log C>h[k □ r □ l, l№ log C>h [k □ r □ l, l □ rQ q

K □ 0, □, N□ 1 (57.28)

NO 1

Ci, r □ N □[ [ [, l □ r0 log Cy [, lQ

K □ 0

And A is a [(N — 1) x (N — 1)] sparse matrix; its first row contains a 1 at the rth column and 0s elsewhere. The kth row contains a -1 at column (k — 1)moduloN, and a 1 at column (k + r — 1)moduloN. Although b(k) requires the computation of the logarithm of a complex quantity, it was shown in [49], that the log — bispectrum can be computed based on the principal argument of the bispectral phase, thus bypassing the need for two-dimensional phase unwrapping. The bispectrum needed in this method can be substi­tuted with a bispectral estimate, obtained with any of the methods outlined in section 57.3. Reduction of the HOS estimate variance can be achieved via a low-rank approach [4], [9], and [8].

The above described reconstruction method can be easily extended to any order spectrum [49]. Matlab code for this method can be found at Http://www. ece. drexel. edu/CSPL.

Parametric Methods

These methods fit a parametric model to the data x(k). For example, let us fit a real autoregressive moving average (ARMA) model to x(k), i. e.,

P q

□ aOUXUkf □ i bk’Lk □ j) (57.29)

I= 0 j □ 0

Where v(k) is a stationary zero-mean nth order white non-Gaussian process, and a(i) and b(j) represent the AR and MA parameters. Let y(k) be as defined in Eq. 57.25. For estimation of the AR parameters

From y(k), equations analogous to the Yule-Walker equations can be derived based on third-order

Cumulants of y(k), i. e.,


□ ak^i, j[0, □□ q, (57.30)

I □ 0



Ak)y (□ i, j)n k j[] □□ q, (57.31)

I ni

Where it was assumed a(0) = 1. Concatenating Eq. 57.31 for t = q + 1, …, q + M, M > 0, and j = q — p, …, q, the matrix equation

Ca □ c, (57.32)

Can be formed, where C and c are a matrix and a vector, respectively, formed by third-order cumulants of the process according to Eq. 57.31, and the vector a contains the AR parameters. If the AR order p is unknown and Eq. 57.32 is formed based on an overestimate of p, the resulting matrix C always has rank p. In this case the AR parameters can be obtained using a low-rank approximation of C [16].

Using the estimated AR parameters, a(i), i = 1,., p, a, pth order filter with transfer function A(z) = 1 + Zp^ a(i’)z-1 can be constructed. Based on the filtered-through A(z) process y(k), i. e., y(k), or otherwise known as the residual time series [16], the MA parameters can be estimated via any MA method [34], for example:

C 3y k, k □

BkLn-^—n k □ 0,1, □ , q (57.33)

C3 Lk,0L

Known as the c(q, k) formula [15].

Practical problems associated with the described approach are sensitivity to model order mismatch, and AR estimation errors that propagate in the estimation of the MA parameters. Additional parametric

Methods can be found in [5], [14], [34], [35], [44]. [63], [64].

Higher-Order Spectral Analysis

FIGURE 57.1 Second-order Volterra system. Linear and quadratic parts are connected in parallel.

Nonlinear Processes

Despite the fact that progress has been established in developing the theoretical properties of nonlinear models, only a few statistical methods exist for detection and characterization of nonlinearities from a finite set of observations. In this section we will consider nonlinear Volterra systems excited by Gaussian stationary inputs. Let y(k) be the response of a discrete time invariant, pth order, Volterra filter whose input is X(k). Then,

YQQ h □□ □ h [, □, qQQ □dQ-xK □□i [] (57.34)

I □, Ci

Where h^Tp…, t) are the Volterra kernels of the system, which are symmetric functions of their argu­ments; for causal systems h^Tp… t) = 0 for any t; < 0.

The output of a second-order Volterra system when the input is zero-mean stationary is

Y&Q K □□ hr 0&□^[□□ □ h[,q□ qQ[□^0 (57.35)

□1 ^2

Equation 57.35 can be viewed as a parallel connection of a linear system K^tJ and a quadratic system h^Tp t2) as illustrated in Fig. 57.1. Let

C2y [0 { □ #[[ my ] (57.36)

Be the cross-covariance of input and output, and

CX^r^Q □QOS □ □, Q& Q my] (57.37)

Be the third-order cross-cumulant sequence of input and output. The system’s linear part can be obtained by [61]

C2y [[

f^n, (57.38)

1 CX [[

N n C3Xxy UP1, □ 2)

H k □□2 □□ ^2 [ (57.39)

Where CXf(ro) and Cxf(m 1, ro2) are the Fourier transforms of cy (t) and tf’y(t1, t2), respectively. It should be noted that the above equations are valid only for Gaussian input signals. More general results assuming non-Gaussian input have been obtained in [22] and [48]. Additional results on particular nonlinear systems have been reported in [11] and [54].

An interesting phenomenon, caused by a second-order nonlinearity, is the quadratic-phase coupling. There are situations where nonlinear interaction between two harmonic components of a process con­tribute to the power of the sum and/or difference frequencies. The signal

X()= AcosQ 1k □*1Qn BcosQ 2k + *2) (57.40)

After passing through the quadratic system:

XQB Јx2 (Be + 0, (57.41)

Contains cosinusoidal terms in (^1, 01), (^2,02), (2^1, 201), (2^2, 202), (^1 + ^2, 01 + 02), (^1 — ^2, 01 — 02). Such a phenomenon that results in phase relations that are the same as the frequency relations is called quadratic phase coupling [28]. Quadratic phase coupling can arise only among harmonically related

Components. Three frequencies are harmonically related when one of them is the sum or difference of

The other two. Sometimes it is important to find out if peaks at harmonically related positions in the power spectrum are in fact phase coupled. Due to phase suppression, the power spectrum is unable to provide an answer to this problem. As an example, consider the process [51]


Xk[ □ cosk k□, ;[ (57.42)


Where ^1 > ^2 > 0, + ^5 > 0, ^3 = ^1 + ^2, + ^5, ^,… ^5 are all independent, uniformly

Distributed random variables over (0, 2n), and + ^5. Among the six frequencies (^1, ^2, ^3) and

(^4, ^5, ^6) are harmonically related, however, only is the result of phase coupling between ^4 and ^5. The power spectrum of this process consists of six impulses at Xi, i, …, 6 (Fig. 57.2), offering no indication whether each frequency component is independent or a result of frequency coupling. On the other hand, the bispectrum of x(k), (evaluated in its principal region) is zero everywhere, except at point (^4, ^5) of the (ro 1, ro2) plane, where it exhibits an impulse (Fig. 57.2b). The peak indicates that only ^4, ^5 are phase coupled.

The biocherence index, defined as

CXX (, □ 2 B


подпись: x]1,ь 2 LB 1 , (57.43)

Has been extensively used in practical situations for the detection and quantification of quadratic phase coupling [28]. The value of the bicoherence index at each frequency pair indicates the degree of coupling

Higher-Order Spectral Analysis

FIGURE 57.2 Quadratic phase coupling. (a) The power spectrum of the process described in Eq. 57.42 cannot determine what frequencies are coupled. (b) The corresponding magnitude bispectrum is zero everywhere in the principle region, except at points corresponding to phase coupled frequencies.

Among the frequencies of that pair. Almost all bispectral estimators can be used in Eq. 57.43. However, estimates obtained based on parametric modeling of the bispectrum have been shown to yield resolution superior [51 and 52] to the ones obtained with conventional methods.

HOS in Biomedical Signal Processing

The applications of HOS on biomedical signals are clustered according to the HOS property they most rely on, i. e., (1) the ability to describe non-Gaussian processes and preserve phase, (2) the Gaussian noise immunity, and (3) the ability to characterize nonlinearities.

In the first class are the works of [1], [3], and [60], where ultrasound imaging distortions are estimated from the ultrasound echo and subsequently compensated for, to improve the diagnostic quality of the image. HOS have been used in modeling the ultrasound radio-frequency (RF) echo [2], where schemes for estimation of resolvable periodicity as well as correlations among nonresolvable scatters have been proposed. The “tissue color,” a quantity that describes the scatterer spatial correlations, and which can be obtained from the HOS of the RF echo, has been proposed [2] as a tissue characterization feature.

The skewness and kurtosis of mammogram images have been proposed in [18] as a tool for detecting microcalcifications in mammograms.

In the second class are methods that process multicomponent biomedical signals, treating one com­ponent as the signal of interest and the rest as noise. HOS have been used in [19] to process lung sounds in order to suppress sound originating from the heart, and in [65] to detect human afferent nerve signals, which are usually observed in very poor signal-to-noise ratio conditions.

Most of the HOS applications in the biomedical area belong to the third class, and usually investigate the formation and the nature of frequencies present in signals through the presence of quadratic phase coupling (QPC). The bispectrum has been applied in EEG signals of the rat during various vigilance states

, where QPC between specific frequencies was observed. QPC changes in auditory evoked potentials of healthy subjects and subjects with Alzheimer’s dementia has been reported [55]. Visual evoked potentials have been analyzed via the bispectrum in [59], [24], and [21]. Bispectral analysis of interactions between electrocerebral activity resulting from stimulation of the left and right visual fields revealed nonlinear interactions between visual fields [57]. The bispectrum has been used in the analysis of electromyogram recruitment [66], and in defining the pattern of summation of muscle fiber twitches in the surface mechanomyogram generation [40]. QPC was also observed in the EEG of humans [6] and [23].

In the sequel, presented in some detail are the application of HOS on improving the resolution of ultrasound images, a topic that has recently attracted a lot of attention.

Ultrasonic imaging is a valuable tool in the diagnosis of tumors of soft tissues. Some of the distinctive features of ultrasound are its safe, non-ionizing acoustic radiation, and its wide availability as a low cost, portable equipment. The major drawback that limits the use of ultrasound images in certain cases, (e. g., breast imaging) is poor resolution. In B-Scan images, the resolution is compromised due to: (1) the finite bandwidth of the ultrasonic transducer, (2) the non-negligible beam width, and (3) phase aberrations and velocity variations arising from acoustic inhomogeneity of tissues themselves. The observed ultra­sound image can be considered as a distorted version of the true tissue information. Along the axial direction the distortion is dominated by the pulse-echo wavelet of the imaging system, while along the lateral direction the distortion is mainly due to finite-width lateral beam profile. Provided that these distortions are known in advance, or non-invasively measurable, their effects can be compensated for in the observed image. With propagation in tissue, however, both axial and lateral distortions change due to the inhomogeneities of the media and the geometry of the imaging system. They also change among different tissue types and individuals. Distortions measure in a simplified setting, e. g. in a water tank, are rarely applicable to clinical images, due to the effects of tissue-dependent components. Therefore, distortions must be estimated based on the backscattered RF data that lead to the B-mode image.

Assuming a narrow ultrasound beam, linear propagation and weak scattering, the ultrasonic RF echo, yt(n), corresponding to the tth axial line in the B-mode image is modeled as [1]:

YtO h(Qfcf([+ WjOz’n 1,2, □, (57.44)

Where n is the discrete time; wt(n) is observation noise; f(n) represents the underlying tissue structure and is referred to as tissue response; and ht(n) represents the axial distortion kernel. Let us assume that: (A1) yt(n) is non-Gaussian; (A2) f(n) is white, non-Gaussian random process; (A3) wt(n) is Gaussian noise uncorrelated with f(n); and (A4) yt(n) is a short enough segment, so that the attenuation effects stay constant over its duration. (Long segments of data can be analyzed by breaking them into a sequence of short segments.) A similar model can be assumed to hold in the lateral direction of the RF image.

Even though the distortion kernels were known to be non-minimum phase [26], their estimation was mostly carried out using second-order statistics (SOS), such as autocorrelation or power spectrum, thereby neglecting Fourier phase information. Phase is important in preserving edges and boundaries in images. It is particularly important in medical ultrasound images where the nature of edges of a tumor provide important diagnostic information about the malignancy. Complex cepstrum-based operations that take phase information into account have been used [60] to estimate distortions. HOS retain phase and in addition, are not as sensitive to singularities as the cepstrum. HOS were used [1] for the first time to estimate imaging distortions from B-scan images. It was demonstrated that the HOS-based distortion estimation and subsequent image deconvolution significantly improved resolution. For the case of breast data, in was demonstrated [3] that deconvolution via SOS-based distortion estimates was not as good as its HOS counterpart.

In the following we present some results of distortion estimation followed by deconvolution of clinical B-scan images of human breast data. The data were obtained using a flat linear array transducer with a nominal center frequency of 7.5 MHz on a clinical imaging system UltraMark-9 Advanced Technology Laboratories. Data were sampled at a rate of 20 MHz. Figure 57.3(A) showspartsof the originalimage, where the the logarithm of the envelope has been used for display purposes. Axial and lateral distortion kernels were estimated from the RF data via the HOS-based non-parametric method outlined in Section 57.4 of this chapter [49], and also via an SOS power cepstrum based method [3]. Each kernel estimate was obtained from a rectangular block of RF data described by (x, y, Nx, Ny) where (x, y) are the co-ordinates of the upper left corner of the block, Nx is its lateral width, and Ny is the axial height. Note that y corresponds to the depth of the upper left corner of the image block from the surface of the transducer. In the following, all dimensions are specified in sample numbers. The size of the images used was 192 x 1024. Assuming that within the same block the axial distortion kernel does not significantly depend on the lateral location of the RF data, all axial RF data in the block can be concatenated to form a longer one-dimensional vector. In both HOS and SOS based axial kernel estimations, it was assumed that Nx = 10, Ny = 128, N = 128. Note that the Ny = 128 samples correspond to 5mm in actual tissue space, as is reasonable to assume that attenuation over such a small distance may be assumed constant. Fifty kernels were estimated from the blocks (x, 400, 10, 128), the x taking values in the range 1.50. Lateral kernels were estimated from the same images with parameters Nx = 192, Ny =10, and N = 64. All the kernels were estimated from the blocks (1, y, 192, 10), the y taking values in the range 400.600 increasing each time by 5. Data from all adjacent lateral image lines in the block were concatenated to make a long one-dimensional vector. Note that Ny = 10 corresponds to a very small depth, 0.4mm in the real tissue space, hence he underlying assumption thatthelateraldistortion kernel doesnotvary much over this depth, shouldbe reasonable.

Figure 57.3(B-c) show that the result of lateral followed by axial deconvolution of the image of Fig. 57.3(A), using respectively HOS — and SOS-based distortion kernel estimates. The deconvolution was performed via the constrained Wiener Filter technique [62]. According to Fig. 57.3, tHe deconvolution in the RF-domain resulted in a significant reduction in speckle size. The speckle size appears smaller in the HOS-deconvolved image than in SOS-deconvolved ones, which is evidence of higher resolution. The resolution improvement can be quantified based on the width of the main lobe of the autocovariance of the image. The axial resolution gain for the HOS-based approach was 1.8 times that of the SOS-based approach. The lateral resolution gain for the HOS-based approach was 1.73 times as much. From a radiologist’s point of view, overall, the HOS image appears to have better spatial resolution than the original as well as the SOS deconvolved images. Conversely, the SOS image seems to have incurred a loss of both axial and lateral resolution.


Parts of this chapter have been based on the book: Nikias, C. L. and Petropulu, A. P., Higher Order Spectra Analysis: A Nonlinear Signal Processing Framework, Prentice Hall, Englewood Cliffs, NJ, 1993.

Major support for this work came from the National Science Foundation under grant MIP-9553227, the Whitaker Foundation, the National Institute of Health under grant 2P01CA52823-07A1, Drexel University and the Laboratoire des Signaux et Systems, CNRS, Universite Paris Sud, Ecole Superieure d’Electric, France.

The author would like to thank Drs. F. Forsberg and E. Conant for providing the ultrasound image and for evaluating the processed image.

Higher-Order Spectral Analysis

FIGURE 57.3 (a) Original image, (b) HOS-based deconvolution, (c) SOS-based deconvolution. In all cases the

Logarithm of the envelope is displayed. Autocovariance plots for (d) original (e) HOS deconvolved and (f) SOS deconvolved RF data.


Abeyratne, U., Petropulu, A. P., and Reid, J. M. “Higher-Order Spectra Based Deconvolution of Ultrasound Images,” IEEE Trans. Ultrason. Ferroelec., andFreq. Con., vol. 42(6): 1064-1075, Novem­ber 1995.

Abeyratne, U. R., Petropulu, A. P., and Reid, J. M. “On Modeling the Tissue Response from Ultrasound Images,” IEEE Trans. Med. Imag., vol. 14(4): 479-490, August 1996.

Abeyratne, U. R., Petropulu, A. P., Golas, T., Reid, J. M., Forsberg, F., and Consant, E. “Blind deconvolution of ultrasound breast images: A comparative study of autocorrelation based versus higher-order statistics based methods,” IEEE Trans. Ultrason. Ferroelec, and Freq. Con., vol. 44(6): 1409-1416, November 1997.

Andre, T. F., Nowak, R. D., Van Veen, B. D. “Low Rank Estimation of Higher Order Statistics,” IEEE Trans. on Sig. Proc., vol. 45(3): 673-685, March 1997.

Alshbelli, S. A., Venetsanopoulos, A. N., and Cetin, A. E. “Cumulant Based Identification Approaches for Nonminimum Phase Systems,” IEEE Trans. Sig. Proc., vol. 41(4): 1576-1588, April 1993.

Barnett, T. P., Johnson, L. C. et al., “Bispectrum analysis of electroencephalogram signals during waking and sleeping,” Science, 172: 401-402, 1971.

Bartlet, H., Lohmann, A. W., and Wirnitzer, B. “Phase and amplitude recovery from bispectra,” Applied Optics, 23(18): 3121-3129, Sept. 1984.

Bradaric, I., and Petropulu, A. P. “Subspace Design of Low Rank Estimators for Higher Order Statistics,” IEEE Trans. Sig. Proc., submitted in 1999.

Bradaric, I., and Petropulu, A. P. “Low Rank Approach in System Identification Using Higher Order Statistics,” 9th IEEE Sig. Processing Workshop on Statistical Signal and Array Processing — SSAP’98, Portland, Oregon, September 1998.

Brillinger, D. R. and Rosenblatt, M. “Computation and Interpretation of kth-order Spectra,” In: Spectral Analysis of Time Series, B Harris, Ed., John Wiley and Sons, New York, NY, 189-232, 1967.

Brillinger, D. R., “The Identification of a Particular Nonlinear Time Series System,” Biometrika, 64(3): 509-515, 1977.

Cohen, F. N. “Modeling of Ultrasound Speckle with Applications in Flaw Detection in Metals,” IEEE Trans. on Signal Processing, 40(3): 624-632, 1992.

Dianat, S. A. and Raghuveer, M. R. “Fast Algorithms for Phase and Magnitude Reconstruction from Bispectra,” Optical Engineering, 29: 504-512, 1990.

Fonollosa, J. A. R. and Vidal, J. “System Identification Using a Linear Combination of Cumulant Slices,” IEEE Trans. Sig. Proc., 41: 2405-2411, 1993.

Giannakis, G. B. “Cumulants: A Powerful Tool in Signal Processing,” Proc. IEEE, 75, 1987.

Giannakis, G. B. and Mendel, J. M., “Cumulant-Based Order Determination of Non-Gaussian ARMA Models,” IEEE Trans. on Acoustics, Speech, and Signal Processing, 38: 1411-1423, 1990.

Giannakis, G. B. and Swami, A. “On Estimating Noncausal Nonminimum Phase ARMA Models of Non-Gaussian Processes,” IEEE Transactions Acoustics, Speech, and Signal Processing, 38(3): 478-495, 1990.

Gurcan, M. N., Yardimci, Y., Cetin, A. E., Ansari, R. “Detection of Microcalcifications in Mammo­grams Using Higher-Order Statistics,” IEEE Sig. Proc. Letters., 4(8), 1997.

Hadjileontiadis, L. and Panas, S. M., “Adaptive Reduction of Heart Sounds from Lung Sounds Using Fourth-Order Statistics,” IEEE Trans. Biomed. Eng., 44(7), 1997.

Haykin, S. Nonlinear Methods of Spectral Analysis, 2nd ed., Berlin, Germany, Springer-Verlag, 1983.

Henning, G. and Husar, P. “Statistical Detection of Visually Evoked Potentials,” IEEE Engineering in Medicine and Biology, 14(4), 386-390, 1995.

Hinich, M. J., “Identification of the Coefficients in a Nonlinear Time Series of the Quadratic Type,” J. of Economics, 30: 269-288, 1985.

Husar, P. J., Leiner, B., et al., “Statistical Methods for Investigating Phase Relations in Stochastic Processes,” IEEE Trans. on Audio and Electroacoustics, 19(1): 78-86, 1971.

Husar, P. and Henning, G. “Bispectrum Analysis of Visually Evoked Potentials,” IEEE Engineering in Medicine and Biology, 16(1), 1997.

Jensen, J. A., “Deconvolution of Ultrasound Images,” Ultrasonic Imaging, 14:1-15, 1992.

Jensen, J. A. and Leeman, S. “Nonparametric Estimation of Ultrasound Pulses,” IEEE Trans. Biomed. Eng., 41(10): 929-936, 1994.

Kay, S. M. Modern Spectral Estimation, Prentice-Hall, Inc., Englewood Cliffs, NJ, 1988.

Kim, Y. C. and Powers, E. J. “Digital Bispectral Analysis of Self-Excited Fluctuation Spectra,” Phys. Fluids, 21(8): 1452-1453, 1978.

Le Roux, J, Coroyer, C., and Rossille, D. “Illustration of the Effects of Sampling on Higher-Order Spectra,” Signal Processing, 36: 375-390, 1994.

Le Roux, J. and Sole, P. “Least-Squared Error Reconstruction of a Deterministic Sampled Signal Fourier Transform Logarithm from its Nth Order Polyspectrum Logarithm,” Signal Processing, 35: 75-81, 1994.

Lii, K. S. and Rosenblatt, M. “Deconvolution and Estimation of Transfer Function Phase and Coefficients for Nongaussian Linear Processes,” The Annals of Statistics, 10: 1195-1208, 1982.

Marple, Jr., S. L. Digital Spectral Analysis with Applications, Prentice-Hall, Inc., Englewood Cliffs, NJ, 1987.

Matsuoka, T. and Ulrych, T. J. “Phase Estimation Using Bispectrum,” Proc. of IEEE, 72: 1403-1411, Oct., 1984.

Mendel, J. M. “Tutorial on Higher-Order Statistics (Spectra) in Signal Processing and System Theory: Theoretical Results and Some Applications,” IEEE Proc., 79: 278-305, 1991.

Nikias, C. L. and Petropulu, A P. Higher-Order Spectra Analysis: A Nonlinear Signal Processing Framework, Prentice Hall Inc., Englewood Cliffs, NJ, 1993.

Nikias, C. L. and Raghuveer, M. R. “Bispectrum Estimation: A Digital Signal Processing Frame­work,” Proceedings IEEE, 75(7): 869-891, 1987.

Ning, T. and Bronzino, J. D. “Bispectral Analysis of the EEG During Various Vigilance States,” IEEE Transactions on Biomedical Engineering, 36(4): 497-499, 1989.

Ning, T. and Bronzino, J. D. “Autogressive and Bispectral Analysis Techniques: EEG Applications,” IEEE Engineering in Medicine and Biology, March 1990.

Oppenheim, A. V. and Schafer, R. W. “Discrete-Time Signal Processing, Prentice-Hall, Englewood Cliffs, NJ, 1989.

Orizio, C., Liberati, D., Locatelli, C., De Grandis, D., Veicsteinas, A., “Surface Mechanomyogram Reflects Muscle Fibres Twitches Summation,” J. of Biomech. 29(4), 1996.

Pan, R. and Nikias, C. L. “The Complex Cepstrum of Higher-Order Cumulants and Nonminimum Phase System Identification,” IEEE Trans. Acoust, Speech, Signal Processing, 36: 186-205, 1988.

Papoulis, A. Probability Random Variables and Stochastic Processes, New York: McGraw-Hill, 1984.

Petropulu, A. P. and Nikias, C. L. “Signal Reconstruction from the Phase of the Bispectrum,” IEEE Trans. Acoust., Speech, Signal Processing, 40: 601-610, 1992.

Petropulu, A. P. “Noncausal Nonminimum Phase ARMA Modeling of Non-Gaussian Processes,” IEEE Trans. on Signal Processing, 43(8): 1946-1954, 1995.

Petropulu, A P. and Abeyratne, U. R. “System Reconstruction from Higher-Order Spectra Slices,” IEEE Trans. Sig. Proc. , 45(9): 2241-2251, 1997.

Petropulu, A. P. “Higher-Order Spectra in Signal Processing,” in Signal Processing Handbook, CRC Press, Boca Raton, FL, 1998.

Petropulu, A. P. and Pozidiz, H. “Phase Estimation from Bispectrum Slices,” IEEE Trans. Sig., Proc. 46(2): 527-531, 1998.

Powers, E. J., Ritz, C. K., et al. “Applications of Digital Polyspectral Analysis to Nonlinear Systems Modeling and Nonlinear Wave Phenomena,” Workshop on Higher-Order Spectral Analysis, 73-77, Vail, CO, 1989.

Pozidis, H. and Petropulu, A. P. “System Reconstruction from Selected Regions of the Discretized Higher-Order Spectrum,” IEEE Trans. Sig. Proc., 46(12): 3360-3377, 1998.

Subba Rao, T. and Gabr, M. M. “An Introduction to Bispectral Analysis and Bilinear Time Series Models,” Lecture Notes in Statistics, 24, New York: Springer-Verlag, 1984.

Raghuveer, M. R. and Nikias, C. L. “Bispectrum Estimation: A Parametric Approach,” IEEE Trans. on Acous., Speech, and Sig. Proc., ASSP 33(5): 1213-1230, 1985.

Raghuveer, M. R. and Nikias, C. L. “Bispectrum Estimation via AR Modeling,” Signal Processing, 10:35-45, 1986.

Rangoussi, M. and Giannakis, G. B. “FIR Modeling Using Log-Bispectra: Weighted Least-Squares Algorithms and Performance Analysis,” IEEE Trans. Circuits and Systems, 38(3): 281-296, 1991.

Rozzario, N. and Papoulis, A. “The Identification of Certain Nonlinear Systems by Only Observing the Output,” Workshop on Higher-Order Spectral Analysis, 73-77, Vail, CO, 1989.

Samar, V. J, Swartz, K. P., et al. “Quadratic Phase Coupling in Auditory Evoked Potentials from Healthy Old Subjects and Subjects with Alzheimer’s Dementia,” IEEE Signal Processing Workshop on Higher-Order Statistics, 361-365, Tahoe, CA, 1993.

Schetzen, M. The Volterra and Wiener Theories on Nonlinear Systems, updated edition, Krieger Publishing Company, Malabar, FL, 1989.

“Bispectral Analysis of Visual Field Interactions,” Proc. of International Conference of the IEEE Engineering in Medicine and Biology, vol. 3, Amsterdam, Netherlands, 1996.

Swami, A. and Mendel, J. M. “ARMA Parameter Estimation Using Only Output Cumulants,” IEEE Trans. Acoust., Speech and Sig. Proc., 38: 1257-1265, 1990.

Tang, Y. and Norcia, A. M. “Coherent Bispectral Analysis of the Steady-State VEP,” Proc. Of Inter­national Conference of the IEEE Engineering in Medicine and Biology, 1995.

Taxt, T. “Comparison of Cepstrum-Based Methods for Radial Blind Deconvolution of Ultrasound Images,” IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control, 44(3), 1997.

Tick, L. J., “The Estimation of Transfer Functions of Quadratic Systems,” Technometrics, 3(4): 562-567, 1961.

Treitel, S. and Lines L. R. “Linear Inverse Theory and Deconvolution,” Geophysics, 47(8): 1153-1159, 1982.

Tugnait, J. “Fitting Non-Causal AR Signal Plus Noise Models to Noisy Non-Gaussian Linear Processes,” IEEE Trans. Automat. Control, 32: 547-552, 1987.

Tugnait, J. “Identification of Linear Stochastic Systems via Second — and Fourth-Order Cumulant Matching,” IEEE Trans. Inform. Theory, 33: 393-407, 1987.

Upshaw, B. “SVD and Higher-Order Statistical Processing of Human Nerve Signals,” Proc. of Int. Conf. of the IEEE Engineering in Medicine and Biology, 1996.

Yana, K., Mizuta, H., and Kajiyama, R. “Surface Electromyogram Recruitment Analysis Using Higher-Order Spectrum,” Proc. of International Conference of the IEEE Engineering in Medicine and Biology, 1995.

Micheli-Tzanakou, E. “ Neural Networks in Biomedical Signal Processing.” The Biomedical Engineering Handbook: Second Edition.

Ed. Joseph D. Bronzino

Boca Raton: CRC Press LLC, 2000

Updated: 13.10.2013 — 07:42