57.1 

57.2 

57.3 
HOS Computation from Real Data Indirect Method • Direct Method 

57.4 
NonParametric Methods • Parametric Methods 

57.5 

Drexeо University 
57.6 
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, geophysics, 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 distribution 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. HigherOrder Spectra (HOS) (of order greater than 2), which are defined in terms of higherorder cumulants of the data, do contain such information [35]. The third order spectrum is commonly referred to as bispectrum and the fourthorder as trispectrum. The power spectrum is also a member of the higherorder spectral class; it is the secondorder spectrum.
HOS consist of higherorder moment spectra, which are defined for deterministic signals, and cumulant 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 higherorder spectra of order greater than 2). Due to this property, HOS are high signaltonoise 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, secondorder statistics (autocorrelation) have been heavily used because they are the result of leastsquares optimization criteria. However, an accurate phase reconstruction 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. Higherorder correlations between input and output can detect and characterize certain nonlinearities [56], and for this purpose several higherorder spectra based methods have been developed.
Definitions and Properties of HOS
In this chapter only random onedimensional processes are considered. The definitions can be easily extended to the twodimensional case [34].
The join moments of order r of the random variables x1, …, xn are given by [42]:
N
Where k1 + …+ kn = r, and O() is their joint characteristic function. The corresponding joint cumulants are defined as:
□ □ 0 • 
Cumxk1, □, xn 
For a stationary discrete time random process x(k), (k denotes discrete time), the moments of order n are given by:
Where Ј{.} denotes expectation. The nth order cumulants are functions of the moments of order up to
N, i. e.,
1st order cumulants:
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, …, Tn1) = 0 for n > 2. In other words, all the information about a Gaussian process is contained in its first and secondorder cumulants. This property can be used to suppress Gaussian noise, or as a measure for nonGaussianity in time series.
[P2] If X(k) is symmetrically distributed, then cX3 (t1, t2) = 0. Thirdorder cumulants suppress not only Gaussian processes, but also all symmetrically distributed processes, such as uniform, Laplace, and BernoulliGaussian.
[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, …, Tn1) = c5n (t1, t2, …, Tn1) + cwn (t1, t2, …, Tn1). 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, …, Tn1) = c5n (t1, t2, …, Tn1), for n > 2. In other words, in higherorder 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, …, Tn1) = mXn (t1, …, Tn1), for n < 3.
Higherorder 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, ron1), exists, and is defined to be the (n1)dimensional Fourier
Transform of the nth order cumulant sequence. In general, CXn (ro1, ro2, ron1) 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 nonGaussian 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 corresponding to higherorder 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 thirdorder case.
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[]
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 thirdorder moments and cumulants are identical, i. e., cx3 (T!, t2) = mx3(T!, t2).
Compute the average cumulants as:
K
Cxx [, T2 □□ — □ m3 [T1, T2 □ (57.10)
K,□:
Obtain the thirdorder spectrum (bispectrum) estimate as the twodimensional 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 < M1, and w(t1, t2) is a twodimensional 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:
Where
L L 
Which is known as the minimum bispectrum bias supremum [36].
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 thirdorder 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 thirdorder spectrum is given as the average over all thirdorder 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 (5719)
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.
Let x(k) be generated by exiting a linear, timeinvariant (LTI), exponentially stable, mixedphase system with frequency response H(ro) (or impulse response h(k)), with an nth order white, zeromean, non — Gaussian stationary process v(k), i. e.,
Where * denotes convolution.
The output Hth order cumulant equals [10]:
C;D,□,□□□:□ hD®nqDhin33 (57.21)
K=0
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 nonGaussian 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= CX1 D □ 2, □ , □ n2. (57.23)
For example, the power spectrum of a nonGaussian linear process can be reconstructed from the bispectrum up to a scalar factor, i. e.,
While for n = 2, Eq. 57.22 provides the magnitude only of H(ro), as a function of CXn(m 1,ro2, …, ron1), 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].
Let
Y&Q x[ Q [ (57.25)
Where X(k) was defined in Eq. 57.20, and w(k) is zeromean Gaussian noise, uncorrelated to X(k). Methods for reconstructing h(k) from y(k) can be divided into two main categories: nonparametric and parametric.
Nonparametric 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, onedimensional 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 polyspectrum 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, …, N1}, {Cy3(m, l + r), m = 0,…, N1} 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 thirdorder spectrum of the system output, as long as the distance between slices, i. e., r and N are coprime 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:
Where
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 twodimensional phase unwrapping. The bispectrum needed in this method can be substituted 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 lowrank 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.
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 zeromean nth order white nonGaussian 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 YuleWalker equations can be derived based on thirdorder
Cumulants of y(k), i. e.,
P
□ ak^i, j[0, □□ q, (57.30)
I □ 0
Or
P
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 thirdorder 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 lowrank 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’)z1 can be constructed. Based on the filteredthrough 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].
FIGURE 57.1 Secondorder Volterra system. Linear and quadratic parts are connected in parallel.
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 □dQxK □□i [] (57.34)
Where h^Tp…, t) are the Volterra kernels of the system, which are symmetric functions of their arguments; for causal systems h^Tp… t) = 0 for any t; < 0.
The output of a secondorder Volterra system when the input is zeromean 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 crosscovariance of input and output, and
CX^r^Q □QOS □ □, Q& Q my] (57.37)
Be the thirdorder crosscumulant sequence of input and output. The system’s linear part can be obtained by [61]
f^n, (57.38)
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 nonGaussian 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 secondorder nonlinearity, is the quadraticphase coupling. There are situations where nonlinear interaction between two harmonic components of a process contribute 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]
6
Xk[ □ cosk k□, ;[ (57.42)
M
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
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
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 nonGaussian 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 radiofrequency (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 component 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 signaltonoise 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, nonionizing 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 BScan images, the resolution is compromised due to: (1) the finite bandwidth of the ultrasonic transducer, (2) the nonnegligible beam width, and (3) phase aberrations and velocity variations arising from acoustic inhomogeneity of tissues themselves. The observed ultrasound image can be considered as a distorted version of the true tissue information. Along the axial direction the distortion is dominated by the pulseecho wavelet of the imaging system, while along the lateral direction the distortion is mainly due to finitewidth lateral beam profile. Provided that these distortions are known in advance, or noninvasively 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 tissuedependent components. Therefore, distortions must be estimated based on the backscattered RF data that lead to the Bmode 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 Bmode 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 nonGaussian; (A2) f(n) is white, nonGaussian 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 nonminimum phase [26], their estimation was mostly carried out using secondorder 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 cepstrumbased 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 Bscan images. It was demonstrated that the HOSbased distortion estimation and subsequent image deconvolution significantly improved resolution. For the case of breast data, in was demonstrated [3] that deconvolution via SOSbased 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 Bscan 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 UltraMark9 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 HOSbased nonparametric 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 coordinates 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 onedimensional 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 onedimensional 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(Bc) show that the result of lateral followed by axial deconvolution of the image of Fig. 57.3(A), using respectively HOS — and SOSbased distortion kernel estimates. The deconvolution was performed via the constrained Wiener Filter technique [62]. According to Fig. 57.3, tHe deconvolution in the RFdomain resulted in a significant reduction in speckle size. The speckle size appears smaller in the HOSdeconvolved image than in SOSdeconvolved 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 HOSbased approach was 1.8 times that of the SOSbased approach. The lateral resolution gain for the HOSbased 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 MIP9553227, the Whitaker Foundation, the National Institute of Health under grant 2P01CA5282307A1, 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.
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. “HigherOrder Spectra Based Deconvolution of Ultrasound Images,” IEEE Trans. Ultrason. Ferroelec., andFreq. Con., vol. 42(6): 10641075, November 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): 479490, 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 higherorder statistics based methods,” IEEE Trans. Ultrason. Ferroelec, and Freq. Con., vol. 44(6): 14091416, 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): 673685, 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): 15761588, April 1993.
Barnett, T. P., Johnson, L. C. et al., “Bispectrum analysis of electroencephalogram signals during waking and sleeping,” Science, 172: 401402, 1971.
Bartlet, H., Lohmann, A. W., and Wirnitzer, B. “Phase and amplitude recovery from bispectra,” Applied Optics, 23(18): 31213129, 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 kthorder Spectra,” In: Spectral Analysis of Time Series, B Harris, Ed., John Wiley and Sons, New York, NY, 189232, 1967.
Brillinger, D. R., “The Identification of a Particular Nonlinear Time Series System,” Biometrika, 64(3): 509515, 1977.
Cohen, F. N. “Modeling of Ultrasound Speckle with Applications in Flaw Detection in Metals,” IEEE Trans. on Signal Processing, 40(3): 624632, 1992.
Dianat, S. A. and Raghuveer, M. R. “Fast Algorithms for Phase and Magnitude Reconstruction from Bispectra,” Optical Engineering, 29: 504512, 1990.
Fonollosa, J. A. R. and Vidal, J. “System Identification Using a Linear Combination of Cumulant Slices,” IEEE Trans. Sig. Proc., 41: 24052411, 1993.
Giannakis, G. B. “Cumulants: A Powerful Tool in Signal Processing,” Proc. IEEE, 75, 1987.
Giannakis, G. B. and Mendel, J. M., “CumulantBased Order Determination of NonGaussian ARMA Models,” IEEE Trans. on Acoustics, Speech, and Signal Processing, 38: 14111423, 1990.
Giannakis, G. B. and Swami, A. “On Estimating Noncausal Nonminimum Phase ARMA Models of NonGaussian Processes,” IEEE Transactions Acoustics, Speech, and Signal Processing, 38(3): 478495, 1990.
Gurcan, M. N., Yardimci, Y., Cetin, A. E., Ansari, R. “Detection of Microcalcifications in Mammograms Using HigherOrder Statistics,” IEEE Sig. Proc. Letters., 4(8), 1997.
Hadjileontiadis, L. and Panas, S. M., “Adaptive Reduction of Heart Sounds from Lung Sounds Using FourthOrder Statistics,” IEEE Trans. Biomed. Eng., 44(7), 1997.
Haykin, S. Nonlinear Methods of Spectral Analysis, 2nd ed., Berlin, Germany, SpringerVerlag, 1983.
Henning, G. and Husar, P. “Statistical Detection of Visually Evoked Potentials,” IEEE Engineering in Medicine and Biology, 14(4), 386390, 1995.
Hinich, M. J., “Identification of the Coefficients in a Nonlinear Time Series of the Quadratic Type,” J. of Economics, 30: 269288, 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): 7886, 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:115, 1992.
Jensen, J. A. and Leeman, S. “Nonparametric Estimation of Ultrasound Pulses,” IEEE Trans. Biomed. Eng., 41(10): 929936, 1994.
Kay, S. M. Modern Spectral Estimation, PrenticeHall, Inc., Englewood Cliffs, NJ, 1988.
Kim, Y. C. and Powers, E. J. “Digital Bispectral Analysis of SelfExcited Fluctuation Spectra,” Phys. Fluids, 21(8): 14521453, 1978.
Le Roux, J, Coroyer, C., and Rossille, D. “Illustration of the Effects of Sampling on HigherOrder Spectra,” Signal Processing, 36: 375390, 1994.
Le Roux, J. and Sole, P. “LeastSquared Error Reconstruction of a Deterministic Sampled Signal Fourier Transform Logarithm from its Nth Order Polyspectrum Logarithm,” Signal Processing, 35: 7581, 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: 11951208, 1982.
Marple, Jr., S. L. Digital Spectral Analysis with Applications, PrenticeHall, Inc., Englewood Cliffs, NJ, 1987.
Matsuoka, T. and Ulrych, T. J. “Phase Estimation Using Bispectrum,” Proc. of IEEE, 72: 14031411, Oct., 1984.
Mendel, J. M. “Tutorial on HigherOrder Statistics (Spectra) in Signal Processing and System Theory: Theoretical Results and Some Applications,” IEEE Proc., 79: 278305, 1991.
Nikias, C. L. and Petropulu, A P. HigherOrder 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 Framework,” Proceedings IEEE, 75(7): 869891, 1987.
Ning, T. and Bronzino, J. D. “Bispectral Analysis of the EEG During Various Vigilance States,” IEEE Transactions on Biomedical Engineering, 36(4): 497499, 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. “DiscreteTime Signal Processing, PrenticeHall, 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 HigherOrder Cumulants and Nonminimum Phase System Identification,” IEEE Trans. Acoust, Speech, Signal Processing, 36: 186205, 1988.
Papoulis, A. Probability Random Variables and Stochastic Processes, New York: McGrawHill, 1984.
Petropulu, A. P. and Nikias, C. L. “Signal Reconstruction from the Phase of the Bispectrum,” IEEE Trans. Acoust., Speech, Signal Processing, 40: 601610, 1992.
Petropulu, A. P. “Noncausal Nonminimum Phase ARMA Modeling of NonGaussian Processes,” IEEE Trans. on Signal Processing, 43(8): 19461954, 1995.
Petropulu, A P. and Abeyratne, U. R. “System Reconstruction from HigherOrder Spectra Slices,” IEEE Trans. Sig. Proc. , 45(9): 22412251, 1997.
Petropulu, A. P. “HigherOrder 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): 527531, 1998.
Powers, E. J., Ritz, C. K., et al. “Applications of Digital Polyspectral Analysis to Nonlinear Systems Modeling and Nonlinear Wave Phenomena,” Workshop on HigherOrder Spectral Analysis, 7377, Vail, CO, 1989.
Pozidis, H. and Petropulu, A. P. “System Reconstruction from Selected Regions of the Discretized HigherOrder Spectrum,” IEEE Trans. Sig. Proc., 46(12): 33603377, 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: SpringerVerlag, 1984.
Raghuveer, M. R. and Nikias, C. L. “Bispectrum Estimation: A Parametric Approach,” IEEE Trans. on Acous., Speech, and Sig. Proc., ASSP 33(5): 12131230, 1985.
Raghuveer, M. R. and Nikias, C. L. “Bispectrum Estimation via AR Modeling,” Signal Processing, 10:3545, 1986.
Rangoussi, M. and Giannakis, G. B. “FIR Modeling Using LogBispectra: Weighted LeastSquares Algorithms and Performance Analysis,” IEEE Trans. Circuits and Systems, 38(3): 281296, 1991.
Rozzario, N. and Papoulis, A. “The Identification of Certain Nonlinear Systems by Only Observing the Output,” Workshop on HigherOrder Spectral Analysis, 7377, 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 HigherOrder Statistics, 361365, 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: 12571265, 1990.
Tang, Y. and Norcia, A. M. “Coherent Bispectral Analysis of the SteadyState VEP,” Proc. Of International Conference of the IEEE Engineering in Medicine and Biology, 1995.
Taxt, T. “Comparison of CepstrumBased 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): 562567, 1961.
Treitel, S. and Lines L. R. “Linear Inverse Theory and Deconvolution,” Geophysics, 47(8): 11531159, 1982.
Tugnait, J. “Fitting NonCausal AR Signal Plus Noise Models to Noisy NonGaussian Linear Processes,” IEEE Trans. Automat. Control, 32: 547552, 1987.
Tugnait, J. “Identification of Linear Stochastic Systems via Second — and FourthOrder Cumulant Matching,” IEEE Trans. Inform. Theory, 33: 393407, 1987.
Upshaw, B. “SVD and HigherOrder 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 HigherOrder Spectrum,” Proc. of International Conference of the IEEE Engineering in Medicine and Biology, 1995.
MicheliTzanakou, E. “ Neural Networks in Biomedical Signal Processing.” The Biomedical Engineering Handbook: Second Edition.
Ed. Joseph D. Bronzino
Boca Raton: CRC Press LLC, 2000