Optimized the use of Durbin algorithm

Author: kagalenko-m-b

src/signalcorr.jl

@@ -549,12 +549,14 @@ function pacf_regress!(r::AbstractMatrix, X::AbstractMatrix, lags::IntegerVector
 end
 
 function pacf_yulewalker!(r::AbstractMatrix, X::AbstractMatrix, lags::IntegerVector, mk::Integer)
-    tmp = Vector{eltype(X)}(undef, mk)
+    p = Vector{eltype(X)}(undef, mk)
+    y = Vector{eltype(X)}(undef, mk)
     for j = 1 : size(X,2)
         acfs = autocor(X[:,j], 1:mk)
+        tmp = durbin!(acfs, p, y)
         for i = 1 : length(lags)
             l = lags[i]
-            r[i,j] = l == 0 ? 1 : l == 1 ? acfs[i] : -durbin!(view(acfs, 1:l), tmp)[l]
+            r[i,j] = l == 0 ? 1 : -p[l]
         end
     end
 end

src/toeplitzsolvers.jl

@@ -1,8 +1,9 @@
 # Symmetric Toeplitz solver
-function durbin!(r::AbstractVector{T}, y::AbstractVector{T}) where T
+function durbin!(r::AbstractVector{T}, p::AbstractVector{T}, y::AbstractVector{T}) where T
     n = length(r)
-    n <= length(y) || throw(DimensionMismatch("Auxiliary vector cannot be shorter than data vector"))
+    n <= length(p) || n <= length(y) || throw(DimensionMismatch("Auxiliary vectors cannot be shorter than data vector"))
     y[1] = -r[1]
+    p[1] = -r[1]
     β = one(T)
     α = -r[1]
     for k = 1:n-1
@@ -19,10 +20,11 @@ function durbin!(r::AbstractVector{T}, y::AbstractVector{T}) where T
         end
         if isodd(k) y[div(k,2)+1] += α*conj(y[div(k,2)+1]) end
         y[k+1] = α
+        p[k+1] = α
     end
-    return y
+    return p
 end
-durbin(r::AbstractVector{T}) where T = durbin!(r, zeros(T, length(r)))
+durbin(r::AbstractVector{T}) where T = durbin!(r, zeros(T, length(r)), zeros(T, length(r)))
 
 function levinson!(r::AbstractVector{T}, b::AbstractVector{T}, x::AbstractVector{T}) where T<:BlasReal
     n = length(b)

test/signalcorr.jl

@@ -103,18 +103,35 @@ c22 = crosscor(x2, x2)
 @test crosscor(x,  x)  ≈ cat([c11 c21], [c12 c22], dims=3)
 
 
-## pacf
-
-pacfr = [-1.598495044296996e-03 - 2.915104118351207e-01im
-         -5.560162016912027e-01 + 2.950837739894279e-01im
-         -2.547001916363494e-02 + 2.326084658014266e-01im
-         -5.427443903358727e-01 + 3.146715147305132e-01im];
-
-@test pacf(x1, 1:4) ≈ pacfr[1:4]
-
-pacfy = [-1.598495044296996e-03 - 2.915104118351207e-01im
-         -5.560162016912027e-01 + 2.950837739894279e-01im
-         -2.547001916363494e-02 + 2.326084658014266e-01im
-         -5.427443903358727e-01 + 3.146715147305132e-01im];
-
-@test pacf(x1, 1:4, method=:yulewalker) ≈ pacfy
+## pacf least squares
+pacf_ls = [-1.598495044296996e-03 - 2.915104118351207e-01im
+           -5.560162016912027e-01 + 2.950837739894279e-01im
+           -2.547001916363494e-02 + 2.326084658014266e-01im
+           -5.427443903358727e-01 + 3.146715147305132e-01im]
+
+@test pacf(x1, 1:4) ≈ pacf_ls[1:4]
+
+## pacf Yule-Walker
+
+function yulewalker_qr(v::AbstractVector)
+    A = toeplitz(v)
+    b = v[2:end]
+    x =-A\b
+end
+function toeplitz(v::AbstractVector{T}) where T
+    N=length(v)
+    A = zeros(T, N - 1, N - 1)
+    for n in 1:N-1
+        A[n, n + 1:end] = conj(v[2:N-n])
+        A[n, 1:n] = reverse(v[1:n])
+    end
+    return A
+end
+# durbin solver 
+acf = autocor(x1)
+p = [yulewalker_qr(acf[1:n])[n-1] for n in 2:length(acf)]
+@test p ≈ StatsBase.durbin(acf[2:end])
+
+pacfy = [];
+
+@test pacf(x1, 1:4, method=:yulewalker) ≈ -p[1:4]