From dd4d79af73c164524afa4c70ee634e117e0e53c7 Mon Sep 17 00:00:00 2001 From: dsweber2 Date: Tue, 30 Nov 2021 17:47:09 -0800 Subject: [PATCH] unused file and better notation --- src/apply.jl | 119 +++++++++++++++++++++--------------------- src/averagingTypes.jl | 0 2 files changed, 59 insertions(+), 60 deletions(-) delete mode 100644 src/averagingTypes.jl diff --git a/src/apply.jl b/src/apply.jl index f57e490..c942808 100644 --- a/src/apply.jl +++ b/src/apply.jl @@ -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) @@ -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 @@ -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 @@ -82,78 +82,78 @@ 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 @@ -161,19 +161,19 @@ 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 @@ -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 @@ -195,9 +195,9 @@ 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 @@ -205,10 +205,10 @@ 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 @@ -216,14 +216,13 @@ function reflect(Y, bt) 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 @@ -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 diff --git a/src/averagingTypes.jl b/src/averagingTypes.jl deleted file mode 100644 index e69de29..0000000