Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Complex signalcorr #742

Open
wants to merge 18 commits into
base: master
Choose a base branch
from
Open

Complex signalcorr #742

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
18 commits
Select commit Hold shift + click to select a range
32a14ae
Get rid of old type aliases (#840)
Dec 6, 2022
e58f9ca
Merge branch 'JuliaStats:master' into master
kagalenko-m-b Feb 14, 2023
3747453
Compute auto- and cross-correlations for complex sequences
kagalenko-m-b Nov 27, 2021
b62800c
Durbin algorithm works correctly
kagalenko-m-b Nov 29, 2021
7103178
Optimized the use of Durbin algorithm
kagalenko-m-b Nov 30, 2021
d32b14d
Reduced allocations in autocov!()
kagalenko-m-b Nov 30, 2021
97d87d3
Update test/signalcorr.jl
kagalenko-m-b Nov 30, 2021
1fff4e7
Apply suggestions from code review
kagalenko-m-b Nov 30, 2021
57b0fab
Add documentation to durbin!() and restore tests
kagalenko-m-b Nov 30, 2021
c06c4d7
Simplify subtracting mean, as suggested by reviewer
kagalenko-m-b Dec 1, 2021
dac79be
Moved new tests with complex data into a separate file
kagalenko-m-b Dec 2, 2021
3384cb1
Apply suggestions from code review
kagalenko-m-b Dec 3, 2021
9f2c943
Limit changes to more general argument types
kagalenko-m-b Dec 4, 2021
cb6ea64
Apply suggestions from code review
Dec 5, 2021
c135edf
Correct the documentation to reflect more general input type
kagalenko-m-b Dec 16, 2021
4f30211
Remove unneeded stuff from tests
May 29, 2022
0bfa47f
Merge branch 'JuliaStats:master' into complex_signalcorr
kagalenko-m-b Mar 17, 2023
e4930f8
Merge branch 'JuliaStats:master' into complex_signalcorr
kagalenko-m-b Mar 28, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
93 changes: 48 additions & 45 deletions src/signalcorr.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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::AbstractVector{<:Real}, x::AbstractMatrix{<:Real}, 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)
Expand Down Expand Up @@ -44,7 +44,7 @@ end
default_autolags(lx::Int) = 0 : default_laglen(lx)

_autodot(x::AbstractVector{<:Union{Float32, Float64}}, lx::Int, l::Int) = dot(x, 1:(lx-l), x, (1+l):lx)
_autodot(x::AbstractVector{<:Real}, 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))
kagalenko-m-b marked this conversation as resolved.
Show resolved Hide resolved


## autocov
Expand All @@ -61,7 +61,7 @@ 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::AbstractVector{<:Real}, x::AbstractVector{<:Real}, lags::AbstractVector{<:Integer}; demean::Bool=true)
function autocov!(r::AbstractVector, x::AbstractVector, lags::AbstractVector{<:Integer}; demean::Bool=true)
lx = length(x)
m = length(lags)
length(r) == m || throw(DimensionMismatch())
kagalenko-m-b marked this conversation as resolved.
Show resolved Hide resolved
Expand All @@ -75,7 +75,7 @@ function autocov!(r::AbstractVector{<:Real}, x::AbstractVector{<:Real}, lags::Ab
return r
end

function autocov!(r::AbstractMatrix{<:Real}, x::AbstractMatrix{<:Real}, lags::AbstractVector{<:Integer}; demean::Bool=true)
function autocov!(r::AbstractMatrix, x::AbstractMatrix, lags::AbstractVector{<:Integer}; demean::Bool=true)
lx = size(x, 1)
ns = size(x, 2)
m = length(lags)
Expand Down Expand Up @@ -110,17 +110,17 @@ 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::AbstractVector{<:Real}, lags::AbstractVector{<:Integer}; 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)
end

function autocov(x::AbstractMatrix{<:Real}, lags::AbstractVector{<:Integer}; 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)
end

autocov(x::AbstractVecOrMat{<:Real}; demean::Bool=true) =
autocov(x::AbstractVecOrMat; demean::Bool=true) =
autocov(x, default_autolags(size(x,1)); demean=demean)

## autocor
Expand All @@ -139,7 +139,7 @@ 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::AbstractVector{<:Real}, x::AbstractVector{<:Real}, lags::AbstractVector{<:Integer}; demean::Bool=true)
function autocor!(r::AbstractVector, x::AbstractVector, lags::AbstractVector{<:Integer}; demean::Bool=true)
lx = length(x)
m = length(lags)
length(r) == m || throw(DimensionMismatch())
Expand All @@ -154,7 +154,7 @@ function autocor!(r::AbstractVector{<:Real}, x::AbstractVector{<:Real}, lags::Ab
return r
end

function autocor!(r::AbstractMatrix{<:Real}, x::AbstractMatrix{<:Real}, lags::AbstractVector{<:Integer}; demean::Bool=true)
function autocor!(r::AbstractMatrix, x::AbstractMatrix, lags::AbstractVector{<:Integer}; demean::Bool=true)
lx = size(x, 1)
ns = size(x, 2)
m = length(lags)
Expand Down Expand Up @@ -191,17 +191,17 @@ 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::AbstractVector{<:Real}, lags::AbstractVector{<:Integer}; 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)
end

function autocor(x::AbstractMatrix{<:Real}, lags::AbstractVector{<:Integer}; 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)
end

autocor(x::AbstractVecOrMat{<:Real}; demean::Bool=true) =
autocor(x::AbstractVecOrMat; demean::Bool=true) =
autocor(x, default_autolags(size(x,1)); demean=demean)


Expand All @@ -213,14 +213,14 @@ 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<:Union{Float32, Float64}}
function _crossdot(x::AbstractVector, y::AbstractVector, lx::Int, l::Int)
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::AbstractVector{<:Real}, y::AbstractVector{<:Real}, 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
Expand All @@ -233,7 +233,7 @@ end
"""
crosscov!(r, x, y, lags; demean=true)

Compute the cross covariance function (CCF) between real-valued vectors or matrices
Compute the cross covariance function (CCF) between vectors or matrices
`x` and `y` at `lags` and store the result in `r`. `demean` specifies whether the
respective means of `x` and `y` should be subtracted from them before computing their
CCF.
Expand All @@ -246,7 +246,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::AbstractVector{<:Real}, x::AbstractVector{<:Real}, y::AbstractVector{<:Real}, lags::AbstractVector{<:Integer}; 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())
Expand All @@ -262,7 +262,7 @@ function crosscov!(r::AbstractVector{<:Real}, x::AbstractVector{<:Real}, y::Abst
return r
end

function crosscov!(r::AbstractMatrix{<:Real}, x::AbstractMatrix{<:Real}, y::AbstractVector{<:Real}, lags::AbstractVector{<:Integer}; 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)
Expand All @@ -282,7 +282,7 @@ function crosscov!(r::AbstractMatrix{<:Real}, x::AbstractMatrix{<:Real}, y::Abst
return r
end

function crosscov!(r::AbstractMatrix{<:Real}, x::AbstractVector{<:Real}, y::AbstractMatrix{<:Real}, lags::AbstractVector{<:Integer}; 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)
Expand All @@ -302,7 +302,7 @@ function crosscov!(r::AbstractMatrix{<:Real}, x::AbstractVector{<:Real}, y::Abst
return r
end

function crosscov!(r::AbstractArray{<:Real,3}, x::AbstractMatrix{<:Real}, y::AbstractMatrix{<:Real}, lags::AbstractVector{<:Integer}; demean::Bool=true)
function crosscov!(r::AbstractArray{<:Any,3}, x::AbstractMatrix, y::AbstractMatrix, lags::AbstractVector{<:Integer}; demean::Bool=true)
lx = size(x, 1)
nx = size(x, 2)
ny = size(y, 2)
Expand Down Expand Up @@ -343,7 +343,7 @@ end
"""
crosscov(x, y, [lags]; demean=true)

Compute the cross covariance function (CCF) between real-valued vectors or
Compute the cross covariance function (CCF) between vectors or
matrices `x` and `y`, optionally specifying the `lags`. `demean` specifies
whether the respective means of `x` and `y` should be subtracted from them
before computing their CCF.
Expand All @@ -356,35 +356,35 @@ 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::AbstractVector{<:Real}, y::AbstractVector{<:Real}, lags::AbstractVector{<:Integer}; 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::AbstractMatrix{<:Real}, y::AbstractVector{<:Real}, lags::AbstractVector{<:Integer}; 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::AbstractVector{<:Real}, y::AbstractMatrix{<:Real}, lags::AbstractVector{<:Integer}; 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::AbstractMatrix{<:Real}, y::AbstractMatrix{<:Real}, lags::AbstractVector{<:Integer}; 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)


## crosscor
"""
crosscor!(r, x, y, lags; demean=true)

Compute the cross correlation between real-valued vectors or matrices `x` and `y` at
Compute the cross correlation between vectors or matrices `x` and `y` at
`lags` and store the result in `r`. `demean` specifies whether the respective means of
`x` and `y` should be subtracted from them before computing their cross correlation.

Expand All @@ -397,7 +397,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::AbstractVector{<:Real}, x::AbstractVector{<:Real}, y::AbstractVector{<:Real}, lags::AbstractVector{<:Integer}; 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())
Expand All @@ -414,7 +414,7 @@ function crosscor!(r::AbstractVector{<:Real}, x::AbstractVector{<:Real}, y::Abst
return r
end

function crosscor!(r::AbstractMatrix{<:Real}, x::AbstractMatrix{<:Real}, y::AbstractVector{<:Real}, lags::AbstractVector{<:Integer}; 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)
Expand All @@ -436,7 +436,7 @@ function crosscor!(r::AbstractMatrix{<:Real}, x::AbstractMatrix{<:Real}, y::Abst
return r
end

function crosscor!(r::AbstractMatrix{<:Real}, x::AbstractVector{<:Real}, y::AbstractMatrix{<:Real}, lags::AbstractVector{<:Integer}; 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)
Expand All @@ -458,7 +458,7 @@ function crosscor!(r::AbstractMatrix{<:Real}, x::AbstractVector{<:Real}, y::Abst
return r
end

function crosscor!(r::AbstractArray{<:Real,3}, x::AbstractMatrix{<:Real}, y::AbstractMatrix{<:Real}, lags::AbstractVector{<:Integer}; demean::Bool=true)
function crosscor!(r::AbstractArray{<:Any, 3}, x::AbstractMatrix, y::AbstractMatrix, lags::AbstractVector{<:Integer}; demean::Bool=true)
lx = size(x, 1)
nx = size(x, 2)
ny = size(y, 2)
Expand Down Expand Up @@ -504,7 +504,7 @@ end
"""
crosscor(x, y, [lags]; demean=true)

Compute the cross correlation between real-valued vectors or matrices `x` and `y`,
Compute the cross correlation between vectors or matrices `x` and `y`,
optionally specifying the `lags`. `demean` specifies whether the respective means of
`x` and `y` should be subtracted from them before computing their cross correlation.

Expand All @@ -517,27 +517,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::AbstractVector{<:Real}, y::AbstractVector{<:Real}, lags::AbstractVector{<:Integer}; 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::AbstractMatrix{<:Real}, y::AbstractVector{<:Real}, lags::AbstractVector{<:Integer}; 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::AbstractVector{<:Real}, y::AbstractMatrix{<:Real}, lags::AbstractVector{<:Integer}; 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::AbstractMatrix{<:Real}, y::AbstractMatrix{<:Real}, lags::AbstractVector{<:Integer}; 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)


Expand All @@ -549,7 +549,7 @@ crosscor(x::AbstractVecOrMat{<:Real}, y::AbstractVecOrMat{<:Real}; demean::Bool=
#
#######################################

function pacf_regress!(r::AbstractMatrix{<:Real}, X::AbstractMatrix{<:Real}, lags::AbstractVector{<:Integer}, mk::Integer)
function pacf_regress!(r::AbstractMatrix, X::AbstractMatrix, lags::AbstractVector{<:Integer}, mk::Integer)
lx = size(X, 1)
tmpX = ones(eltype(X), lx, mk + 1)
for j = 1 : size(X,2)
Expand All @@ -561,19 +561,22 @@ function pacf_regress!(r::AbstractMatrix{<:Real}, X::AbstractMatrix{<:Real}, lag
for i = 1 : length(lags)
l = lags[i]
sX = view(tmpX, 1+l:lx, 1:l+1)
r[i,j] = l == 0 ? 1 : (cholesky!(sX'sX, Val(false)) \ (sX'view(X, 1+l:lx, j)))[end]
r[i,j] = l == 0 ? 1 : (cholesky!(sX'sX) \ (sX'view(X, 1+l:lx, j)))[end]
end
end
r
end

function pacf_yulewalker!(r::AbstractMatrix{<:Real}, X::AbstractMatrix{T}, lags::AbstractVector{<:Integer}, mk::Integer) where T<:Union{Float32, Float64}
tmp = Vector{T}(undef, mk)
function pacf_yulewalker!(r::AbstractMatrix, X::AbstractMatrix, lags::AbstractVector{<:Integer}, mk::Integer)
T = eltype(X)
p = Vector{T}(undef, mk)
y = Vector{T}(undef, mk)
for j = 1 : size(X,2)
acfs = autocor(X[:,j], 1:mk)
durbin!(acfs, y, p)
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
Expand All @@ -590,7 +593,7 @@ using the Yule-Walker equations.

`r` must be a matrix of size `(length(lags), size(x, 2))`.
"""
function pacf!(r::AbstractMatrix{<:Real}, X::AbstractMatrix{T}, lags::AbstractVector{<:Integer}; method::Symbol=:regression) where T<:Union{Float32, Float64}
function pacf!(r::AbstractMatrix, X::AbstractMatrix, lags::AbstractVector{<:Integer}; method::Symbol=:regression)
lx = size(X, 1)
m = length(lags)
minlag, maxlag = extrema(lags)
Expand All @@ -611,7 +614,7 @@ end
"""
pacf(X, lags; method=:regression)

Compute the partial autocorrelation function (PACF) of a real-valued vector
Compute the partial autocorrelation function (PACF) of a vector
or matrix `X` at `lags`. `method` designates the estimation method. Recognized
values are `:regression`, which computes the partial autocorrelations via successive
regression models, and `:yulewalker`, which computes the partial autocorrelations
Expand All @@ -621,11 +624,11 @@ If `x` is a vector, return a vector of the same length as `lags`.
If `x` is a matrix, return a matrix of size `(length(lags), size(x, 2))`,
where each column in the result corresponds to a column in `x`.
"""
function pacf(X::AbstractMatrix{<:Real}, lags::AbstractVector{<:Integer}; method::Symbol=:regression)
function pacf(X::AbstractMatrix, lags::AbstractVector{<:Integer}; method::Symbol=:regression)
out = Matrix{float(eltype(X))}(undef, length(lags), size(X,2))
pacf!(out, float(X), lags; method=method)
end

function pacf(x::AbstractVector{<:Real}, lags::AbstractVector{<:Integer}; method::Symbol=:regression)
function pacf(x::AbstractVector, lags::AbstractVector{<:Integer}; method::Symbol=:regression)
vec(pacf(reshape(x, length(x), 1), lags, method=method))
end
Loading