Skip to content

Commit

Permalink
unused file and better notation
Browse files Browse the repository at this point in the history
  • Loading branch information
dsweber2 committed Dec 1, 2021
1 parent 086686d commit dd4d79a
Show file tree
Hide file tree
Showing 2 changed files with 59 additions and 60 deletions.
119 changes: 59 additions & 60 deletions src/apply.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -19,11 +19,11 @@
assumption is that the sampling rate is 2kHz.
"""
function cwt(Y::AbstractArray{T,N}, cWav::CWT, daughters, fftPlans = 1) where {N, T}
@assert typeof(N)<:Integer
function cwt(Y::AbstractArray{T,N}, cWav::CWT, daughters, fftPlans = 1) where {N,T}
@assert typeof(N) <: Integer
# vectors behave a bit strangely, so we reshape them
if N==1
Y= reshape(Y,(length(Y), 1))
if N == 1
Y = reshape(Y, (length(Y), 1))
end
n1 = size(Y, 1)

Expand All @@ -35,8 +35,8 @@ function cwt(Y::AbstractArray{T,N}, cWav::CWT, daughters, fftPlans = 1) where {N
x̂, fromPlan = prepSignalAndPlans(x, cWav, fftPlans)
# If the vector isn't long enough to actually have any other scales, just
# return the averaging
if nScales <= 0 || size(daughters,2) == 1
daughters = daughters[:,1:1]
if nScales <= 0 || size(daughters, 2) == 1
daughters = daughters[:, 1:1]
nScales = 1
end

Expand All @@ -61,8 +61,8 @@ function cwt(Y::AbstractArray{T,N}, cWav::CWT, daughters, fftPlans = 1) where {N
wave = permutedims(wave, [1, ndims(wave), (2:(ndims(wave)-1))...])
ax = axes(wave)
wave = wave[1:n1, ax[2:end]...]
if N==1
wave = dropdims(wave, dims=3)
if N == 1
wave = dropdims(wave, dims = 3)
end

return wave
Expand All @@ -82,98 +82,98 @@ end
# yes | both | fft
# no | rfft | fft
# Analytic on Real input
function prepSignalAndPlans(x::AbstractArray{T}, cWav::CWT{W,S,WaTy, N, true}, fftPlans) where {T <: Real, W,S,WaTy,N}
function prepSignalAndPlans(x::AbstractArray{T}, cWav::CWT{W,S,WaTy,N,true}, fftPlans) where {T<:Real,W,S,WaTy,N}
# analytic wavelets that are being applied on real inputs
if fftPlans isa Tuple{<:AbstractFFTs.Plan{<:Real}, <:AbstractFFTs.Plan{<:Complex}}
if fftPlans isa Tuple{<:AbstractFFTs.Plan{<:Real},<:AbstractFFTs.Plan{<:Complex}}
# they handed us the right kind of thing, so no need to make new ones
x̂ = fftPlans[1] * x
fromPlan = fftPlans[2]
else
toPlan = plan_rfft(x,1)
toPlan = plan_rfft(x, 1)
x̂ = toPlan * x
fromPlan = plan_fft(x,1)
fromPlan = plan_fft(x, 1)
end
return x̂, fromPlan
end

# Non-analytic on Real input
function prepSignalAndPlans(x::AbstractArray{T}, cWav::CWT{W,S,WaTy, N, false}, fftPlans) where {T <: Real, W,S,WaTy,N}
function prepSignalAndPlans(x::AbstractArray{T}, cWav::CWT{W,S,WaTy,N,false}, fftPlans) where {T<:Real,W,S,WaTy,N}
# real wavelets that are being applied on real inputs
if fftPlans isa AbstractFFTs.Plan{<:Real}
# they handed us the right kind of thing, so no need to make new ones
x̂ = fftPlans * x
fromPlan = fftPlans
else
fromPlan = plan_rfft(x,1)
fromPlan = plan_rfft(x, 1)
x̂ = fromPlan * x
end
return x̂, fromPlan
end
# complex input
function prepSignalAndPlans(x::AbstractArray{T}, cWav, fftPlans) where {T <: Complex, W,S,WaTy,N}
function prepSignalAndPlans(x::AbstractArray{T}, cWav, fftPlans) where {T<:Complex}
# real wavelets that are being applied on real inputs
if fftPlans isa AbstractFFTs.Plan{<:Complex}
# they handed us the right kind of thing, so no need to make new ones
x̂ = fftPlans * x
fromPlan = fftPlans
else
fromPlan = plan_fft(x,1)
fromPlan = plan_fft(x, 1)
x̂ = fromPlan * x
end
return x̂, fromPlan
end

# analytic on real data with an averaging function
function analyticTransformReal!(wave, daughters, x̂, fftPlan, averagingType::Union{Father, Dirac})
function analyticTransformReal!(wave, daughters, x̂, fftPlan, ::Union{Father,Dirac})
outer = axes(x̂)[2:end]
n1 = size(x̂, 1)
isSourceEven = mod(size(wave,1)+1,2)
isSourceEven = mod(size(wave, 1) + 1, 2)
# the averaging function isn't analytic, so we need to do both positive and
# negative frequencies
@views tmpWave = x̂ .* daughters[:,1]
@views wave[(n1+1):end, outer..., 1] = reverse(conj.(tmpWave[2:end-isSourceEven, outer...]), dims=1)
@views tmpWave = x̂ .* daughters[:, 1]
@views wave[(n1+1):end, outer..., 1] = reverse(conj.(tmpWave[2:end-isSourceEven, outer...]), dims = 1)
@views wave[1:n1, outer..., 1] = tmpWave
@views wave[:, outer..., 1] = fftPlan \ (wave[:, outer..., 1]) # averaging
for j in 2:size(daughters,2)
@views wave[1:n1, outer..., j] = x̂ .* daughters[:,j]
for j = 2:size(daughters, 2)
@views wave[1:n1, outer..., j] = x̂ .* daughters[:, j]
wave[:, outer..., j] = fftPlan \ (wave[:, outer..., j]) # wavelet transform
end
end

# analytic on complex data with an averaging function
function analyticTransformComplex!(wave, daughters, x̂, fftPlan, averagingType::Union{Father, Dirac})
function analyticTransformComplex!(wave, daughters, x̂, fftPlan, ::Union{Father,Dirac})
outer = axes(x̂)[2:end]
n1 = size(daughters, 1)
isSourceEven = mod(size(wave,1)+1,2)
isSourceEven = mod(size(wave, 1) + 1, 2)
# the averaging function isn't analytic, so we need to do both positive and
# negative frequencies
@views positiveFreqs = x̂[1:n1, outer...] .* daughters[:,1]
@views negativeFreqs = x̂[(n1-isSourceEven+1):end, outer...] .* reverse(conj.(daughters[2:end,1]))
@views positiveFreqs = x̂[1:n1, outer...] .* daughters[:, 1]
@views negativeFreqs = x̂[(n1-isSourceEven+1):end, outer...] .* reverse(conj.(daughters[2:end, 1]))
@views wave[(n1-isSourceEven+1):end, outer..., 1] = negativeFreqs
@views wave[1:n1, outer..., 1] = positiveFreqs
@views wave[:, outer..., 1] = fftPlan \ (wave[:, outer..., 1]) # averaging
for j in 2:size(daughters,2)
@views wave[1:n1, outer..., j] = x̂[1:n1, outer...] .* daughters[:,j]
for j = 2:size(daughters, 2)
@views wave[1:n1, outer..., j] = x̂[1:n1, outer...] .* daughters[:, j]
@views wave[:, outer..., j] = fftPlan \ (wave[:, outer..., j]) # wavelet transform
end
end

function analyticTransformComplex!(wave, daughters, x̂, fftPlan, averagingType)
outer = axes(x̂)[2:end]
n1 = size(x̂, 1)
for j in 1:size(daughters,2)
@views wave[1:n1, outer..., j] = x̂[1:n1, outer...] .* daughters[:,j]
for j = 1:size(daughters, 2)
@views wave[1:n1, outer..., j] = x̂[1:n1, outer...] .* daughters[:, j]
@views wave[:, outer..., j] = fftPlan \ (wave[:, outer..., j]) # wavelet transform
end
end

# analytic on real data without an averaging function
function analyticTransformReal!(wave, daughters, x̂, fftPlan, averagingType::NoAve)
function analyticTransformReal!(wave, daughters, x̂, fftPlan, ::NoAve)
outer = axes(x̂)[2:end]
n1 = size(x̂, 1)
# the no averaging version
for j in 1:size(daughters,2)
wave[1:n1, outer..., j] = x̂ .* daughters[:,j]
for j = 1:size(daughters, 2)
wave[1:n1, outer..., j] = x̂ .* daughters[:, j]
wave[:, outer..., j] = fftPlan \ (wave[:, outer..., j]) # wavelet transform
end
end
Expand All @@ -183,8 +183,8 @@ function otherwiseTransform!(wave::AbstractArray{<:Real}, daughters, x̂, fromPl
# real wavelets on real data: that just makes sense
outer = axes(x̂)[2:end]
n1 = size(x̂, 1)
for j in 1:size(daughters, 2)
@views tmp = x̂ .* daughters[:,j]
for j = 1:size(daughters, 2)
@views tmp = x̂ .* daughters[:, j]
@views wave[:, outer..., j] = fromPlan \ tmp # wavelet transform
end
end
Expand All @@ -195,35 +195,34 @@ function otherwiseTransform!(wave::AbstractArray{<:Complex}, daughters, x̂, fro
outer = axes(x̂)[2:end]
n1 = size(daughters, 1)
isSourceEven = mod(size(fromPlan, 1) + 1, 2)
for j in 1:size(daughters, 2)
@views wave[1:n1, outer..., j] = @views x̂[1:n1,outer...] .* daughters[:,j]
@views wave[n1-isSourceEven+1:end, outer..., j] = x̂[n1-isSourceEven+1:end,outer...] .* reverse(conj.(daughters[2:end,j]))
for j = 1:size(daughters, 2)
@views wave[1:n1, outer..., j] = @views x̂[1:n1, outer...] .* daughters[:, j]
@views wave[n1-isSourceEven+1:end, outer..., j] = x̂[n1-isSourceEven+1:end, outer...] .* reverse(conj.(daughters[2:end, j]))
@views wave[:, outer..., j] = fromPlan \ (wave[:, outer..., j]) # wavelet transform
end
end

function reflect(Y, bt)
n1 = size(Y, 1)
if typeof(bt) <: ZPBoundary
base2 = ceil(Int, log2(n1)); # power of 2 nearest to N
x = cat(Y, zeros(2^(base2)-n1, size(Y)[2:end]...), dims=1)
base2 = ceil(Int, log2(n1)) # power of 2 nearest to N
x = cat(Y, zeros(2^(base2) - n1, size(Y)[2:end]...), dims = 1)
elseif typeof(bt) <: SymBoundary
x = cat(Y, reverse(Y,dims = 1), dims = 1)
x = cat(Y, reverse(Y, dims = 1), dims = 1)
else
x = Y
end
return x
end


function cwt(Y::AbstractArray{T}, c::CWT{W}; varArgs...) where {T<:Number, S<:Real, V<: Real,
W<:WaveletBoundary}
daughters,ω = computeWavelets(size(Y, 1), c; varArgs...)
function cwt(Y::AbstractArray{T}, c::CWT{W}; varArgs...) where {T<:Number,W<:WaveletBoundary}
daughters, ω = computeWavelets(size(Y, 1), c; varArgs...)
return cwt(Y, c, daughters)
end

cwt(Y::AbstractArray{T}, w::ContWaveClass; varargs...) where {T<:Number, S<:Real, V<:Real} = cwt(Y,CWT(w); varargs...)
cwt(Y::AbstractArray{T}) where T<:Real = cwt(Y,Morlet())
cwt(Y::AbstractArray{T}, w::ContWaveClass; varargs...) where {T<:Number} = cwt(Y, CWT(w); varargs...)
cwt(Y::AbstractArray{T}) where {T<:Real} = cwt(Y, Morlet())

abstract type InverseType end
struct DualFrames <: InverseType end
Expand All @@ -244,35 +243,35 @@ Return the inverse continuous wavelet transform, computed using the simple dual
Return the inverse continuous wavelet transform, computed using the canonical dual frame ``\\tilde{\\widehat{ψ}} = \\frac{ψ̂_n(ω)}{∑_n\\|ψ̂_n(ω)\\|^2}``. The algorithm is to compute the cwt again, but using the canonical dual frame; consiquentially, it is the most computationally intensive of the three algorithms, and typically the best behaved. Will be numerically unstable if the high frequencies of all of the wavelets are too small however, and tends to fail spectacularly in this case.
"""
function icwt(res::AbstractArray, cWav::CWT, inverseStyle::PenroseDelta)
Ŵ = computeWavelets(size(res,1), cWav)[1]
function icwt(res::AbstractArray, cWav::CWT, ::PenroseDelta)
Ŵ = computeWavelets(size(res, 1), cWav)[1]
β = computeDualWeights(Ŵ, cWav)
testDualCoverage(β, Ŵ)
compXRecon = sum(res .* β, dims=2)
imagXRecon = irfft(im*rfft(imag.(compXRecon),1), size(compXRecon,1)) # turns out the dual frame for the imaginary part is rather gross in the time domain
compXRecon = sum(res .* β, dims = 2)
imagXRecon = irfft(im * rfft(imag.(compXRecon), 1), size(compXRecon, 1)) # turns out the dual frame for the imaginary part is rather gross in the time domain
return imagXRecon + real.(compXRecon)
end

function icwt(res::AbstractArray, cWav::CWT, inverseStyle::NaiveDelta)
Ŵ = computeWavelets(size(res,1), cWav)[1]
β = computeNaiveDualWeights(Ŵ, cWav, size(res,1))
function icwt(res::AbstractArray, cWav::CWT, ::NaiveDelta)
Ŵ = computeWavelets(size(res, 1), cWav)[1]
β = computeNaiveDualWeights(Ŵ, cWav, size(res, 1))
testDualCoverage(β, Ŵ)
compXRecon = sum(res .* β, dims=2)
imagXRecon = irfft(im*rfft(imag.(compXRecon),1), size(compXRecon,1)) # turns out the dual frame for the imaginary part is rather gross in the time domain
compXRecon = sum(res .* β, dims = 2)
imagXRecon = irfft(im * rfft(imag.(compXRecon), 1), size(compXRecon, 1)) # turns out the dual frame for the imaginary part is rather gross in the time domain
return imagXRecon + real.(compXRecon)
end

function icwt(res::AbstractArray, cWav::CWT, inverseStyle::DualFrames)
Ŵ = computeWavelets(size(res,1), cWav)[1]
function icwt(res::AbstractArray, cWav::CWT, ::DualFrames)
Ŵ = computeWavelets(size(res, 1), cWav)[1]
canonDualFrames = explicitConstruction(Ŵ)
testDualCoverage(canonDualFrames)
convolved = cwt(res, cWav, canonDualFrames)
ax = axes(convolved)
@views xRecon = sum(convolved[:,i,i,ax[4:end]...] for i=1:size(Ŵ,2))
@views xRecon = sum(convolved[:, i, i, ax[4:end]...] for i = 1:size(Ŵ, 2))
return xRecon
end

icwt(Y::AbstractArray, w::ContWaveClass; varargs...) where {S<:Real, T<:Real, V<:Real} = icwt(Y,CWT(w), PenroseDelta(); varargs...)
icwt(Y::AbstractArray; varargs...) = icwt(Y,Morlet(), PenroseDelta(); varargs...)
icwt(Y::AbstractArray, w::ContWaveClass; varargs...) = icwt(Y, CWT(w), PenroseDelta(); varargs...)
icwt(Y::AbstractArray; varargs...) = icwt(Y, Morlet(), PenroseDelta(); varargs...)

# CWT (continuous wavelet transform directly) TODO: direct if sufficiently small
Empty file removed src/averagingTypes.jl
View file
Empty file.

0 comments on commit dd4d79a

Please sign in to comment.