CN104793177B - Microphone array direction-finding method based on least square method - Google Patents

Abstract

Based on the microphone array direction-finding method of least square method, comprise the following steps：The microphone array reception space far field sound-source signal being made up of N number of array element, obtains the sampled data of microphone array；Obtain the arrival delay inequality that sound-source signal reaches each array element of microphone array；Construction microphone array apart from difference vector；The location matrix of construction microphone array；Judge whether the array element of microphone array is generally aligned in the same plane, the two-dimensional directional vector of sound-source signal is if it is calculated using least square method, so as to obtain azimuth estimate and the angle of pitch estimate of sound-source signal；The unit direction vector estimate of sound-source signal is otherwise calculated using Constrained Least Square Methods, so as to obtain azimuth estimate and the angle of pitch estimate of sound-source signal.The angle of arrival of sound-source signal is solved and is converted into the corresponding direction vector of solution by the present invention, is solved the problems, such as the angle estimation of broadband sound source signal, is reduced the difficulty of angle estimation, and improve the accuracy of angle estimation.

Description

Microphone array direction-finding method based on least square method

Technical field

The invention belongs to signal processing technology field, more particularly to a kind of side for realizing voice signal based on microphone array Parallactic angle and the method for pitch angle measurement, for solving the real-time direction finding problem of wideband speech signal, can be applicable to communication or radar The fields such as radiation source direction finding are realized Deng based on array antenna.

Background technology

The sound signal collecting device that microphone array is made up of the one group of microphone that puts according to ad-hoc location, according to Practical application, microphone array signals process mainly solve be extract from the output result of array useful signal or for Carry out parameter Estimation.Microphone array direction finding has obtained extensive concern as the important research content of Speech processing.

The general principle for auditory localization being carried out based on reaching time-difference is the range difference that make use of sound source to reach different array elements Relevant with respect to the deflection of array with sound source, therefore also commonly referred to as range difference direction-finding method.As microphone array is general It is to be made up of multiple array elements, therefore how is become with the angle estimation performance that improves sound source using the metrical information of all array elements and work as The emphasis of front microphone array research.In prior art, the phase difference between signal is received using interferometer measurement difference array element, Again the incident angle for obtaining corresponding signal is resolved by ambiguity solution and phase difference.As voice signal is broadband signal, therefore very Hardly possible directly carries out the angle estimation of voice signal using traditional arrowband phase-interferometer method, if utilization space Power estimation side Method, it is necessary to transform to frequency domain, and carry out sub-band processing, not only realizes complicated, operand greatly, and cannot ensure that full frequency band is processed Uniformity.

Content of the invention

It is an object of the invention to provide a kind of direction-finding method of the microphone array for effectively improving real-time and accuracy.

To achieve these goals, the present invention takes following technical solution：

Based on the microphone array direction-finding method of least square method, comprise the following steps：The microphone being made up of N number of array element Array received space far-field sound-source signal, N >=3,

Step one, the sampled data of acquisition microphone array；

Step 2, acquisition sound-source signal reach the arrival delay inequality of each array element of microphone array；

Step 3, construction microphone array apart from difference vector r；

Microphone array apart from difference vector r=[Δ r₂,…,Δr_i,…,Δr_N]^T, wherein, Δ r_i=c τ_iBelieve for sound source Number reach i-th array element of microphone array range difference, τ_iThe arrival delay inequality of i-th array element is reached for sound-source signal, and c is sound Sound spread speed, i=2 ..., N；

Step 4, the location matrix A of construction microphone array；

Wherein, x_n、y_n、z_nFor the position coordinates of n-th array element, n=1 ..., N；

The unit direction vector of sound-source signalThe unit of sound-source signal Direction vector and the relation between difference vector are：A Θ=r, wherein, θ is the azimuth of sound-source signal,For sound-source signal The angle of pitch, []^TRepresent transposition operation；

Step 5, judge whether the array element of microphone array is generally aligned in the same plane, if it is execution step six, otherwise hold Row step 7；

Step 6, using least square method calculate sound-source signal two-dimensional directional vector, so as to obtain the side of sound-source signal Parallactic angle estimate and angle of pitch estimate；

When the array element of microphone array is generally aligned in the same plane, the unit direction vector Θ of location matrix A and sound-source signal divides Two-dimensional position matrix A 2 and two-dimensional directional vector Θ are not deteriorated to₂, according toCalculate the orientation of sound-source signal Angular estimation valueWith angle of pitch estimate

μ in formula₁For two-dimensional directional vector Θ₂Section 1 element, μ₂For two-dimensional directional vector Θ₂Section 2 element；

Step 7, using Constraint least square algorithm calculate sound-source signal unit direction vector estimateSo as to obtain The azimuth estimate and angle of pitch estimate of sound-source signal；

According toCalculate the azimuth estimate of sound-source signalWith angle of pitch estimate

μ in formula₁For unit direction vector estimateSection 1 element, μ₂For unit direction vector estimate's Section 2 element, μ₃For unit direction vector estimateSection 3 element, I be unit matrix, λ for estimate parameter.

The inventive method further scheme is, by checking that whether each row of location matrix A be in the step 5 Zero judging whether array element is generally aligned in the same plane, and is classified as zero if any one, then the array element of explanation microphone array is distributed in same flat On face, the array element of otherwise explanation microphone array is not distributed across on same plane.

The inventive method further scheme is, in the step 5 by check location matrix A whether sequency spectrum sentencing Whether disconnected array element is generally aligned in the same plane, and such as location matrix A is unsatisfactory for sequency spectrum, then the array element of explanation microphone array is distributed in same In one plane, the array element of otherwise explanation microphone array is not distributed across on same plane.

The inventive method further scheme is that the estimation parameter lambda passes through equation r^TA(A^TA+λI)^-2A^TR=1 is asked Solution, calculation procedure are as follows：

To matrix A^TA carries out feature decomposition, i.e.,Wherein, Λ=diag (γ₁, γ₂,γ₃) it is characterized value matrix, U=[u₁,u₂,u₃] corresponding feature matrix for it, N_dimFor the columns of location matrix A, N_dim=3, u_kIt is characterized value γ_kCorresponding characteristic vector, k=1,2,3, ()^HRepresent transposed complex conjugate operation；

According toAndObtain a_k= p_kq_k；

Estimate that parameter lambda passes through polynomial equation c₆λ⁶+c₅λ⁵+c₄λ⁴+c₃λ³+c₂λ²+c₁λ+c₀=0 is solved, multinomial Coefficient formulas in equation are as follows：

c₆=1,

c₅=2 γ₁+2γ₂+2γ₃,

Using polynomial rooting algorithm, equation root is tried to achieve according to above coefficient value, using root minimum for wherein modulus value as Required estimation parameter lambda.

The angle of arrival of sound-source signal is solved and is converted into the corresponding direction vector of solution by the inventive method, and then has obtained sound Source direction vector and sound source reach the linear equation between each array element range difference, and the angle for efficiently solving broadband sound source signal is estimated Meter problem；By required direction angular dimensions is converted into direction vector, it is achieved that be linear equation by non-linear equation, drop The low difficulty of angle estimation, while improve the accuracy of angle estimation, is redundancy observation bar when number of microphone is more Angle estimation under part provides effective processing method；Least square is asked by using the constant modulus property of Sounnd source direction vector Solution method enters row constraint, improves the estimated accuracy of azimuth and the angle of pitch.

Description of the drawings

In order to be illustrated more clearly that the embodiment of the present invention or technical scheme of the prior art, below will be to embodiment or existing Having in technology description needs the accompanying drawing for using to do simple introduction, it should be apparent that, drawings in the following description are only the present invention Some embodiments, for those of ordinary skill in the art, on the premise of not paying creative work, can also basis These accompanying drawings obtain other accompanying drawings.

Fig. 1 is the flow chart of the inventive method；

Fig. 2 is to calculate the distribution map for estimating parameter in emulation experiment using distinct methods；

Fig. 3 is the angle estimation error comparison diagram under different root of a polynomial；

Fig. 4 a to Fig. 4 d is respectively maximum distance error to angle-of- arrival estimation result influence curve comparison diagram；

Fig. 5 a to Fig. 5 d is respectively array element annular radii to angle estimated result influence curve comparison diagram；

Fig. 6 a to Fig. 6 d is respectively the otherness of array Z axis coordinate to angle estimated result influence curve comparison diagram；

Fig. 7 a to Fig. 7 d is respectively the evaluated error distribution map using Least Square Method angle of arrival；

Fig. 8 a to Fig. 8 d is respectively the evaluated error distribution map for estimating angle of arrival using Constraint least square algorithm.

Specific embodiment

In order to above and other objects of the present invention, feature and advantage can be become apparent from, the embodiment of the present invention cited below particularly, And coordinate appended diagram, it is described below in detail.

With reference to the flow chart that Fig. 1, Fig. 1 are the inventive method.The step of the inventive method, is as follows：It is made up of N number of array element Microphone array reception space far field sound-source signal, N are array number, N >=3,

Step one, the sampled data of acquisition microphone array；

Step 2, acquisition sound-source signal reach the arrival delay inequality of each array element of microphone array；

It is each that each channel data in the microphone array sampled data obtained using step one estimates that sound-source signal is reached Delay inequality between array element, according to microphone array structure choice reference array element, is obtained using Time Delay Estimation Method such as broad sense cross-correlation The arrival delay inequality that sound-source signal is reached between other array element and reference array element is taken, reference array element, ginseng are chosen according to array structure Array element is examined as the 1st array element, sound-source signal reaches arrival delay inequality τ of reference array element (the 1st array element)₁=0；

Step 3, construction microphone array apart from difference vector r；

Microphone array apart from difference vector r=[Δ r₂,…,Δr_i,…,Δr_N]^T, wherein, Δ r_i=c τ_iBelieve for sound source Number reach microphone array in i-th array element range difference, τ_iThe arrival delay inequality of i-th array element is reached for sound-source signal, and c is Sound propagation velocity, i=2 ..., N；

Step 4, the location matrix A of construction microphone array；

Wherein, x_n、y_n、z_nFor the position coordinates of n-th array element, n=1 ..., N, the position coordinates of array element is according to microphone The arrange parameter of array determines；

The unit direction vector of sound-source signal For unit The Section 1 element of direction vector Θ with regard to azimuth and the functional relation of the angle of pitch,Swear for unit direction Amount Θ Section 2 element with regard to azimuth and the functional relation of the angle of pitch,Section 3 for unit direction vector Θ Element with regard to azimuth and the functional relation of the angle of pitch, the unit direction vector of sound-source signal and between difference vector Relation is：A Θ=r；Wherein, θ is the azimuth of sound-source signal,For the angle of pitch of sound-source signal, []^TRepresent transposition operation, Azimuth be projection of the incident direction of arrival on X/Y plane and X-axis forward direction between angle, its span for (- π, π], bow The elevation angle is the angle of incident direction of arrival and X/Y plane, and its span is [0, pi/2]；

Step 5, judge whether the array element of microphone array is generally aligned in the same plane, if it is execution step six, otherwise hold Row step 7；

Judge whether the array element of microphone array is generally aligned in the same plane, by checking each row of location matrix A can be No is zero to be judged, is classified as zero if any one, then explanation array element is distributed the battle array that microphone array is otherwise described at grade Unit is not distributed across on same plane；Or by checking location matrix A whether sequency spectrum, if order r (A) of location matrix A= 2, that is, be unsatisfactory for sequency spectrum, then the array element of explanation microphone array is distributed at grade, otherwise, i.e. during r (A)=3, position Matrix A sequency spectrum, then illustrate that the array element of microphone array is not distributed across on same plane；Selected according to different array distribution Different angle estimating methods；

Step 6, using least square method calculate sound-source signal two-dimensional directional vector, so as to obtain the side of sound-source signal Parallactic angle estimate and angle of pitch estimate；

When the array element of microphone array is generally aligned in the same plane, then there is one-dimensional coordinate identical, location matrix A respective column Element be 0, now the location matrix A and unit direction vector Θ of sound-source signal deteriorates to two-dimensional position matrix A respectively₂With Two-dimensional directional vector Θ₂, the present embodiment illustrated so that Z axis coordinate is identical as an example, then

Accordingly, when X-axis coordinate is identical, Work as Y When axial coordinate is identical,

According toCalculate the azimuth estimate of sound-source signalWith angle of pitch estimate

μ in formula₁For two-dimensional directional vector Θ₂Section 1 element, μ₂For two-dimensional directional vector Θ₂Section 2 element, R is for microphone array apart from difference vector, A₂For the two-dimensional position matrix of microphone array, μ in the present embodiment₁And μ₂Counterparty The expression formula of parallactic angle and the angle of pitch is respectively：

Step 7, using Constraint least square algorithm calculate sound-source signal unit direction vector estimateSo as to obtain The azimuth estimate and angle of pitch estimate of sound-source signal；

When each array element is not in same plane, location matrix A now is sequency spectrum, i.e. r (A)=3, in order to enter one Step improves the direction finding precision of sound-source signal, enters row constraint using the constant modulus property of sound-source signal direction vector, using formulaSolution is optimized to direction vector, obtains unit direction vector estimate

According toCalculate the azimuth estimate of sound-source signalWith angle of pitch estimate

Or

I in formula is unit matrix, μ₁For unit direction vector estimateSection 1 element, μ₂Swear for unit direction Amount estimateSection 2 element, μ₃For unit direction vector estimateSection 3 element, μ₁, μ₂And μ₃Corresponding azimuth It is respectively with the expression formula of the angle of pitch： R for microphone array away from Deviation vector, A are the location matrix of microphone array, and λ is for estimating parameter.

Angle of pitch estimateAvailable formulaOr formulaCalculate, also can adopt simultaneously With both of the aforesaid formula calculate the angle of pitch, then average as last angle of pitch estimate, with obtain more accurate, can The estimated result for leaning on.

The estimation parameter lambda of the present invention can pass through equation r^TA(A^TA+λI)^-2A^TR=1 is solved, and calculation procedure is as follows：

To matrix A^TA carries out feature decomposition, i.e.,Wherein, Λ=diag (γ₁, γ₂,γ₃) it is characterized value matrix, U=[u₁,u₂,u₃] corresponding feature matrix for it, N_dimFor the columns of location matrix A, As A is sequency spectrum, therefore N_dim=3, u_kIt is characterized value γ_kCorresponding characteristic vector, k=1,2,3, ()^HRepresent transposition Complex conjugate operation；

According toAndObtain a_k= p_kq_k, k=1,2,3, p_kFor the kth item element of vector p, q_kKth item element for vector q；

The estimation parameter lambda of the present invention passes through polynomial equation c₆λ⁶+c₅λ⁵+c₄λ⁴+c₃λ³+c₂λ²+c₁λ+c₀=0 is asked Solution, in polynomial equation, each coefficient formulas are as follows：

c₆=1,

c₅=2 γ₁+2γ₂+2γ₃,

In above formula, γ_kIt is characterized k-th characteristic value of value matrix Λ, k=1,2,3, using polynomial rooting algorithm, root Go up coefficient value according to this and equation root is tried to achieve, as polynomial equation has six root λ₁、λ₂、λ₃、λ₄、λ₅、λ₆, will wherein modulus value be most Little root is used as required estimation parameter lambda.

The effect of the present invention can be further illustrated by following simulation result：

Experiment one,

Simulated conditions are as follows：It is uniform rings battle array that microphone array is classified as array element in X/Y plane, the Z axis coordinate tool of adjacent array element There is certain difference, the error in emulation is the wherein error by addition in the distance between target and array element aberration measurements Maximum be R/10, array number N=8, the X of all array elements, Y-coordinate be distributed on the annulus of R=0.2m, and the coordinate of Z axis is not With wherein odd indexed array element Z=0, and even number sequence number array element Z=R/5.The azimuth of sound-source signal and the angle of pitch for (70 °, 30°).

Fig. 2 is shown and is calculated the distribution map for estimating parameter lambda using distinct methods.Wherein be labeled as Fun-Opt+be employing The constraints estimation parameter lambda that obtains of search, be labeled as the zero of Err-Opt be using estimating that global search minimum error method is obtained Meter parameter lambda, be labeled as Roots × be using the inventive method polynomial rooting method obtain estimation parameter lambda.

Figure it is seen that when calculating optimum λ based on polynomial rooting method, the modulus value of R6 is minimum, and distance (0,0) is most Closely, in conjunction with Fig. 3, it should choose the minimum root of modulus value as required optimum λ value.But based on global search orientation and the angle of pitch The optimum λ that evaluated error is minimum and searches is plural number, and has one with the λ obtained using constraints and polynomial rooting method Fixed difference, this are caused due to both optimization aim differences.And the λ based on constraints acquisition and polynomial rooting method Otherness between the λ of acquisition is caused by the step-length in search procedure or search precision, is wherein based on polynomial rooting The optimum λ that method is obtained is most accurate, and speed is faster.

Fig. 3 is the angle estimation error comparison diagram under different root of a polynomial, the curve weight of CLSMerr and CLSMfun in Fig. 3 Close, from figure 3, it can be seen that based on different root of a polynomial as λ angle estimation error, the 6th root, that is, modulus value is minimum The corresponding error of root minimum.

Experiment two,

The arrangement of the microphone array array element of experiment two is identical with experiment one.Experiment two is minimum with constraint by least square method Square law is contrasted.Due to when the Z axis coordinate of array element is inconsistent, being estimated to be two formula and can utilize of the angle of pitch, because This, the method calculated using inverse cosine function is labeled as-F1, and the method calculated using arcsin function is labeled as-F2, will Both are labeled as-MF12 at the average method of result of calculation, are labeled as LSM using the curve of least square method, minimum using constraint The curve of square law is labeled as CLSM.

Fig. 4 a is maximum distance error azimuthal estimated result influence curve comparison diagram, and Fig. 4 b to Fig. 4 d is respectively maximum Range error affects contrast curve chart to angle of pitch estimated result.Modulus value be can be seen that from Fig. 4 a to Fig. 4 d and constrain azimuthal Have greatly improved with angle of pitch estimated result tool, wherein improved most substantially, from the emulation based on arcsin function result of calculation As a result it can also be seen that, three in Constraint least square algorithm kind angle of pitch estimated result is identical, and using anti-remaining in least square method The pitching angular accuracy highest that string function is calculated.

Fig. 5 a to Fig. 5 d is respectively array element annular radii to angle estimated result influence curve comparison diagram.From Fig. 5 a to Fig. 5 d As can be seen that with the increase of annular radii, the azimuth and angle of pitch evaluated error of least square and constraint least square by Gradual change is little, and the estimated result of Constraint least square algorithm azimuthal and the angle of pitch has certain improvement, is superior to least square method Azimuthal and the estimated result of the angle of pitch, the result that is wherein estimated using arcsin function method improve most obvious.

Fig. 6 a to Fig. 6 d is respectively the otherness of array Z axis coordinate to angle estimated result influence curve comparison diagram, from figure 6a to Fig. 6 d can be seen that to be increased with the otherness of Z axis coordinate, and evaluated error is tapered into, and modulus value constraint azimuthal Estimated result has certain improvement, and pitching angular estimation is improved than larger, especially when array Z axis otherness is smaller, utilizes Arcsine method estimated result improves most obvious.

Fig. 7 a to Fig. 7 d is respectively the evaluated error distribution map using Least Square Method angle of arrival.Fig. 8 a to Fig. 8 d divides It is not the evaluated error distribution map for estimating angle of arrival using Constraint least square algorithm.By contrast it is found that modulus value is constrained not Azimuthal estimated accuracy can be only improved, the estimated accuracy for the angle of pitch is improved most substantially, compared to least-squares algorithm, Azimuth and angle of pitch estimated accuracy all have a distinct increment, and especially the error of two kinds of angle of pitch method of estimation acquisitions is distributed almost Equally.

By experimental verification, not only direction finding speed is fast for the inventive method, and high precision, especially when utilization constraint is weighted most During little least square method, with higher azimuth and pitching angular estimation performance, positioning and the tracking field of voice signal is can be applicable to Scape.

The above, is only presently preferred embodiments of the present invention, not does any pro forma restriction to the present invention, though So the present invention is disclosed above with preferred embodiment, but is not limited to the present invention, any is familiar with this professional technology people Member, in the range of without departing from technical solution of the present invention, when the technology contents using the disclosure above make a little change or modification For the Equivalent embodiments of equivalent variations, as long as being the content without departing from technical solution of the present invention, the technical spirit of the foundation present invention Any simple modification, equivalent variations and the modification made by above example, all still falls within the range of technical solution of the present invention.

Claims (4)

1. the microphone array direction-finding method based on least square method, it is characterised in that comprise the following steps：By N number of array element group The microphone array reception space far field sound-source signal for becoming, N >=3,

Step one, the sampled data of acquisition microphone array；

Step 2, acquisition sound-source signal reach the arrival delay inequality of each array element of microphone array；

Step 3, construction microphone array apart from difference vector r；

Microphone array apart from difference vector r=[Δ r₂,…,Δr_i,…,Δr_N]^T, wherein, Δ r_i=c τ_iArrive for sound-source signal Reach the range difference of i-th array element of microphone array, τ_iThe arrival delay inequality of i-th array element is reached for sound-source signal, and c is passed for sound Broadcast speed, i=2 ..., N；

Step 4, the location matrix A of construction microphone array；

Wherein, x_n、y_n、z_nFor the position coordinates of n-th array element, n=1 ..., N；

The unit direction vector of sound-source signalThe unit direction of sound-source signal Vector and the relation between difference vector are：A Θ=r, wherein, θ is the azimuth of sound-source signal,Bowing for sound-source signal The elevation angle, []^TRepresent transposition operation；

Step 5, judge whether the array element of microphone array is generally aligned in the same plane, if it is execution step six, otherwise execute step Rapid seven；

Step 6, using least square method calculate sound-source signal two-dimensional directional vector, so as to obtain the azimuth of sound-source signal Estimate and angle of pitch estimate；

When the array element of microphone array is generally aligned in the same plane, the unit direction vector Θ of location matrix A and sound-source signal is moved back respectively Turn to two-dimensional position matrix A₂With two-dimensional directional vector Θ₂, according toThe azimuth for calculating sound-source signal is estimated EvaluationWith angle of pitch estimate

$\widehat{\mathrm{\θ}}=a\tan \left(\frac{{\mathrm{\μ}}_2}{{\mathrm{\μ}}_1}\right),$

μ in formula₁For two-dimensional directional vector Θ₂Section 1 element, μ₂For two-dimensional directional vector Θ₂Section 2 element；

Step 7, using Constraint least square algorithm calculate sound-source signal unit direction vector estimateSo as to obtain sound source The azimuth estimate and angle of pitch estimate of signal；

According toCalculate the azimuth estimate of sound-source signalWith angle of pitch estimate

μ in formula₁For unit direction vector estimateSection 1 element, μ₂For unit direction vector estimateSecond Item element, μ₃For unit direction vector estimateSection 3 element, I be unit matrix, λ for estimate parameter.

2. the microphone array direction-finding method based on least square method according to claim 1, it is characterised in that：The step It is whether zero judging whether array element is generally aligned in the same plane by checking each row of location matrix A in rapid 5, is classified as if any one Zero, then at grade, the array element of otherwise explanation microphone array is not distributed across for the array element distribution of explanation microphone array On same plane.

3. the microphone array direction-finding method based on least square method according to claim 1, it is characterised in that：The step In rapid 5 by check location matrix A whether sequency spectrum judging whether array element is generally aligned in the same plane, such as location matrix A is unsatisfactory for Sequency spectrum, then the array element of explanation microphone array be distributed at grade, the array element of otherwise explanation microphone array is not point Cloth is at grade.

4. the microphone array direction-finding method based on least square method according to claim 1, it is characterised in that：Described estimate Meter parameter lambda passes through equation r^TA(A^TA+λI)^-2A^TR=1 is solved, and calculation procedure is as follows：

To matrix A^TA carries out feature decomposition, i.e.,Wherein,

Λ=diag (γ₁,γ₂,γ₃) it is characterized value matrix, U=[u₁,u₂,u₃] corresponding feature matrix for it, N_dimFor The columns of location matrix A, N_dim=3, u_kIt is characterized value γ_kCorresponding characteristic vector, k=1,2,3, ()^HRepresent transposition again altogether Yoke is operated；

According to $r^{T}A\·U\stackrel{\mathrm{\Δ}}{=}{p}^T=[p_1,p_2,p_3]$ And $U^H\·A^{T}r\stackrel{\mathrm{\Δ}}{=}q=[q_1,q_2,q_3{]}^T,$ Obtain a_k=p_kq_k；

Estimate that parameter lambda passes through polynomial equation c₆λ⁶+c₅λ⁵+c₄λ⁴+c₃λ³+c₂λ²+c₁λ+c₀=0 is solved, polynomial equation In coefficient formulas as follows：

c₆=1,

c₅=2 γ₁+2γ₂+2γ₃,

$\begin{array}{c}{c}_4=({\mathrm{\γ}}_1^2+{\mathrm{\γ}}_2^2+{\mathrm{\γ}}_3^2+{4\mathrm{\γ}}_1{\mathrm{\γ}}_2+4{\mathrm{\γ}}_1{\mathrm{\γ}}_3+4{\mathrm{\γ}}_2{\mathrm{\γ}}_3)\\ -(a_1+a_2+a_3)\end{array},$

$\begin{array}{c}{c}_3=(2{\mathrm{\γ}}_1{\mathrm{\γ}}_2^2+2{\mathrm{\γ}}_1{\mathrm{\γ}}_3^2+2{\mathrm{\γ}}_1^2{\mathrm{\γ}}_2+2{\mathrm{\γ}}_1^2{\mathrm{\γ}}_3+2{\mathrm{\γ}}_2{\mathrm{\γ}}_3^2+2{\mathrm{\γ}}_2^2{\mathrm{\γ}}_3+8{\mathrm{\γ}}_1{\mathrm{\γ}}_2{\mathrm{\γ}}_3)\\ -[a_1(2{\mathrm{\γ}}_2+2{\mathrm{\γ}}_3)+a_2(2{\mathrm{\γ}}_1+2{\mathrm{\γ}}_3)+a_3(2{\mathrm{\γ}}_1+2{\mathrm{\γ}}_2)]\end{array},$

$\begin{array}{c}{c}_2=({\mathrm{\γ}}_1^2{\mathrm{\γ}}_2^2+{\mathrm{\γ}}_1^2{\mathrm{\γ}}_3^2+{\mathrm{\γ}}_2^2+{\mathrm{\γ}}_3^2+4{\mathrm{\γ}}_1{\mathrm{\γ}}_2^2{\mathrm{\γ}}_3+4{\mathrm{\γ}}_1^2{\mathrm{\γ}}_2{\mathrm{\γ}}_3)\\ -[a_1({\mathrm{\γ}}_2^2+{\mathrm{\γ}}_3^2+4{\mathrm{\γ}}_2{\mathrm{\γ}}_3)+a_2({\mathrm{\γ}}_1^2+{\mathrm{\γ}}_3^2+4{\mathrm{\γ}}_1{\mathrm{\γ}}_3)+a_3({\mathrm{\γ}}_1^2+{\mathrm{\γ}}_2^2+4{\mathrm{\γ}}_1{\mathrm{\γ}}_2)]\end{array},$

$\begin{array}{c}{c}_1=(2{\mathrm{\γ}}_1^2{\mathrm{\γ}}_2^2{\mathrm{\γ}}_3+2{\mathrm{\γ}}_1^2{\mathrm{\γ}}_2{\mathrm{\γ}}_3^2+2{\mathrm{\γ}}_1{\mathrm{\γ}}_2^2{\mathrm{\γ}}_3^2)\\ -[a_1(2{\mathrm{\γ}}_2{\mathrm{\γ}}_3^2+2{\mathrm{\γ}}_2^2{\mathrm{\γ}}_3)+a_2(2{\mathrm{\γ}}_1{\mathrm{\γ}}_3^2+2{\mathrm{\γ}}_1^2{\mathrm{\γ}}_3)+a_3(2{\mathrm{\γ}}_1{\mathrm{\γ}}_2^2+2{\mathrm{\γ}}_1^2{\mathrm{\γ}}_2)]\end{array},$

$c_0={\mathrm{\γ}}_1^2{\mathrm{\γ}}_2^2{\mathrm{\γ}}_3^2-[a_1\left({\mathrm{\γ}}_2^2{\mathrm{\γ}}_3^2\right)+a_2\left({\mathrm{\γ}}_1^2{\mathrm{\γ}}_3^2\right)+a_3\left({\mathrm{\γ}}_1^2{\mathrm{\γ}}_2^2\right)];$

Using polynomial rooting algorithm, equation root is tried to achieve according to above coefficient value, using root minimum for wherein modulus value as required Estimation parameter lambda.