Compute auto- and cross-correlations for complex sequences

Author: kagalenko-m-b

src/signalcorr.jl

@@ -14,7 +14,7 @@
 default_laglen(lx::Int) = min(lx-1, round(Int,10*log10(lx)))
 check_lags(lx::Int, lags::AbstractVector) = (maximum(lags) < lx || error("lags must be less than the sample length."))
 
-function demean_col!(z::RealVector, x::RealMatrix, j::Int, demean::Bool)
+function demean_col!(z::AbstractVector, x::AbstractMatrix, j::Int, demean::Bool)
     T = eltype(z)
     m = size(x, 1)
     @assert m == length(z)
@@ -43,8 +43,7 @@ end
 
 default_autolags(lx::Int) = 0 : default_laglen(lx)
 
-_autodot(x::AbstractVector{<:RealFP}, lx::Int, l::Int) = dot(x, 1:(lx-l), x, (1+l):lx)
-_autodot(x::RealVector, lx::Int, l::Int) = dot(view(x, 1:(lx-l)), view(x, (1+l):lx))
+_autodot(x::AbstractVector, lx::Int, l::Int) = dot(view(x, 1:(lx-l)), view(x, (1+l):lx))
 
 
 ## autocov
@@ -61,34 +60,32 @@ where each column in the result will correspond to a column in `x`.
 
 The output is not normalized. See [`autocor!`](@ref) for a method with normalization.
 """
-function autocov!(r::RealVector, x::RealVector, lags::IntegerVector; demean::Bool=true)
+function autocov!(
+    r::AbstractVector,
+    x::AbstractVector,
+    lags::AbstractVector{<:Integer},
+    z=zeros(typeof(zero(eltype(x)) / 1), length(x));
+    demean::Bool=true
+    )
     lx = length(x)
     m = length(lags)
     length(r) == m || throw(DimensionMismatch())
     check_lags(lx, lags)
-
     T = typeof(zero(eltype(x)) / 1)
-    z::Vector{T} = demean ? x .- mean(x) : x
+    z .= x
+    demean && z .= z .- mean(z)
     for k = 1 : m  # foreach lag value
         r[k] = _autodot(z, lx, lags[k]) / lx
     end
     return r
 end
 
-function autocov!(r::RealMatrix, x::RealMatrix, lags::IntegerVector; demean::Bool=true)
-    lx = size(x, 1)
-    ns = size(x, 2)
-    m = length(lags)
-    size(r) == (m, ns) || throw(DimensionMismatch())
-    check_lags(lx, lags)
-
-    T = typeof(zero(eltype(x)) / 1)
-    z = Vector{T}(undef, lx)
-    for j = 1 : ns
-        demean_col!(z, x, j, demean)
-        for k = 1 : m
-            r[k,j] = _autodot(z, lx, lags[k]) / lx
-        end
+function autocov!(
+    r::AbstractMatrix, x::AbstractMatrix, lags::AbstractVector{<:Integer}; demean::Bool=true
+    )
+    z = zeros(typeof(zero(eltype(x))/1), size(x, 1))
+    for n in 1:size(x, 2)
+        autocov!(view(r, :, n), view(x, :, n), lags, z; demean)
     end
     return r
 end
@@ -110,18 +107,18 @@ When left unspecified, the lags used are the integers from 0 to
 
 The output is not normalized. See [`autocor`](@ref) for a function with normalization.
 """
-function autocov(x::RealVector, lags::IntegerVector; demean::Bool=true)
+function autocov(x::AbstractVector, lags::AbstractVector{<:Integer}; demean::Bool=true)
     out = Vector{float(eltype(x))}(undef, length(lags))
-    autocov!(out, x, lags; demean=demean)
+    autocov!(out, x, lags; demean)
 end
 
-function autocov(x::RealMatrix, lags::IntegerVector; demean::Bool=true)
+function autocov(x::AbstractMatrix, lags::AbstractVector{<:Integer}; demean::Bool=true)
     out = Matrix{float(eltype(x))}(undef, length(lags), size(x,2))
-    autocov!(out, x, lags; demean=demean)
+    autocov!(out, x, lags; demean)
 end
 
-autocov(x::AbstractVecOrMat{<:Real}; demean::Bool=true) =
-    autocov(x, default_autolags(size(x,1)); demean=demean)
+autocov(x::AbstractVecOrMat; demean::Bool=true) =
+    autocov(x, default_autolags(size(x,1)); demean)
 
 ## autocor
 
@@ -139,36 +136,27 @@ where each column in the result will correspond to a column in `x`.
 The output is normalized by the variance of `x`, i.e. so that the lag 0
 autocorrelation is 1. See [`autocov!`](@ref) for the unnormalized form.
 """
-function autocor!(r::RealVector, x::RealVector, lags::IntegerVector; demean::Bool=true)
-    lx = length(x)
-    m = length(lags)
-    length(r) == m || throw(DimensionMismatch())
-    check_lags(lx, lags)
-
-    T = typeof(zero(eltype(x)) / 1)
-    z::Vector{T} = demean ? x .- mean(x) : x
-    zz = dot(z, z)
-    for k = 1 : m  # foreach lag value
-        r[k] = _autodot(z, lx, lags[k]) / zz
-    end
+function autocor!(
+    r::AbstractVector,
+    x::AbstractVector,
+    lags::AbstractVector{<:Integer},
+    z=zeros(typeof(zero(eltype(x)) / 1), length(x));
+    demean::Bool=true
+    )
+    autocov!(view(r, 1:1), x, 0:0, z; demean)
+    zz = r[1]
+    autocov!(r, x, lags, z; demean)
+    r ./= zz
     return r
 end
 
-function autocor!(r::RealMatrix, x::RealMatrix, lags::IntegerVector; demean::Bool=true)
-    lx = size(x, 1)
-    ns = size(x, 2)
-    m = length(lags)
-    size(r) == (m, ns) || throw(DimensionMismatch())
-    check_lags(lx, lags)
-
-    T = typeof(zero(eltype(x)) / 1)
-    z = Vector{T}(undef, lx)
-    for j = 1 : ns
-        demean_col!(z, x, j, demean)
-        zz = dot(z, z)
-        for k = 1 : m
-            r[k,j] = _autodot(z, lx, lags[k]) / zz
-        end
+function autocor!(
+    r::AbstractMatrix, x::AbstractMatrix, lags::AbstractVector{<:Integer}; demean::Bool=true
+    )
+    T = typeof(zero(eltype(x))/1)
+    z = Vector{T}(undef, size(x, 1))
+    for n in 1:size(x, 2)
+        autocor!(view(r, :, n), view(x, :, n), lags, z; demean)
     end
     return r
 end
@@ -191,18 +179,18 @@ When left unspecified, the lags used are the integers from 0 to
 The output is normalized by the variance of `x`, i.e. so that the lag 0
 autocorrelation is 1. See [`autocov`](@ref) for the unnormalized form.
 """
-function autocor(x::RealVector, lags::IntegerVector; demean::Bool=true)
+function autocor(x::AbstractVector, lags::AbstractVector{<:Integer}; demean::Bool=true)
     out = Vector{float(eltype(x))}(undef, length(lags))
-    autocor!(out, x, lags; demean=demean)
+    autocor!(out, x, lags; demean)
 end
 
-function autocor(x::RealMatrix, lags::IntegerVector; demean::Bool=true)
+function autocor(x::AbstractMatrix, lags::AbstractVector{<:Integer}; demean::Bool=true)
     out = Matrix{float(eltype(x))}(undef, length(lags), size(x,2))
-    autocor!(out, x, lags; demean=demean)
+    autocor!(out, x, lags; demean)
 end
 
-autocor(x::AbstractVecOrMat{<:Real}; demean::Bool=true) =
-    autocor(x, default_autolags(size(x,1)); demean=demean)
+autocor(x::AbstractVecOrMat; demean::Bool=true) =
+    autocor(x, default_autolags(size(x,1)); demean)
 
 
 #######################################
@@ -213,14 +201,7 @@ autocor(x::AbstractVecOrMat{<:Real}; demean::Bool=true) =
 
 default_crosslags(lx::Int) = (l=default_laglen(lx); -l:l)
 
-function _crossdot(x::AbstractVector{T}, y::AbstractVector{T}, lx::Int, l::Int) where {T<:RealFP}
-    if l >= 0
-        dot(x, 1:(lx-l), y, (1+l):lx)
-    else
-        dot(x, (1-l):lx, y, 1:(lx+l))
-    end
-end
-function _crossdot(x::RealVector, y::RealVector, lx::Int, l::Int)
+function _crossdot(x::AbstractVector, y::AbstractVector, lx::Int, l::Int)
     if l >= 0
         dot(view(x, 1:(lx-l)), view(y, (1+l):lx))
     else
@@ -246,7 +227,7 @@ three-dimensional array of size `(length(lags), size(x, 2), size(y, 2))`.
 
 The output is not normalized. See [`crosscor!`](@ref) for a function with normalization.
 """
-function crosscov!(r::RealVector, x::RealVector, y::RealVector, lags::IntegerVector; demean::Bool=true)
+function crosscov!(r::AbstractVector, x::AbstractVector, y::AbstractVector, lags::AbstractVector{<:Integer}; demean::Bool=true)
     lx = length(x)
     m = length(lags)
     (length(y) == lx && length(r) == m) || throw(DimensionMismatch())
@@ -262,7 +243,7 @@ function crosscov!(r::RealVector, x::RealVector, y::RealVector, lags::IntegerVec
     return r
 end
 
-function crosscov!(r::RealMatrix, x::RealMatrix, y::RealVector, lags::IntegerVector; demean::Bool=true)
+function crosscov!(r::AbstractMatrix, x::AbstractMatrix, y::AbstractVector, lags::AbstractVector{<:Integer}; demean::Bool=true)
     lx = size(x, 1)
     ns = size(x, 2)
     m = length(lags)
@@ -282,7 +263,7 @@ function crosscov!(r::RealMatrix, x::RealMatrix, y::RealVector, lags::IntegerVec
     return r
 end
 
-function crosscov!(r::RealMatrix, x::RealVector, y::RealMatrix, lags::IntegerVector; demean::Bool=true)
+function crosscov!(r::AbstractMatrix, x::AbstractVector, y::AbstractMatrix, lags::AbstractVector{<:Integer}; demean::Bool=true)
     lx = length(x)
     ns = size(y, 2)
     m = length(lags)
@@ -302,7 +283,7 @@ function crosscov!(r::RealMatrix, x::RealVector, y::RealMatrix, lags::IntegerVec
     return r
 end
 
-function crosscov!(r::AbstractArray{<:Real,3}, x::RealMatrix, y::RealMatrix, lags::IntegerVector; demean::Bool=true)
+function crosscov!(r::AbstractArray{U,3} where U, x::AbstractMatrix, y::AbstractMatrix, lags::AbstractVector{<:Integer}; demean::Bool=true)
     lx = size(x, 1)
     nx = size(x, 2)
     ny = size(y, 2)
@@ -356,27 +337,27 @@ When left unspecified, the lags used are the integers from
 
 The output is not normalized. See [`crosscor`](@ref) for a function with normalization.
 """
-function crosscov(x::RealVector, y::RealVector, lags::IntegerVector; demean::Bool=true)
+function crosscov(x::AbstractVector, y::AbstractVector, lags::AbstractVector{<:Integer}; demean::Bool=true)
     out = Vector{float(Base.promote_eltype(x, y))}(undef, length(lags))
     crosscov!(out, x, y, lags; demean=demean)
 end
 
-function crosscov(x::RealMatrix, y::RealVector, lags::IntegerVector; demean::Bool=true)
+function crosscov(x::AbstractMatrix, y::AbstractVector, lags::AbstractVector{<:Integer}; demean::Bool=true)
     out = Matrix{float(Base.promote_eltype(x, y))}(undef, length(lags), size(x,2))
     crosscov!(out, x, y, lags; demean=demean)
 end
 
-function crosscov(x::RealVector, y::RealMatrix, lags::IntegerVector; demean::Bool=true)
+function crosscov(x::AbstractVector, y::AbstractMatrix, lags::AbstractVector{<:Integer}; demean::Bool=true)
     out = Matrix{float(Base.promote_eltype(x, y))}(undef, length(lags), size(y,2))
     crosscov!(out, x, y, lags; demean=demean)
 end
 
-function crosscov(x::RealMatrix, y::RealMatrix, lags::IntegerVector; demean::Bool=true)
+function crosscov(x::AbstractMatrix, y::AbstractMatrix, lags::AbstractVector{<:Integer}; demean::Bool=true)
     out = Array{float(Base.promote_eltype(x, y)),3}(undef, length(lags), size(x,2), size(y,2))
     crosscov!(out, x, y, lags; demean=demean)
 end
 
-crosscov(x::AbstractVecOrMat{<:Real}, y::AbstractVecOrMat{<:Real}; demean::Bool=true) =
+crosscov(x::AbstractVecOrMat, y::AbstractVecOrMat; demean::Bool=true) =
     crosscov(x, y, default_crosslags(size(x,1)); demean=demean)
 
 
@@ -397,7 +378,7 @@ three-dimensional array of size `(length(lags), size(x, 2), size(y, 2))`.
 The output is normalized by `sqrt(var(x)*var(y))`. See [`crosscov!`](@ref) for the
 unnormalized form.
 """
-function crosscor!(r::RealVector, x::RealVector, y::RealVector, lags::IntegerVector; demean::Bool=true)
+function crosscor!(r::AbstractVector, x::AbstractVector, y::AbstractVector, lags::AbstractVector{<:Integer}; demean::Bool=true)
     lx = length(x)
     m = length(lags)
     (length(y) == lx && length(r) == m) || throw(DimensionMismatch())
@@ -414,7 +395,7 @@ function crosscor!(r::RealVector, x::RealVector, y::RealVector, lags::IntegerVec
     return r
 end
 
-function crosscor!(r::RealMatrix, x::RealMatrix, y::RealVector, lags::IntegerVector; demean::Bool=true)
+function crosscor!(r::AbstractMatrix, x::AbstractMatrix, y::AbstractVector, lags::AbstractVector{<:Integer}; demean::Bool=true)
     lx = size(x, 1)
     ns = size(x, 2)
     m = length(lags)
@@ -436,7 +417,7 @@ function crosscor!(r::RealMatrix, x::RealMatrix, y::RealVector, lags::IntegerVec
     return r
 end
 
-function crosscor!(r::RealMatrix, x::RealVector, y::RealMatrix, lags::IntegerVector; demean::Bool=true)
+function crosscor!(r::AbstractMatrix, x::AbstractVector, y::AbstractMatrix, lags::AbstractVector{<:Integer}; demean::Bool=true)
     lx = length(x)
     ns = size(y, 2)
     m = length(lags)
@@ -458,7 +439,7 @@ function crosscor!(r::RealMatrix, x::RealVector, y::RealMatrix, lags::IntegerVec
     return r
 end
 
-function crosscor!(r::AbstractArray{<:Real,3}, x::RealMatrix, y::RealMatrix, lags::IntegerVector; demean::Bool=true)
+function crosscor!(r::AbstractArray{U,3} where U, x::AbstractMatrix, y::AbstractMatrix, lags::AbstractVector{<:Integer}; demean::Bool=true)
     lx = size(x, 1)
     nx = size(x, 2)
     ny = size(y, 2)
@@ -517,27 +498,27 @@ When left unspecified, the lags used are the integers from
 The output is normalized by `sqrt(var(x)*var(y))`. See [`crosscov`](@ref) for the
 unnormalized form.
 """
-function crosscor(x::RealVector, y::RealVector, lags::IntegerVector; demean::Bool=true)
+function crosscor(x::AbstractVector, y::AbstractVector, lags::AbstractVector{<:Integer}; demean::Bool=true)
     out = Vector{float(Base.promote_eltype(x, y))}(undef, length(lags))
     crosscor!(out, x, y, lags; demean=demean)
 end
 
-function crosscor(x::RealMatrix, y::RealVector, lags::IntegerVector; demean::Bool=true)
+function crosscor(x::AbstractMatrix, y::AbstractVector, lags::AbstractVector{<:Integer}; demean::Bool=true)
     out = Matrix{float(Base.promote_eltype(x, y))}(undef, length(lags), size(x,2))
     crosscor!(out, x, y, lags; demean=demean)
 end
 
-function crosscor(x::RealVector, y::RealMatrix, lags::IntegerVector; demean::Bool=true)
+function crosscor(x::AbstractVector, y::AbstractMatrix, lags::AbstractVector{<:Integer}; demean::Bool=true)
     out = Matrix{float(Base.promote_eltype(x, y))}(undef, length(lags), size(y,2))
     crosscor!(out, x, y, lags; demean=demean)
 end
 
-function crosscor(x::RealMatrix, y::RealMatrix, lags::IntegerVector; demean::Bool=true)
+function crosscor(x::AbstractMatrix, y::AbstractMatrix, lags::AbstractVector{<:Integer}; demean::Bool=true)
     out = Array{float(Base.promote_eltype(x, y)),3}(undef, length(lags), size(x,2), size(y,2))
     crosscor!(out, x, y, lags; demean=demean)
 end
 
-crosscor(x::AbstractVecOrMat{<:Real}, y::AbstractVecOrMat{<:Real}; demean::Bool=true) =
+crosscor(x::AbstractVecOrMat, y::AbstractVecOrMat; demean::Bool=true) =
     crosscor(x, y, default_crosslags(size(x,1)); demean=demean)
 
 

test/rcall_test.jl

@@ -0,0 +1,29 @@
+using RCall
+
+_r_ccf(x, y, lm) = rcopy(rcall(:ccf, x, y, plot=false, lag=lm, type="covariance"))[:acf]
+_r_acf(x, lm) = rcopy(rcall(:acf, x, plot=false, lag=lm, type="covariance"))[:acf]
+_r_pacf(x, lm) = rcopy(rcall(:acf, x, plot=false, lag=lm))[:acf]
+
+function _ccf(is_auto::Bool, x::AbstractVector{<:Real}...)
+    lx = size(x[1], 1)
+    lag_max =  min(lx - 1, round(Int, 10*log10(lx)))
+    cor = reverse(dropdims(_r_ccf(x[1], x[end], lag_max), dims = (2, 3)))
+    is_auto ? cor[lx:end] : cor
+end
+
+r_autocov(x::AbstractVector) = r_crosscov(x, x, true)
+
+function r_crosscov(x::AbstractVector, y::AbstractVector, is_auto=false)
+    cv(x,y) = _ccf(is_auto, x, y)
+    cc_re = cv(real.(x), real.(y)) + cv(imag.(x), imag.(y))
+    cc_im = cv(real.(x), imag.(y)) - cv(imag.(x), real.(y))
+    cc_re + im*cc_im
+end
+
+r_autocor(x::AbstractVector) = r_crosscor(x, x, true)
+
+function r_crosscor(x::AbstractVector, y::AbstractVector, is_auto=false)
+    cc = r_crosscov(x, y, is_auto)
+    sc(z) = _r_ccf(real.(z), real.(z), 0) + _r_ccf(imag.(z), imag.(z), 0)
+    cc/sqrt(sc(x)*sc(y))
+end