SIViP DOI 10.1007/s11760-016-0888-3
ORIGINAL PAPER
Bispectrum estimation using a MISO autoregressive model A. Tanju Erdem1 · Ali Ö. Ercan1
Received: 6 March 2015 / Revised: 17 September 2015 / Accepted: 14 March 2016 © Springer-Verlag London 2016
Abstract Bispectra are third-order statistics that have been used extensively in analyzing nonlinear and non-Gaussian data. Bispectrum of a process can be computed as the Fourier transform of its bicumulant sequence. It is in general hard to obtain reliable bicumulant samples at high lags since they suffer from large estimation variance. This paper proposes a novel approach for estimating bispectrum from a small set of given low lag bicumulant samples. The proposed approach employs an underlying MISO system composed of stable and causal autoregressive components. We provide an algorithm to compute the parameters of such a system from the given bicumulant samples. Experimental results show that our approach is capable of representing non-polynomial spectra with a stable underlying system model, which results in better bispectrum estimation than the leading algorithm in the literature. Keywords Bispectrum estimation · Bicumulant sequence · MISO autoregressive system · System identification
1 Introduction Bicumulant sequence and its frequency counterpart bispectra are third-order statistics that have been used extensively in analyzing nonlinear and non-Gaussian data [1]. In particular, they have been applied to system identification [2], data mining [3], and ECG signal analysis [4]. As opposed
B
A. Tanju Erdem
[email protected] Ali Ö. Ercan
[email protected]
1
Electrical and Electronics Engineering Department, Ozyegin University, Istanbul, Turkey
to autocorrelation sequence and power spectra, these thirdorder statistics retain information about the non-linearities and the Fourier phase of the underlying process responsible for generating the data. For example, in some applications such as the ones listed above, it is necessary to find out whether harmonically related frequency components are phase coupled, which might have resulted from quadratic non-linearities of the underlying system [1]. In these cases, power spectrum is not sufficient since it does not contain any phase information. On the other hand, bispectrum preserves the phase information, enabling the detection of such coupling. This useful property, however, comes at a price: Although it is possible to model an arbitrary power spectrum via a single-input single-output (SISO) linear system driven by white noise (a well-known result due to the power spectrum factorization theorem), modeling of an arbitrary bispectrum using a SISO linear system is not always possible [5]. In fact, it has been shown that the set of bispectra that can be modeled by a SISO linear system has Lebesque measure zero [6]. One can find a SISO linear system to approximately model an arbitrary bispectrum in the least squares sense [7]; however, in general, one needs a non-SISO system to exactly model an arbitrary bispectrum. It was previously shown that a system with multiplicity (SWM) can be used to exactly model a polynomial bispectrum [8]. An SWM is a multiple-input single-output (MISO) linear system whose output is defined as the sum of the outputs of a finite number of SISO moving average (MA) systems, all driven by mutually independent non-Gaussian white noise processes. Being composed of MA systems, however, a SWM cannot be used to model non-polynomial bispectra. Bicumulant values at high lags suffer from large estimation variance [1,9]. Thus, there have been several attempts [1, 10,11] to model non-polynomial bispectra from bicumulant
123
SIViP 16
17
18
19
20
21
17
11
12
13
14
15
20
18
12
7
8
9
10
14
19
19
13
8
4
5
6
9
13
18
20
14
9
5
2
3
5
8
12
17 16
5 4 3 2
n
1 0
21
15
10
6
3
1
2
4
7
11
−1
20
14
9
5
2
3
5
8
12
17
−2
19
13
8
4
5
6
9
13
18
−3
18
12
7
8
9
10
14
19
−4
17
11
12
13
14
15
20
−5
16
17
18
19
20
21
−5
−4
−3
−2
−1
0 m
1
2
3
4
5
Fig. 1 Illustration of the sixfold symmetry of bicumulant samples for p = 5. The (m, n) points in the non-redundant region O p are enumerated by an index printed in boldface. The points outside of O p are marked with the index of the point in O p that has identical bicumulant value
samples at low lags. Neither of these methods obtain a bispectrum that can match all given bicumulant samples and provide an underlying process model at the same time. Both [1] and [10] employ a single autoregressive (AR) system. Hence, they cannot guarantee a match for the given bicumulant sequence due to the theoretical limitation of a linear system in modeling an arbitrary bicumulant sequence as mentioned above. Our previous work [11] obtains a bispectrum that can match the given bicumulant samples; however, it fails to provide an underlying stable system model. In this paper, we introduce a novel MISO AR system model that can be used to estimate a non-polynomial bispectrum form a given set of bicumulant samples. In Sect. 2, we present the problem definition. Then the proposed MISO AR model that will be used for bispectrum estimation is given in Sect. 3. The system identification and the bispectrum estimation methodologies are presented in Sect. 4. Examples are provided via extensive simulations in Sect. 5 to demonstrate the performance of the proposed approach. The paper is concluded in Sect. 6.
Fig. 2 A pth-order MISO AR system consisting of p SISO AR subsystems. Each sub-system i is parametrized by gain βi and an FIR impulse response h i (m), which is ci,m if 1 ≤ m ≤ i and h i (m) = 0 otherwise
r (m, n) = r (n, m) = r (m − n, −n).
(2)
As a result of (2), bicumulant samples in the so-called nonredundant region O p = {(m, n): 0 ≤ n ≤ m ≤ p} completely represent the whole bicumulant sequence. This non-redundant region is illustrated in Fig. 1, where the (m, n) points in O p , p = 0, 1, 2, . . . , 5, are enumerated by an index printed in boldface. For example, O2 has the bicumulant samples indexed 1, 2, . . . , 6, whereas O3 has the bicumulant samples indexed 1, 2, . . . , 10. The points outside of O p are marked with the index of the point in O p that has identical bicumulant value due to (2). The term sixfold symmetry corresponds to the symmetry of the six triangles, the boundaries of which are represented by the dotted lines in the figure. Thus, the problem that we are addressing in this paper is to find a system model that matches a given set of bicumulant samples over O p , and using the system model, estimates the following bicumulant samples in O p˜ , for any desired p˜ > p in order to estimate the bispectrum.
3 A MISO autoregressive model for bicumulants 2 Problem definition The bicumulant sequence r (m, n) of a third-order stationary signal x(m) is defined as [1] r (m, n) = E[x(k)x(k + m)x(k + n)]
(1)
In order to estimate bispectrum based only on bicumulant samples in O p , we propose to use the following MISO AR system, which is composed of sum of p SISO AR subsystems as illustrated in Fig. 2, where x(m) =
p
xi (m) + β0 0 (m),
(3)
i=1
where E[·] denotes the expectation operator. Due to the symmetry of the above definition in m and n, and the stationarity of the signal x(m), the bicumulant sequence turns out to have a sixfold symmetry [5] which can be compactly expressed by the following equations
123
and xi (m) =
i j=1
ci, j xi (m − j) + βi i (m).
(4)
SIViP
Above, i (m) are mutually independent zero-mean nonGaussian white noise process with E[i3 (m)] = 1. In the figure, the FIR impulse response h i (m) = ci,m if 1 ≤ m ≤ i, and it is zero otherwise. Using the definitions (1) and (3), the bicumulant sequence r (m, n) of x(m) is obtained as follows r (m, n) = E [x(k)x(k + m)x(k + n)] p =E xi (k) + β0 0 (k)
(5)
i=1
·
j=1 p
x (k + n) + β0 0 (k + n)
(6)
=1
=E
⎧ p p p ⎨ ⎩
xi (k)x j (k + m)x (k + n)
i=1 j=1 =1
+ β0 0 (k + n)
p p
+ β0 0 (k + m)
xi (k)x j (k + m)
i=1 j=1 p p
xi (k)x (k + n)
i=1 =1
+ β02 0 (k
+ m)0 (k + n)
p
xi (k)
i=1
+ β0 0 (k)
p p
x j (k + m)x (k + n)
j=1 =1
+ β02 0 (k)0 (k
+ n)
p
+ β02 0 (k)0 (k + m)
x j (k + m)
j=1 p
(a)
=
i
(10)
ci, j E [xi (k)xi (k + m − j)xi (k + n)]
j=1
+ E [βi i (k + m)xi (k)xi (k + n)] =
i
ci, j ri (m − j, n) + βi3 δ(m, n).
(11) (12)
j=1
⎡ ⎤ p ·⎣ x j (k + m) + β0 0 (k + m)⎦
ri (m, n) = E[xi (k)xi (k + m)xi (k + n)]
Above, step (a) uses (4) for x(k+m), and the last term is equal to βi3 δ(m, n) since m ≥ n ≥ 0 in the non-redundant region O p ; thus, i (k + m) is independent of xi (k) and xi (k + n) unless m = n = 0. The bispectum can be found as the discrete-time Fourier transform of r (m, n) given in Eqs. (9) and (12). This model has the following advantages compared to the other approaches in the literature. References [1] and [10] also try to match the bicumulant samples with an underlying AR system; however, their system has p + 1 free parameters, whereas the bicumulant has ( p + 1)( p + 2)/2 samples. Therefore, the bicumulant samples can only be approximately matched. On the other hand, [11] matches the bicumulant with enough free parameters; however, there is no underlying system in their model, which might result in unstable bicumulants when it is used for estimating the unobserved samples. The system model presented above has ( p + 1)( p + 2)/2 free parameters, consisting of the subsystem AR coefficients ci, j , j = 1, . . . , i; i = 1, . . . , p and gains βi , i = 0, . . . , p. Thus, it has the potential to exactly match the given bicumulant. In addition, if the subsystems are chosen to be stable, the estimation will also be stable. In the next section, we will propose a method to compute the parameters of our system to achieve these goals.
x (k + n)
=1
+ β03 0 (k)0 (k + m)0 (k + n) p xi (k)xi (k + m)xi (k + n) =E i=1
+ β03 0 (k)0 (k + m)0 (k + n) =
p
ri (m, n) + β03 δ(m, n),
(7)
(8) (9)
i=1
where δ(m, n) = 1 if m = n = 0, and it is 0 otherwise. Above, we used the fact that the noise process 0 (k) and the processes xi (k) are zero mean and mutually independent of each other. The bicumulant sequence ri (m, n) of each AR process xi (m) can be found as follows:
4 System identification and bicumulant estimation method Given the bicumulant samples in O p , our strategy is to first identify the parameters of a stable MISO AR system given in Sect. 3, whose bicumulant matches the given bicumulant parameters in O p . Then using these system parameters and (12), we estimate the unknown bicumulant samples in O p˜ , where p˜ > p. The bispectrum then can be found by taking the Fourier transform of the estimated bicumulant. Since we are interested in modeling practical systems, it is safe to assume that all the sub-systems in Fig. 2 are stable and causal. This means that the poles of each sub-system must lie inside the unit circle. Therefore, it makes sense to solve for the poles of the sub-systems and convert the resulting poles to the model coefficients ci, j . To find the poles and the gains
123
SIViP
of the sub-systems that match the given bicumulant samples, we formulate the following optimization problem minimize r − g(x, β) over x, β subject to 0 ≤ x ≤ xu
(13) (14) (15)
Above, vector r contains the given bicumulant samples r (m, n) ∈ O p that we would like to model, which are ordered in r according to the indexing shown in Fig. 1. x is a vector containing information on the pole locations, β is a vector containing the gains βi , g(x, β) is a mapping from given poles and gains to the bicumulant samples ordered in the same way as r, xu is an upper bound on x, and 0 is a vector of zeros. Each sub-system i has i poles; thus, a total of p( p + 1)/2 poles are needed for a pth-order MISO AR system. Since ci, j are real, the poles come in complex conjugate pairs. For an ith-order AR sub-system, we assume i/2 complex conjugate pairs if i is even, and (i − 1)/2 complex conjugate pairs and one real pole if i is odd. Note that a conjugate pair with zero imaginary parts results in two more real, albeit repeated, poles. This gives us the following definition for x T x = xμT xθT ,
(16)
where xμ and xθ contain the magnitudes and the phases of the poles, respectively. With the given definition of x and the stability and causality requirement, the optimization constraints become 0 ≤ xμ ≤ α1,
(17)
0 ≤ xθ ≤ π 1,
(18)
where 0 and 1 are a vector of zeros and ones, respectively. The reason for having 0 < α < 1 in (17) is as follows. Constraining the magnitudes to 1 might result in poles on the unit circle, which violates the stability criterion; thus, they should be constrained to a value less than one. One might want to keep the magnitudes further away from the unit circle to have more damped results in the estimation of the further bicumulant samples, in which case, α can be chosen even smaller. Thus, Eqs. (17) and (18) give the upper bound xu in (15): T xu = α1T π 1T .
(19)
The mapping g from x and β to the bicumulant samples is described in Algorithm 1. In the algorithm, the function “get_poles” picks the poles of sub-system i only, and the function “poles_2_coefficients” converts the poles to system coefficients by first finding the polynomial with its roots
123
Algorithm 1 Function g(x, β) Inputs: x (poles’ info for the sub-systems); β (all gains). Outputs: rˆ (Modeled bicumulant samples of the MISO AR system) T 1: rˆ ← β03 0 . . . 0 2: for i = 1, 2, . . . , p do 3: xi ← get_poles(x, i) 4: {ci, j }ij=1 ← poles_2_coefficients(xi ) 5: Ci ← compute_C({ci, j }ij=1 ) T 6: ri ← Ci−1 βi3 0 . . . 0 7: r˜ i ← extend(ri , p) 8: rˆ ← rˆ + r˜ i 9: end for
equal to the given poles, which is equal to y i − ci,1 y i−1 · · · − ci,i , and reading off the system coefficients {ci, j }ij=1 from the polynomial. It is possible to write the family of equations in (12) for m = 0, 1, . . . , i and n = 0, 1, . . . , m, in matrix form, Ci ri = bi , where T (20) bi = βi3 0 . . . 0 T ri = ri (0, 0) ri (1, 0) ri (1, 1) ri (2, 0) . . . ri (i, i) . (21) The function “compute_C” exactly performs the construction of this matrix Ci ∈ Rk×k , where k = (i + 1)(i + 2)/2. From this, the computation of ri simply follows as ri = Ci−1 bi . Also note that support of ri is Oi , however, that of rˆ is O p . Thus, the function “extend” extrapolates each ri to the support O p using (12), so that rˆ and r˜ i can be summed in step 1 of Algorithm 1. Both of these operations also utilize the symmetry in (2) (Fig. 1), whenever the indices (m − j, n) ∈ / Oi , to replace them with the corresponding indices in Oi . For example, for p = 3 and i = 2, the matrix C2 is given as follows: ⎡
1 ⎢−c2,1 ⎢ ⎢ 0 C2 = ⎢ ⎢−c2,2 ⎢ ⎣ 0 0
0 1 −c2,1 −c2,1 −c2,2 0
−c2,1 −c2,2 1 0 −c2,1 0
0 0 0 1 0 −c2,2
0 0 −c2,2 0 1 −c2,1
⎤ −c2,2 0 ⎥ ⎥ 0 ⎥ ⎥. 0 ⎥ ⎥ 0 ⎦ 1 (22)
The extrapolation of r2 to r˜ 2 is performed as follows: r˜ 2 [ j] = r2 [ j], for j = 1, 2, . . . , 6
(23)
r˜ 2 [7] = r2 (3, 0) = c2,1 r2 (2, 0) + c2,2 r2 (1, 0) = c2,1 r˜ 2 [4] + c2,2 r˜ 2 [2],
(24)
r˜ 2 [8] = r2 (3, 1) = c2,1r2 (2, 1) + c2,2 r2 (1, 1) = c2,1 r˜ 2 [5] + c2,2 r˜ 2 [3],
(25)
SIViP
r˜ 2 [9] = r2 (3, 2) = c2,1r2 (2, 2) + c2,2 r2 (2, 1) = c2,1 r˜ 2 [6] + c2,2 r˜ 2 [5],
(26)
Table 1 Original and estimated model parameters of a second-order MISO AR system at SNR = 30 dB β0
β1
β2
Original
1.780
1.396
−1.594
Estimated
1.797
1.401
−1.599
c1,1
c2,1
c2,2
r˜ 2 [10] = r2 (3, 3) = c2,1 r2 (3, 2) + c2,2 r2 (3, 1) = c2,1 r˜ 2 [9] + c2,2 r˜ 2 [8].
(27)
Once the system parameters are found, given the known bicumulant samples in O p , we use the exact same steps that are exemplified in (23)–(27) in order to estimate the unknown bicumulant samples in O p˜ , where p˜ > p. Next, we provide some examples to demonstrate the performance of our approach.
5 Examples In this section, we provide the extensive simulation results of two experiments. In the first experiment, the goal is to validate the system identification approach presented above. In the second experiment, the bispectrum estimation performance of the proposed method is investigated. In the first experiment, we generated a second-order MISO AR system by choosing the gains and stable poles of the subsystems randomly. The gains, the pole magnitudes and pole angles are chosen randomly in the intervals [−2, 2], [0, 0.99] and [0, π ], respectively. Then we computed the resulting bicumulant samples of the MISO AR system with additive white Gaussian noise corresponding to 30 dB SNR. Then we estimated the MISO AR system’s parameters ci, j and βi using our approach with the given noisy bicumulant samples and we compared them with the original parameters. Since more than one underlying MISO AR system may result in the same bicumulant, we solved the optimization problem in (13)–(15) multiple times and selected the one that approximates the original ci, j and βi values best. This is because this first experiment is to validate that an underlying stable MISO AR system may be found that approximates the given bicumulant closely using our approach, and existence of other systems that also approximate the bicumulant well does not weaken our approach. However, in the following experiment, when the bispectrum estimation performance of the approach is tested, only the results that minimize the objective in (13) are reported. The results of the experiment above are reported in Table 1. Here, μ1,1 refers to the magnitude of the real pole of the firstorder AR sub-system, μ2,1 and θ2,1 refer to the magnitude and positive phase of the complex conjugate pole pair of the second-order AR sub-system. As can be seen from the table, our approach matched the system coefficients well, with some minor differences due to the additive noise. Then we repeated the same experiment for a third-order MISO AR system. The results of this experiment are given in Table 2. Again, the estimated system parameters from the noisy bicumulant samples are close to the original parameters.
Original
0.485
−0.852
−0.189
Estimated
0.463
−0.848
−0.180
μ1,1
μ2,1
θ2,1
Original
0.485
0.435
2.935
Estimated
0.462
0.424
3.141
Next, we iterated the experiment above many times for second- and third-order systems and additive noise levels corresponding to 20, 30, and 40 dB SNR. For an individual experiment , the RMSE performance of estimation of βi and ci, j are computed as follows RMSE() βi RMSE() ci, j
= =
p+1 1 (βi () − βˆi ())2 , p+1
(28)
i=1
2 (ci, j () − cˆi, j ())2 , (29) p( p + 1) p
i
i=1 j=1
where βi () and ci, j () are the actual system parameters at experiment , and βˆi (), cˆi, j () are their estimates respectively. For each system order (second or third) and SNR level (20, 30, or 40 dB) pair, we repeated the experiment 60 times and the median RMSE performances are reported in Table 3. To compare these results, we also computed the sample standard deviations of randomly generated MISO AR systems, which are 1.16 for βi values, 0.62 for ci, j values when p = 2, and 0.65 when p = 3. The estimates’ RMSE performances are much better than these sample standard deviations. Since the method fits the system parameters to match the given noisy bicumulant sequence, the performance degrades as the SNR decreases. However, since our approach inherently assumes a stable system, it does not overfit to the noise. Next, we tested our approach using bicumulant data that are estimated from observed data, rather than bicumulant data that are obtained by direct computation. For this purpose, we first randomly generated stable secondand third-order MISO AR systems. Then we generated a zero-mean, non-Gaussian i.i.d. noise processes i (m) with E[i2 (m)] = E[i3 (m)] = 1 for 100,000 time steps. To obtain one sample of a process, we followed the method of summing 100 random variables, each of which are mixtures of uniform random variables [12]. Then the generated random
123
SIViP Table 2 Original and estimated model parameters of a third-order MISO AR system at SNR = 30 dB
Table 3 Simulation results for the first experiment
β0
β1
β2
β3
0.035
−1.528
1.505
−1.178
−0.054
−1.488
1.476
−1.197
c1,1
c2,1
Original
0.019
Estimated
0.016
Original Estimated
c2,2
c3,1
c3,2
c3,3
0.967
−0.653
0.995
−0.673
−0.530
0.284
0.415
−0.472
0.357
0.369
μ1,1
μ2,1
θ2,1
μ3,1
μ3,2
θ3,2
Original
0.403
Estimated
0.342
0.277
2.267
0.176
0.580
2.270
0.561
2.136
0.528
0.545
2.691
SNR
20 dB
30 dB
40 dB
AR system order
2
3
2
3
2
3
Median RMSE for βi
0.379
0.581
0.315
0.406
0.134
0.295
Median RMSE for ci, j
0.212
0.411
0.153
0.267
0.116
0.228
Table 4 Simulation results for the first experiment using estimated bicumulants from observed data AR system order
2
3
Median RMSE for βi
0.217
0.368
Median RMSE for ci, j
0.193
0.296
processes are passed through the MISO AR system to obtain the observed data x(k). Then to estimate the bicumulant, which is defined in (1), the following average is computed. T− p−1 1 x(k)x(k + m)x(k + n), rˆ (m, n) = T−p
(30)
k=0
where T = 100, 000 is the total number of time steps and p is the order of the MISO AR system ( p = 2 or 3). Since this method of obtaining the bicumulant already contains estimation error, no noise is added. Then following the approach given above, we estimated the system parameters ci, j and βi and computed the RMSE performance of the estimation over 60 runs. The median RMSE performances are given in Table 4. These RMSE values fall within the range of RMSE values in Table 3, validating the conclusions drawn using direct computation of bicumulants. To test the bicumulant estimation performance of our approach, we first construct a MISO system composed of three tenth-order MA sub-systems as shown in Fig. 3. This system also has a non-polynomial bispectrum; however, underlying system model is quite different than our proposed MISO AR model. This choice is made in order to demonstrate
123
Fig. 3 A MISO MA system consisting of three SISO MA sub-systems of order 10. Each sub-system i has a tenth-order FIR impulse response h i (m)
the performance of our approach with an arbitrary underlying system with non-polynomial bispectrum. We computed the bicumulant samples in O10 of the system in Fig. 3. Assuming that only samples in O3 are observed with additive white Gaussian noise corresponding to SNR = 30 dB, the goal in this experiment is to estimate the remaining samples of the bicumulant in O10 . Given the bicumulant estimate, bispectrum estimate can be found by taking its Fourier transform. We first identified a third-order MISO AR system to match the given noisy bicumulant samples in O3 using our approach. Then as an estimate of the remaining samples, we extrapolated the bicumulant samples of the resulting MISO AR system to O10 using (9) and (12). We also used the extrapolation approach in [11] using the same bicumulant samples in O3 and compared the resulting samples in O10 for both approaches with the originals. The results of this experiment are reported below. Due to space constraints, we show the samples in O6 only. Below, (31), (32), and (33) are the original (noiseless) bicumulant samples, those estimated using our approach, and the approach in [11], respectively. The
SIViP
representation in these equations is designed to reflect the triangular non-redundant region in Fig. 1. The vertical line is to mark the region used for identifying the model, namely O3 . r (m, n) ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ = −1.243 ⎪ ⎪ ⎪ ⎪ −0.429 −0.707 ⎪ ⎪ ⎪ ⎪ ⎪ 0.064 0.668 1.371 ⎪ ⎪ ⎪ ⎩ 1.789 −0.170 −2.684 0.683
SNR
20 dB
30 dB
40 dB
Method
Prop.
[11]
Prop.
[11]
Prop.
[11]
RMSEr
7.738
76.341
5.849
48.255
4.298
84.794
..
−0.673 . −1.039 0.924 . . . −0.027 −0.114 −0.084 . . . −0.271 0.353 0.013 . . . 0.338 0.081 0.372 . . . −0.350 −0.567 −0.309 . . . −1.284 −0.563 −1.491 . . .
(31) rˆ (m, n) ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ = −0.918 ⎪ ⎪ ⎪ ⎪ −0.722 −0.723 ⎪ ⎪ ⎪ ⎪ ⎪ 0.120 −0.092 1.138 ⎪ ⎪ ⎪ ⎩ 1.792 −0.133 −2.730 1.027
Table 5 Simulation results for the second experiment
. 0.029 . . −0.045 −0.010 . . . −0.132 −0.145 0.103 . . . −0.602 −0.164 0.022 . . . −0.269 −0.077 0.091 . . . 0.008 0.415 −0.268 . . . −0.839 1.012 −0.493 . . .
(32) rˆ[11] (m, n) ⎧ . ⎪ ⎪ ⎪ 21.139 . . ⎪ ⎪ ⎪ ⎪ ⎪ −12.323 10.967 . . . ⎪ ⎪ ⎪ ⎪ ⎨ 2.878 −0.290 −7.050 . . . = −1.288 −1.016 −0.441 0.223 . . . ⎪ ⎪ ⎪ ⎪ ⎪ −0.496 −0.640 2.513 10.282 −14.293 ... ⎪ ⎪ ⎪ ⎪ ⎪ 0.074 0.652 1.374 −6.630 −6.756 46.613 ... ⎪ ⎪ ⎩ 1.792 −0.131 −2.739 0.701 11.128 −16.568 −62.919 . . .
Table 6 Simulation results for the second experiment using estimated bicumulants from observed data
Method
Proposed
[11]
RMSEr
8.0405
94.5686
be as good as [11] in O3 , it does not have such an instability problem outside of O3 . Since practical systems that one would like to model are expected to exhibit stable behavior, the performance of the approach presented here should perform better in general. To show that the example presented above is not a seldom case, we performed the same experiment for SNR levels of 20, 30, and 40 dB, repeating each SNR case 60 times. For each experiment, we computed the RMSE performance of both our approach and the approach in [11] and reported the median RMSEs in Table 5. As seen here, the approach in [11] has the instability problem most of the time resulting in very high median estimation RMSE. Finally, similar to the previous experiment, we repeated the above experiment with the bicumulant data that are estimated using observed data. We did this experiment 60 times, and median RMSE performances of our approach and the approach in [11] are given in Table 6. These results support the results in Table 5, demonstrating that the proposed approach estimates the non-polynomial bispectra better due to underlying stable system model.
(33)
6 Conclusion Since the RMS errors in bicumulant and bispectrum estimation are equivalent due to Parseval’s theorem, we report the RMS errors for bicumulant estimation as a performance metric. For this example, the RMS error in the estimation of the bicumulant samples in O10 is 0.803 for our approach and 301.45 for the approach in [11]. As mentioned before, the approach in [11] can match the given noisy bicumulant samples in O3 exactly. However, since there is no underlying system model here, this approach may overfit to the noise and the estimated bicumulant in O10 may perform very poorly, especially when the estimated bicumulant sequence is not stable. This is exactly what is observed above. The estimated bicumulant samples for the approach in [11] in (33) start to grow outside of O3 , reaching as high as 1488.2 at the further samples of the sequence not shown above. However, as our approach inherently assumes an underlying stable system model, even though the matching performance might not
We introduced a MISO AR system model to estimate nonpolynomial bispectra from a finite set of given bicumulant samples. We presented an algorithm using numerical optimization to compute the parameters of such a stable system. The proposed model has the same number of parameters as the bicumulant samples in a triangular non-redundant region; thus, it allows matching non-polynomial bispectra in an inherently stable manner. We demonstrated via examples and extensive simulations the performance of our method using bicumulant data that are both computed directly with additive measurement noise with different SNR levels and bicumulant data that are estimated using observed data. We compared the performance of the proposed method to the leading algorithm in the literature. We observed that the leading algorithm can match the given bicumulant samples exactly; however, this might lead to overfitting and unsta-
123
SIViP
ble bispectrum estimation. On the other hand, the proposed approach is shown to not have this instability problem due to the underlying stable system model, as well as being capable of representing non-polynomial spectra, resulting in better estimation. Compliance with ethical standards Conflict of interest The authors declare that they have no conflict of interest.
References 1. Raghuveer, M.R., Nikias, C.L.: Bispectrum estimation: a parametric approach. IEEE Trans. Acoust. Speech Signal Process. 33(5), 1213–1230 (1985) 2. Chen, B., Petropulu, A.P.: Frequency domain blind MIMO system identification based on second and higher order statistics. IEEE Trans. Signal Process. 49(8), 1677–1688 (2001) 3. Mookiah, M.R.K., Acharya, U.R., Lim, C.M., Petznick, A., Suri, J.S.: Data mining technique for automated diagnosis of glaucoma using higher order spectra and wavelet energy features. Knowledge-Based Syst. 33, 73–82 (2012). doi:10.1016/j.knosys. 2012.02.010 4. Martis, R.J., Acharya, U.R., Mandana, K., Ray, A., Chakraborty, C.: Cardiac decision making using higher order spectra. Biomed. Signal Process. Control 8(2), 193–203 (2013). doi:10.1016/j.bspc. 2012.08.004
123
5. Tekalp, A.M., Erdem, A.T.: Higher-order spectrum factorization in one and two dimensions with applications in signal modeling and nonminimum phase system identification. IEEE Trans. Acoust. Speech Signal Process. 37(10), 1537–1549 (1989) 6. Erdem, A.T., Tekalp, A.M.: On the measure of the set of factorizable polynomial bispectra. IEEE Trans. Acoust. Speech Signal Process. 38(9), 1637–1639 (1990) 7. Roux, J.L., Sole, P., Tekalp, A.M., Erdem, A.T.: Tekalp-Erdem estimator gives the LS estimate for Fourier phase and log-Fourier modulus. IEEE Trans. Signal Process. 41(4), 1705–1707 (1993) 8. Erdem, A.T., Tekalp, A.M., Chang, M.M.: Modeling arbitrary polynomial bispectra in one and two dimensions. IEEE Trans. Signal Process. 40(4), 823–833 (1992) 9. Basumallick, N., Narasimhan, S.V.: Improved bispectrum estimation based on modified group delay. Signal Image Video Process. 6(2), 273–286 (2012). doi:10.1007/s11760-011-0217-9 10. Raghuveer, M.R., Dianat, S.A.: On the existence of autoregressive models for third-order cumulant matching. IEEE Trans. Acoust. Speech Signal Process. 37(12), 1931–1938 (1989) 11. Erdem, A.T., Tekalp, A.M.: Matching-extrapolation of bicumulants of one-D signals using two-D AR modeling. IEEE Int. Conf. Acoust. Speech. Signal Process. 5, 3453–3456 (1991) 12. Ermak, D.L., Nasstrom, J.S.: A method for generating skewed random numbers using two overlapping uniform distributions. Tech. Rep. Lawrence Livermore National Laboratory (1995)
