Table of Contents Author Guidelines Submit a Manuscript
International Journal of Antennas and Propagation
Volume 2016, Article ID 5063450, 9 pages
http://dx.doi.org/10.1155/2016/5063450
Research Article

Azimuth/Elevation Directional Finding with Automatic Pair Matching

Department of Electrical Engineering, Prince Mohamed Bin Fahd University, Al Khobar, Saudi Arabia

Received 11 April 2016; Revised 21 September 2016; Accepted 15 November 2016

Academic Editor: Aboulnasr Hassanien

Copyright © 2016 Nizar Tayem. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract

We addressed the problem of two-dimensional (2D) direction-of-arrival (DOA) elevation and azimuth angles estimation for multiple uncorrelated signals using L-shaped antenna array configuration. The key points of the proposed method are the following: () it obtains azimuth and elevation angles through construction of three cross-correlation matrices from the collected data of the received signals; this implies that the noise reduces significantly in the reconstructed data matrices; () it derives a parallel factor analysis (PARAFAC) model and applies trilinear least squares method to avoid pair matching problem between 2D DOA azimuth and elevation angles for multiple sources; () it does not require spectral peak searching; and () it has better 2D DOA estimation compared with signal parameters via rotational invariance technique and fourth-order signal parameters via rotational invariance technique. Simulation results demonstrate the estimation accuracy and the effectiveness of the proposed method.

1. Introduction

The problem of estimating the two-dimensional direction-of-arrival (2D DOA) azimuth and elevation angles of multiple incident sources plays an important role in many practical applications in the fields of wireless communication and signal processing such as radar, sonar, and multiple-input-multiple-output (MIMO) systems. Several 2D DOA estimation methods have been proposed in the literature considering different geometries of the antenna arrays such as circular array [1], parallel uniform linear array [2], uniform rectangular array [3], and L-shaped array [4, 5]. In particular, the L-shaped array which is composed of two orthogonal uniform linear arrays with one placed on the -axis and the other on the -axis has received considerable attention due to its geometric configuration, ease of implementation, use of less number of antennas, and higher estimation accuracy compared with other geometrical arrays.

Conventional methods of DOA estimation rely on the decomposition of the observation space into a signal subspace and a noise subspace. Two widely used subspace techniques, MUSIC [6] and ESPRIT [7], are computationally complex as they require either eigenvalue decomposition (EVD) of the sample covariance matrix or the singular value decomposition (SVD) of the received data matrix. Improved techniques which are simpler and less complex have been reported in the literature [611] which do not rely on either EVD or SVD. However, some of these methods [6, 7] suffer from heavy losses of the array aperture and encounter an estimation failure problem. To avoid aperture loss problem many cumulant-based methods [1217] have been proposed. But these methods are computationally intensive and require parameter pairing.

The 2D DOA methods using L-shaped array [3, 4] fail to resolve the pair matching problem for multiple incident sources which result in incorrect 2D estimation of azimuth and elevation angles and hence severe performance degradation. Many schemes have been proposed to solve the problem of pair matching [11, 1821]. These pair matching techniques have high computational cost and complexity and they do not always provide accurate pairing results.

A parallel factor analysis (PARAFAC) [22] model is a method that transfers low rank matrix decomposition to three-way arrays (TWAs); it was first introduced in psychometrics and flow injection analysis. It has been widely used in subspace domain in array signal processing and wireless communication areas [23, 24]. Several methods have applied PARAFAC model to alleviate the problem of parameter pairing (or pair matching) for 2D DOA estimation for azimuth and elevation angles [2529]. A trilinear decomposition-based blind 2D DOA estimation algorithm employing PARAFAC data model for parallel shaped array has been proposed in [17] to achieve automatic pair matching. One drawback of this method is that it requires high number of snapshots and has high computational complexity.

In this paper, we propose a novel 2D DOA estimation method that employs an L-shaped uniform antenna array based on a new computationally efficient cross-correlation with automatic pair matching based on PARAFAC data model. The constructed data consists of three cross-correlation matrices which contain information about azimuth and elevation angles of multiple uncorrelated signals that are spatially independent; this implies that the noise significantly reduces in the constructed data matrices. Further, we derive a PARAFAC model and apply trilinear least squares method to avoid pair matching problem. Compared with existing methods, the proposed method works with less than a hundred snapshots and has very accurate estimation and has lower computational complexity.

The rest of the paper is structured as follows. The signal model and proposed method are presented and developed in Section 2. Analysis of computational complexity is presented in Section 3. Simulation results are presented in Section 4. Conclusions are drawn in Section 5.

2. Signal Model and Proposed Method

The geometry of the proposed L-shaped array is shown in Figure 1. The distance between adjacent antenna elements is where with being the wavelength of the incident waveform. The arrays are divided into four subarrays: , , , and . Each linear subarray consists of antenna elements. The element located at point is considered as a reference element. Hence, the total number of antenna elements in both -axis and -axis is . Consider narrowband noncoherent sources in far-field of the antenna array.

Figure 1: Configuration of array geometry for 2D DOA azimuth and elevation angles.

collected signal vectors received at , , , and subarrays are defined as follows:where represents the snapshot index and superscript represents the transpose operation. The received signal vectors in (1) can be represented as follows:where the array manifold vector for th source , , and , . is the complex envelope vector of incident sources and at snapshot , , and are the Gaussian white noise vectors of zero mean and variance . and are array manifold matrices.

The matrices in (3) and in (5) are diagonal matrices containing information about the azimuth angle and elevation angle which can be presented as follows:

In the proposed method, we first construct three cross-correlation matrices ( subarray, subarray), ( subarray, subarray), and ( subarray, subarray) as follows:where the superscript represents the conjugate and transpose operations. Note that and are spatially independent of each other. Therefore, and , where matrix has a dimension of with all entries zero. Additive noises are not correlated with incident signals. For uncorrelated sources , and are wide-sense stationary sequences. As a consequence, the correlation matrix of the signal sources is a diagonal matrix where its entries represent the power of signal sources. The matrices since both matrices are diagonal.

The cross-correlation matrices in (10) are concatenated to form a new data matrix as follows:Parallel factor (PARAFAC) model is applied on the data in (11) along with trilinear least squares method to jointly estimate the correct pair of azimuth and elevation angles for each signal source.

Using the definition of PARAFAC model, the outer product of the three vectors , and can be expressed in a third-order tensor form as with typical elements defined as . can be expressed as a sum of tensor product:where , , and are th columns of the following load matrices , , and . These matrices for a given PARAFAC model can be defined as follows:

The PARAFAC decomposition in (12) can also be represented in another matrix form of 3D tensor matrix and can be represented using three slice matrices , , as follows:where is the Khatri-Rao product based on column wise Kronecker products and is Kronecker product.

The PARAFAC decomposition in (14) is considered to be essentially unique to arbitrary permutation and scaling under the condition ([31, 32]):where , , denote the maximum number of arbitrary linearly independent columns of matrices , , and , respectively. The arbitrary permutation and scaling implies that there exists a triplet matrix related to as follows:where is a permutation matrix and are arbitrary diagonal matrices satisfying .

On the basis of PARAFAC theorem, three-way array (TWA) can be constructed using (11) as follows:

Let ; according to the definition of Khatri-Rao product the matrix in (17) can be transformed as follows:where represents the row vector data built from the diagonal elements of the diagonal matrix . The uniqueness of (18), (19), and (20) will be guaranteed if the following inequality holds:

For different DOAs and independent sources and have Vandermonde structure with minimum rank equal to the number of sources, and also is a nonsingular matrix whose rank equals the number of incident sources. This implies that the minimum rank of ; for multiple incident sources and, therefore, (22) will be always guaranteed.

One of the methods of solving PARAFAC model in (18), (19), and (20) is trilinear alternative least squares (TALS) approach [24, 3133]. TALS method can be applied to solve the matrices , and and then estimate the azimuth and elevation angles. There are three basic steps behind TALS method: (a) update one of the matrices , and each time using least squares (LS) method, (b) continue updating the remaining matrices based on the LS results from the previous estimation step, and (c) repeat previous steps (a) and (b) until convergence of the LS cost function. The detailed procedure of estimating , and using TALS is as follows.

Define the cost functions for finding the matrices , and aswhere stands for Frobenius norm. Given the estimation of matrices and , the matrix can be found from (23) as follows:where represents the pseudoinverse of matrix . can be also obtained by minimizing the cost function in (24) and keeping and fixed.

Similarly, an estimation of matrix can be obtained as

The process in (27), (28), and (29) will continue until matrices , , and converge. PARAFAC along with TALS method can be initialized to speed up the convergence by exploiting an ESPRIT method idea on the concatenated data formed in (11).

The estimated matrix in (27) will be sufficient for 2D DOA estimation of azimuth and elevation angles. The TALS method guarantees the convergence but it is slow. For fast implementation of alternative least squares method to solve PARAFAC model in (23), (24), and (25), the COMFAC algorithm is employed which speeds up the least square fitting by utilizing a compressed version of the three-way data into smaller matrix dimensions. COMFAC MATLAB function will be used to estimate the matrices , , and (as described in [24]) as follows:where is the input data, is the number of sources, represents the number of iterations, and are the output identification matrices.

Now the diagonal matrices in (3) and are estimated from the identification matrix as follows:where is the estimated th entry of the diagonal matrix . Similarly, th diagonal entry of can be obtained as follows:

The azimuth angle and elevation angle for th source can be estimated as follows:

The estimated azimuth and elevation angles in (33) are for th source. In case of multiple sources, the following pairs are obtained: . The estimated matrices , , and have the same column permutation matrix. This implies automatic pair matching since th column of the steering matrix corresponds to th column of matrix .

The procedure of the 2D DOA proposed method is summarized as follows.

Step 1. Construct , and subarrays according to (2), (3), (4), and (5) from multiple snapshots of the array data ,  for.

Step 2. Estimate the cross-correlation matrices , , and from multiple snapshots in (10) as follows:

Step 3. Concatenate the estimated cross-correlation matrices according to (34).

Step 4. Construct the three-way array (TWA) according to (17).

Step 5. Apply alternative least squares method to estimate , , and which minimize the cost function in (23), (24), and (25).

Step 6. Repeat Step 5 until convergence.

Step 7. Obtain the estimated diagonal matrices and from the identification matrix .

Step 8. Estimate the 2D DOA azimuth and elevation angles according to (33).

3. Analysis of the Computational Complexity of the Proposed Method

The computational complexity of the proposed method is compared with that of 2D DOA fourth-order cumulant method [17] and novel 2D DOA with L-shaped array [30]. For total snapshots, number of antenna elements, number of iterations, and number of sources, considering major processing operations like forming the sample covariance or cross-correlation matrices and applying the alternative least squares method, the total computational complexity of the proposed method is in the order of . The complexity of the novel L-shaped array method is in the order of and the complexity of the fourth-order cumulant method is in the order of . Upon comparison, it can be seen that the proposed method requires slightly less computational complexity compared to the method in [17] and significantly less computational complexity compared with the method in [30].

4. Simulation Results

The performance of the proposed method is presented in the section and compared with the novel L-shaped method in [30] and cumulant-based method in [17]. The performance is measured in terms of root mean square error (RMSE) for the azimuth and elevation angles estimation. We consider 21 antenna elements in total for the first three experiments. The distance between the adjacent elements is taken to be half the wave length of the incoming signal, and the number of uncorrelated sources is taken as (, and ). Several simulation experiments have been conducted to evaluate the performance of the proposed method. The RMSE for the joint 2D DOA estimation azimuth and elevation angles is defined as follows:where represents the source index, represents the expectation value of a random variable , and are the pair of the estimated elevation and azimuth angles.

In the first experiment, we consider two uncorrelated sources with direction-of-arrival azimuth and elevation angles , SNR range is set from −5 to 30 dB, and the number of snapshots is 200. Monte-Carlo trials of 500 are used. The RMSE values for source 1 and source 2 are shown in Figures 2 and 3 versus SNR for both novel L-shaped and cumulant methods and compared with the proposed method. We observe that the proposed method has better performance which is indicated through lower RMSE especially at low SNR. For a given RMSE value of 0.15 degrees for source 1 and 0.25 degrees for source 2, it is clear that the proposed method is 5 dB better compared with novel L-shaped for source 1 and around 2.5 dB for source 2. It is also 7.5 and 8 dB better when compared with cumulant method for source 1 and source 2, respectively. We observe also that proposed method has better performance compared with [17, 30] even at low SNR. It can be deduced from Figures 2 and 3 that the performance of the proposed method will be affected due to noise. For example, for the proposed method the RMSE at SNR of −5 dB for source 1 is about 1.3 degrees and for source 2 is about 1.5 degrees.

Figure 2: Comparison of RMSE of joint azimuth-elevation versus SNR for source 1 located at using proposed method, novel L-shaped, and cumulant method.
Figure 3: Comparison of RMSE of joint azimuth-elevation versus SNR for source 2 located at using proposed method, novel L-shaped, and cumulant method.

In the second experiment, four uncorrelated sources with DOAs at , , and ( are considered and the number of snapshots is set to 200. Monte-Carlo trials of 200 are conducted. SNR is set to 10 dB. Figures 4, 5, and 6 illustrate joint azimuth and elevation scatter diagrams of 2D DOA estimation for the proposed method, novel L-shaped, and cumulant method, respectively. It is shown that the four incoming sources can be clearly observed by all methods. The proposed method in Figure 4 gives better and precise estimation compared to the other two methods in [17, 30]. It is also observed that the cumulant method in Figure 6 has the worst azimuth and elevation estimation for the four sources.

Figure 4: Scatter plot for azimuth and elevation for four uncorrelated sources at , , , and ( at SNR = 10 dB using proposed method.
Figure 5: Scatter plot for azimuth and elevation for four uncorrelated sources at , , , and ( at SNR = 10 dB using novel L-shaped method [30].
Figure 6: Scatter plot for azimuth and elevation for four uncorrelated sources at , , , and ( at SNR = 10 dB using cumulant method [17].

In the third experiment, the effect of the number of snapshots on the performance of the proposed method is evaluated. We consider two uncorrelated sources with DOAs at , SNR is set to 10 dB, and the number of snapshots range is set from 100 to 600. Monte-Carlo trials of 500 are used. The average RMSE of azimuth and elevation angles estimation versus the number of snapshots for the two sources is shown in Figures 7 and 8. From these figures, we observe that RMSE of joint azimuth and elevation angles for source 1 and source 2 decrease with increasing number of snapshots. We can also clearly see that the proposed method has higher estimation accuracy compared to the methods in [17, 30].

Figure 7: RMSE of azimuth-elevation versus number of snapshots at SNR = 10 dB for source 1 located at using proposed method, novel L-shaped, and cumulant method.
Figure 8: RMSE of azimuth-elevation versus number of snapshots at SNR = 10 dB for source 2 located at using proposed method, novel L-shaped, and cumulant method.

In the fourth experiment, the effect of the number of antennas on the performance of the proposed method is evaluated. We consider a single source located at , SNR set to 6 dB, and the number of snapshots set to 300. Monte-Carlo trials of 400 are used. The average RMSE of azimuth and elevation angles estimation versus the number of antennas is shown in Figure 9. From the figure, we observe that RMSE of joint azimuth and elevation angles decrease with increasing number of antennas. We can also clearly see that the proposed method has better performance compared to the methods in [17, 30].

Figure 9: RMSE of azimuth-elevation versus number of antennas at SNR = 6 dB for single source located at using proposed method, novel L-shaped, and cumulant method.

5. Conclusions

We have proposed a new method for 2D DOA azimuth and elevation angles estimation using L-shaped array. The proposed method has lower complexity and better performance compared with existing methods since constructed data matrices from cross-correlation are almost free of noise. PARAFAC model is derived for automatic pair matching of azimuth and elevation angles for multiple incident sources. In addition, the proposed method does not require spectral peak searching.

Competing Interests

The author declares that there is no conflict of interests regarding the publication of this paper.

References

  1. Y. Wu and H. C. So, “Simple and accurate two-dimensional angle estimation for a single source with uniform circular array,” IEEE Antennas and Wireless Propagation Letters, vol. 7, pp. 78–80, 2008. View at Publisher · View at Google Scholar · View at Scopus
  2. Y. Wu, G. Liao, and H. C. So, “A fast algorithm for 2-D direction-of-arrival estimation,” Signal Processing, vol. 83, no. 8, pp. 1827–1831, 2003. View at Publisher · View at Google Scholar · View at Scopus
  3. P. Heidenreich, A. M. Zoubir, and M. Rubsamen, “Joint 2-D DOA estimation and phase calibration for uniform rectangular arrays,” IEEE Transactions on Signal Processing, vol. 60, no. 9, pp. 4683–4693, 2012. View at Publisher · View at Google Scholar · View at Scopus
  4. N. Tayem and H. M. Kwon, “L-shape 2-dimensional arrival angle estimation with propagator method,” IEEE Transactions on Antennas and Propagation, vol. 53, no. 5, pp. 1622–1630, 2005. View at Publisher · View at Google Scholar · View at Scopus
  5. Y. Hua, T. K. Sarkar, and D. D. Weiner, “An L-shaped array for estimating 2-D directions of wave arrival,” IEEE Transactions on Antennas and Propagation, vol. 39, no. 2, pp. 143–146, 1991. View at Publisher · View at Google Scholar · View at Scopus
  6. S. Marcos, A. Marsal, and M. Benidir, “The propagator method for source bearing estimation,” Signal Processing, vol. 42, no. 2, pp. 121–138, 1995. View at Publisher · View at Google Scholar · View at Scopus
  7. J. E. F. Del Rio and M. F. Câtedra-Pérez, “The matrix pencil method for two-dimensional direction of arrival estimation employing an l-shaped array,” IEEE Transactions on Antennas and Propagation, vol. 45, no. 11, pp. 1693–1694, 1997. View at Publisher · View at Google Scholar · View at Scopus
  8. F.-J. Chen, S. Kwong, and C.-W. Kok, “ESPRIT-like two-dimensional DOA estimation for coherent signals,” IEEE Transactions on Aerospace and Electronic Systems, vol. 46, no. 3, pp. 1477–1484, 2010. View at Publisher · View at Google Scholar · View at Scopus
  9. Y. Zhang, Z. Ye, X. Xu, and J. Cui, “Estimation of two-dimensional direction-of-arrival for uncorrelated and coherent signals with low complexity,” IET Radar, Sonar and Navigation, vol. 4, no. 4, pp. 507–519, 2010. View at Publisher · View at Google Scholar · View at Scopus
  10. N. Tayem, M. Omer, and A. Abul Hussain, “DOA estimation method using R matrix of the QR factorized data and its prototype implementation on NI-PXI platform,” in Proceedings of the 33rd Annual IEEE Military Communications Conference (MILCOM '14), pp. 333–337, Baltimore, Md, USA, October 2014. View at Publisher · View at Google Scholar · View at Scopus
  11. T.-H. Liu and J. M. Mendel, “Azimuth and elevation direction finding using arbitrary array geometries,” IEEE Transactions on Signal Processing, vol. 46, no. 7, pp. 2061–2065, 1998. View at Publisher · View at Google Scholar · View at Scopus
  12. B. Porat and B. Friedlander, “Direction finding algorithms based on high-order statistics,” IEEE Transactions on Signal Processing, vol. 39, no. 9, pp. 2016–2024, 1991. View at Publisher · View at Google Scholar · View at Scopus
  13. P. Chevalier, L. Albera, A. Ferréol, and P. Comon, “On the virtual array concept for higher order array processing,” IEEE Transactions on Signal Processing, vol. 53, no. 4, pp. 1254–1271, 2005. View at Publisher · View at Google Scholar · View at MathSciNet · View at Scopus
  14. M. C. Dogan and J. M. Mendel, “Applications of cumulants to array Processing-Part I: aperture extension and array calibration,” IEEE Transactions on Signal Processing, vol. 43, no. 5, pp. 1200–1216, 1995. View at Publisher · View at Google Scholar
  15. P. Chevalier and A. Férréol, “On the virtual array concept for the fourth-order direction finding problem,” IEEE Transactions on Signal Processing, vol. 47, no. 9, pp. 2592–2595, 1999. View at Publisher · View at Google Scholar · View at Scopus
  16. P. Chevalier, A. Ferréol, and L. Albera, “High-resolution direction finding from higher order statistics: the 2q-MUSIC algorithm,” IEEE Transactions on Signal Processing, vol. 54, no. 8, pp. 2986–2997, 2006. View at Publisher · View at Google Scholar · View at Scopus
  17. J. Liang, “Joint azimuth and elevation direction finding using cumulant,” IEEE Sensors Journal, vol. 9, no. 4, pp. 390–398, 2009. View at Publisher · View at Google Scholar
  18. A. J. van der Veen, P. B. Ober, and E. F. Deprettere, “Azimuth and elevation computation in high resolution DOA estimation,” IEEE Transactions on Signal Processing, vol. 40, no. 7, pp. 1828–1832, 1992. View at Publisher · View at Google Scholar · View at Scopus
  19. A. L. Swindlehurst and T. Kailath, “Azimuth/elevation direction finding using regular array geometries,” IEEE Transactions on Aerospace and Electronic Systems, vol. 29, no. 1, pp. 145–156, 1993. View at Publisher · View at Google Scholar · View at Scopus
  20. S. Kikuchi, H. Tsuji, and A. Sano, “Pair-matching method for estimating 2-D angle of arrival with a cross-correlation matrix,” IEEE Antennas and Wireless Propagation Letters, vol. 5, no. 1, pp. 35–40, 2006. View at Publisher · View at Google Scholar · View at Scopus
  21. R. Bro, N. D. Sidiropoulos, and G. B. Giannakis, “Optimal joint azimuth-elevation and signal-array response estimation using parallel factor analysis,” in Proceedings of the 32nd Asilomar Conference on Signals, Systems & Computers, pp. 1594–1598, November 1998. View at Scopus
  22. R. Bro, “PARAFAC: tutorials and applications,” Chemometrics and Intelligent Laboratory Systems, vol. 38, no. 2, pp. 149–171, 1997. View at Publisher · View at Google Scholar
  23. N. D. Sidiropoulos, R. Bro, and G. B. Giannakis, “Parallel factor analysis in sensor array processing,” IEEE Transactions on Signal Processing, vol. 48, no. 8, pp. 2377–2388, 2000. View at Publisher · View at Google Scholar · View at Scopus
  24. N. D. Sidiropoulos, G. B. Giannakis, and R. Bro, “Blind PARAFAC receivers for DS-CDMA systems,” IEEE Transactions on Signal Processing, vol. 48, no. 3, pp. 810–823, 2000. View at Publisher · View at Google Scholar · View at Scopus
  25. Z.-S. Qi, Y. Guo, and B.-H. Wang, “Blind direction-of-arrival estimation algorithm for conformal array antenna with respect to polarisation diversity,” IET Microwaves, Antennas and Propagation, vol. 5, no. 4, pp. 433–442, 2011. View at Publisher · View at Google Scholar · View at Scopus
  26. L. Zou, J. Lasenby, and Z. He, “Direction and polarisation estimation using polarised cylindrical conformal arrays,” IET Signal Processing, vol. 6, no. 5, pp. 395–403, 2012. View at Publisher · View at Google Scholar · View at MathSciNet · View at Scopus
  27. W. Si, L. Wan, L. Liu, and Z. Tian, “Fast estimation of frequency and 2-D doas for cylindrical conformal array antenna using state-space and propagator method,” Progress in Electromagnetics Research, vol. 137, pp. 51–71, 2013. View at Publisher · View at Google Scholar · View at Scopus
  28. R. A. Harshman, “Foundation of the PARAFAC procedure: model and conditions for an ‘explanatory’multi-mode factor analysis,” UCLA Working Papers in Phonetics, vol. 16, pp. 1–84, 1970. View at Google Scholar
  29. L. Wan, W. Si, L. Liu, Z. Tian, and N. Feng, “High accuracy 2D-DOA estimation for conformal array using PARAFAC,” International Journal of Antennas and Propagation, vol. 2014, Article ID 394707, 14 pages, 2014. View at Publisher · View at Google Scholar · View at Scopus
  30. Z. Xiaofei, L. Jianfeng, and X. Lingyun, “Novel two-dimensional DOA estimation with L-shaped array,” EURASIP Journal on Advances in Signal Processing, vol. 2011, article 50, 2011. View at Publisher · View at Google Scholar
  31. A. Stegeman and N. D. Sidiropoulos, “On Kruskal's uniqueness condition for the Candecomp/Parafac decomposition,” Linear Algebra and Its Applications, vol. 420, no. 2-3, pp. 540–552, 2007. View at Publisher · View at Google Scholar · View at MathSciNet · View at Scopus
  32. J. B. Kruskal, “Three-way arrays: rank and uniqueness of trilinear decompositions, with application to arithmetic complexity and statistics,” Linear Algebra and Its Applications, vol. 18, no. 2, pp. 95–138, 1977. View at Google Scholar · View at MathSciNet
  33. R. Bro, “PARAFAC. Tutorial and applications,” Chemometrics and Intelligent Laboratory Systems, vol. 38, no. 2, pp. 149–171, 1997. View at Publisher · View at Google Scholar · View at Scopus