Wilson-Burg spectral factorization

Spectral factorization constructs a minimum-phase signal from its spectrum. The algorithm, suggested by Wilson (1969), approaches this problem directly with Newton's iterative method. In a Z-transform notation, Wilson's method implies solving the equation  

$$S(Z) = A(Z)\bar{A}(1/Z)$$ (79)

for a given spectrum S(Z) and unknown minimum-phase signal A(Z) with an iterative linearization

$$\begin{eqnarray} \nonumber S(Z) & = & A_t(Z)\bar A_t(1/Z)+ A_t( Z)[\bar A_{t+1}... ...bar A_{t+1}(1/Z) + \bar A_t(1/Z) A_{t+1} - A_t(Z)\bar A_t(1/Z) \;,\end{eqnarray}$$ (80)

where A_t(Z) denotes the signal estimate at iteration t. Starting from some initial estimate A₀(Z), such as A₀(Z)=1, one iteratively solves the linear equation () for the updated signal A_t+1(Z). Wilson 1969 presents a rigorous proof that iteration () operates with minimum-phase signals provided that the initial estimate A₀(Z) is minimum-phase.

Burg (1998, personal communication) recognized that dividing both sides of equation () by $\bar A_t(1/Z) A_t(Z)$ leads to a particularly convenient form, where the terms on the left are symmetric, and the two terms on the right are correspondingly strictly causal and anticausal:  

$$1 \ +\ {S(Z) \over \bar A_t(1/Z)\ A_t(Z)} = {A_{t+1}(Z) \over A_t(Z)} \ +\ {\bar A_{t+1}(1/Z) \over \bar A_t(1/Z)}$$ (81)

Equation () leads to the Wilson-Burg algorithm, which accomplishes spectral factorization by a recursive application of convolution (polynomial multiplication) and deconvolution (polynomial division). The algorithm proceeds as follows:

1.

Compute the left side of equation () using forward and adjoint polynomial division.

2.

Abandon negative lags, to keep only the causal part of the signal, and also keep half of the zero lag. This gives us A_t+1(Z)/A_t(Z).

3.

Multiply out (convolve) the denominator A_t(Z). Now we have the desired result A_t+1(Z).

4.

Iterate until convergence.

 

iter a₀ a₁ a₂ a₃

1.000000 0.000000 0.000000 0.000000

1 36.523964 23.737839 6.625787 0.657103

2 26.243151 25.726116 8.471050 0.914951

3 24.162354 25.991493 8.962727 0.990802

4 24.001223 25.999662 9.000164 0.999200

5 24.000015 25.999977 9.000029 0.999944

6 23.999998 26.000002 9.000003 0.999996

7 23.999998 26.000004 9.000001 1.000000

8 23.999998 25.999998 9.000000 1.000000

9 24.000000 26.000000 9.000000 1.000000

An example of the Wilson-Burg convergence is shown in Table  on a simple 1-D signal. The autocorrelation S(Z) in this case is $1334 + 867 \left(Z + 1/Z\right) + 242 \left(Z^2 + 1/Z^2\right) + 24 \left(Z^3 + 1/Z^3\right)$, and the corresponding minimum-phase signal is A(Z) = (2+Z)(3+Z)(4+Z) = 24 + 26 Z + 9 Z² + Z³. A quadratic rate of convergence is visible from the table. The convergence slows down for signals whose polynomial roots are close to the unit circle Wilson (1969).