From 0bd8baa9950a05ed3e6e88274805347a036eb7fc Mon Sep 17 00:00:00 2001 From: nHackel Date: Fri, 9 Aug 2024 21:28:56 +0200 Subject: [PATCH 1/7] Init handsfree kaczmarz port --- src/Algorithms/Algorithms.jl | 3 +- .../HandsFreeAlgorithms.jl | 1 + .../HandsFreeAlgorithms/HandsFreeKaczmarz.jl | 206 ++++++++++++++++++ 3 files changed, 209 insertions(+), 1 deletion(-) create mode 100644 src/Algorithms/HandsFreeAlgorithms/HandsFreeAlgorithms.jl create mode 100644 src/Algorithms/HandsFreeAlgorithms/HandsFreeKaczmarz.jl diff --git a/src/Algorithms/Algorithms.jl b/src/Algorithms/Algorithms.jl index 3c05491..eb768a7 100644 --- a/src/Algorithms/Algorithms.jl +++ b/src/Algorithms/Algorithms.jl @@ -1,2 +1,3 @@ include("SinglePatchAlgorithms/SinglePatchAlgorithms.jl") -include("MultiPatchAlgorithms/MultiPatchAlgorithms.jl") \ No newline at end of file +include("MultiPatchAlgorithms/MultiPatchAlgorithms.jl") +include("HandsFreeAlgorithms/HandsFreeAlgorithms.jl") \ No newline at end of file diff --git a/src/Algorithms/HandsFreeAlgorithms/HandsFreeAlgorithms.jl b/src/Algorithms/HandsFreeAlgorithms/HandsFreeAlgorithms.jl new file mode 100644 index 0000000..436e7be --- /dev/null +++ b/src/Algorithms/HandsFreeAlgorithms/HandsFreeAlgorithms.jl @@ -0,0 +1 @@ +include("HandsFreeKaczmarz.jl") \ No newline at end of file diff --git a/src/Algorithms/HandsFreeAlgorithms/HandsFreeKaczmarz.jl b/src/Algorithms/HandsFreeAlgorithms/HandsFreeKaczmarz.jl new file mode 100644 index 0000000..02577bf --- /dev/null +++ b/src/Algorithms/HandsFreeAlgorithms/HandsFreeKaczmarz.jl @@ -0,0 +1,206 @@ +export HandsFreeKaczmarzState +mutable struct HandsFreeKaczmarzState{T, Tc <: Union{T, Complex{T}}, vecTc <: AbstractArray{Tc}, vecT <: AbstractArray{T}, R <: AbstractRegularization} <: RegularizedLeastSquares.AbstractSolverState{Kaczmarz} + # Usual Kaczmarz + u::vecTc + x::vecTc + vl::vecTc + regs::Vector{R} + εw::vecT + τl::Tc + αl::Tc + iteration::Int64 + # HandsFree + freqSNRs::vecT + curvatures::vecT + norm_cl::vecT + resid::vecT + iterbounds::Tuple{Int64,Int64} + stoppingParas::Tuple{Float64,Float64} + expected_iters::Int64 + wanted_iters::Int64 + SNRbounds::Tuple{Float64,Float64} +end + +function RegularizedLeastSquares.Kaczmarz(A, freqSNRs + ; reg = L2Regularization(real(eltype(A))(5.0)) # encodes start value for λ + , normalizeReg::AbstractRegularizationNormalization = NoNormalization() + , iterbounds::Tuple{Int64,Int64}=(1,50) + , stoppingParas::Tuple{Float64,Float64}=(0.25,2.0) + , expected_iters::Int64=0 + , wanted_iters::Int64=0 + , SNRbounds::Tuple{Float64,Float64}=(60.0,1.5) + ) + + T = real(eltype(A)) + + min_iter, max_iter = iterbounds + + # Prepare regularization terms + reg = isa(reg, AbstractVector) ? reg : [reg] + idx = findsink(L2Regularization, reg) + regs = AbstractRegularization[] + if length(idx) == max_iter # User provided enough L2Regularization terms + regs = reg[idx] + elseif length(idx) == 1 + startλ = λ(reg[idx]) + regs = [L2Regularization(startλ / (1 + (0.2 * i - 0.2)^5)) for i in 1:max_iter] + else + error("HandsFreeKaczmarz requires either one or as many as the upper itterbound L2Regularization terms, found $(length(idx))") + end + deleteat!(reg, idx) + regs = RegularizedLeastSquares.normalize(Kaczmarz, normalizeReg, regs, A, nothing) + + # Tikhonov matrix is only valid with NoNormalization or SystemMatrixBasedNormalization + if λ(first(regs)) isa AbstractVector + error("Tikhonov matrix for Handsfree Kaczmarz is not yet implemented") + end + + indices = findsinks(AbstractProjectionRegularization, reg) + other = AbstractRegularization[reg[i] for i in indices] + deleteat!(reg, indices) + if length(reg) == 1 + push!(other, reg[1]) + elseif length(reg) > 1 + error("Kaczmarz does not allow for more than one additional regularization term, found $(length(reg))") + end + other = identity.(other) + + # setup denom and rowindex + denom, rowindex = inithandsfree(A) + rowIndexCycle = collect(1:length(rowindex)) + + M,N = size(A) + + u = zeros(eltype(A),M) + x = zeros(eltype(A),N) + vl = zeros(eltype(A),M) + εw = zeros(T, max_iter) + τl = zero(eltype(A)) + αl = zero(eltype(A)) + + state = HandsFreeKaczmarzState(u, x, vl, regs, εw, τl, αl, 0, freqSNRs, zeros(T, max_iter), zeros(T, max_iter), zeros(T, max_iter), iterbounds, stoppingParas, expected_iters, wanted_iters, SNRbounds) + + return Kaczmarz(A, nothing, other, denom, rowindex, rowIndexCycle, + false, 0, T[], false, + Int64(0), normalizeReg, max_iter, state) +end + +function inithandsfree(A) + T = real(eltype(A)) + denom = T[] + rowindex = Int64[] + + for i = 1:size(A, 1) + s² = rownorm²(A,i) + if s²>0 + push!(denom,s²) # only compute rownorm2, λ is computed during iterations + push!(rowindex,i) + end + end + return denom, rowindex +end + +function RegularizedLeastSquares.init!(solver::Kaczmarz, state::HandsFreeKaczmarzState{T, Tc, vecTc, vecT}, b::vecTc; x0 = 0) where {T, Tc, vecTc, vecT} + solver.reg = RegularizedLeastSquares.normalize(solver, solver.normalizeReg, solver.reg, solver.A, b) + state.regs = RegularizedLeastSquares.normalize(solver, solver.normalizeReg, state.regs, solver.A, b) + + #if solver.shuffleRows || solver.randomized + # Random.seed!(solver.seed) + #end + #if solver.shuffleRows + # shuffle!(solver.rowIndexCycle) + #end + + # start vector + state.x .= x0 + state.vl .= 0 + + state.resid[:] .= zero(T) + state.norm_cl[:] .= zero(T) + state.curvatures[:] .= zero(T) + + + state.u .= b + state.εw .= sqrt.(λ.(state.regs)) + state.iteration = 1 +end + +function iterate(solver::Kaczmarz{matT, Nothing}, state::HandsFreeKaczmarzState{T, Tc, vecTc, vecT} = solver.state) where {matT, T, Tc, vecTc, vecT} + if RegularizedLeastSquares.done(solver,state) return nothing end + iteration = state.iteration + + usedIndices = filterFrequencies(solver, state) + + for (i, row) in enumerate(usedIndices) + RegularizedLeastSquares.iterate_row_index(solver, state, solver.A, row, i, iteration) + end + + for r in solver.reg + prox!(r, state.x) + end + + state.norm_cl[iteration] = norm(real(state.x)) + state.resid[iteration] = norm(- sqrt(λ(state.regs[iteration])) * state.vl[usedIndices]) / length(usedIndices) + + if iteration == 1 + state.curvatures[iteration] = 0 + elseif iteration == 2 + dcdr = (state.norm_cl[iteration]-state.norm_cl[iteration-1]) / (state.resid[iteration]-state.resid[iteration-1]) + state.curvatures[iteration] = ( (dcdr - 0) / (state.resid[iteration]-state.resid[iteration-1]) ) / ((1 + dcdr^2)^(3/2)) + else + dcdr_old = (state.norm_cl[iteration-1] - state.norm_cl[iteration-2]) / (state.resid[iteration-1] - state.resid[iteration-2]) + dcdr = (state.norm_cl[iteration] - state.norm_cl[iteration-1]) / (state.resid[iteration] - state.resid[iteration-1]) + state.curvatures[iteration] = ( (dcdr - dcdr_old) / (state.resid[iteration]-state.resid[iteration-1]) ) / ((1 + dcdr^2)^(3/2)) + end + + state.iteration += 1 + return state.x, state +end + +function RegularizedLeastSquares.iterate_row_index(solver::Kaczmarz, state::HandsFreeKaczmarzState, A, row, index, iteration) + state.τl = dot_with_matrix_row(A,state.x,row) + state.αl = 1/(solver.denom[index] + state.ɛw[iteration]^2) * (state.u[row]-state.τl-state.ɛw[iteration]*state.vl[row]) + kaczmarz_update!(A,state.x,row,state.αl) + state.vl[row] += state.αl*state.ɛw[iteration] +end + + +function filterFrequencies(::Kaczmarz, state::HF) where {T, Tc, HF <: HandsFreeKaczmarzState{T, Tc}} + iteration = state.iteration + threshl = [j > state.SNRbounds[2] ? T(j) : T(state.SNRbounds[2]) for j in (state.SNRbounds[1]) / (1 + (0.28 * iteration - 0.28)^2)][1] + indexl = sort(findall(x -> x > threshl, state.freqSNRs)) # potentially drop sorting + return indexl +end + +function RegularizedLeastSquares.done(::Kaczmarz,state::HF) where HF <: HandsFreeKaczmarzState + iteration = state.iteration + if state.expected_iters != 0 && iteration >= round(Int,state.expected_iters*3/2) + tmp=state.expected_iters + @info "Would like to go more than $iteration iterations, but expect $tmp iterations." + if !(state.wanted_iters < state.expected_iters) + state.wanted_iters = iteration+2 + end + return true + elseif iteration >= 2 && iteration >= state.iterbounds[1] && state.curvatures[iteration] > state.stoppingParas[1] * state.norm_cl[1] && abs(state.curvatures[iteration-1]*state.stoppingParas[2]) < abs(state.curvatures[iteration]) + if iteration >= round(Int,state.expected_iters*1/2) + tmp=state.expected_iters + @info "Stopped after $iteration iterations. Expected $tmp iterations." + state.wanted_iters = iteration + return true + else + tmp = state.expected_iters + tmp2 = round(Int,state.expected_iters*2/5) + state.expected_iters = tmp2 + state.wanted_iters = iteration + @info "Would like to stop after $iteration iterations, but expect $tmp iterations. Update expected iterations for this reco to $tmp2." + return false + end + elseif iteration >= state.iterbounds[2] + tmp = state.expected_iters + @info "Stopped at the max iter-bound of $iteration iterations. Expected $tmp iterations." + state.wanted_iters = iteration + return true + else + return false + end +end \ No newline at end of file From 0196a8031f5f325054fe0793b572de568c2fb6ad Mon Sep 17 00:00:00 2001 From: nHackel Date: Sat, 10 Aug 2024 18:00:18 +0200 Subject: [PATCH 2/7] Fix bug in hands free row selection --- src/Algorithms/HandsFreeAlgorithms/HandsFreeKaczmarz.jl | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/Algorithms/HandsFreeAlgorithms/HandsFreeKaczmarz.jl b/src/Algorithms/HandsFreeAlgorithms/HandsFreeKaczmarz.jl index 02577bf..eb70a19 100644 --- a/src/Algorithms/HandsFreeAlgorithms/HandsFreeKaczmarz.jl +++ b/src/Algorithms/HandsFreeAlgorithms/HandsFreeKaczmarz.jl @@ -131,7 +131,8 @@ function iterate(solver::Kaczmarz{matT, Nothing}, state::HandsFreeKaczmarzState{ usedIndices = filterFrequencies(solver, state) - for (i, row) in enumerate(usedIndices) + for i in usedIndices + row = solver.rowindex[i] RegularizedLeastSquares.iterate_row_index(solver, state, solver.A, row, i, iteration) end @@ -168,7 +169,7 @@ end function filterFrequencies(::Kaczmarz, state::HF) where {T, Tc, HF <: HandsFreeKaczmarzState{T, Tc}} iteration = state.iteration threshl = [j > state.SNRbounds[2] ? T(j) : T(state.SNRbounds[2]) for j in (state.SNRbounds[1]) / (1 + (0.28 * iteration - 0.28)^2)][1] - indexl = sort(findall(x -> x > threshl, state.freqSNRs)) # potentially drop sorting + indexl = sort(findall(x -> x > threshl, state.freqSNRs)) return indexl end From 8cde5bbb0e17422865ac9a2ef1a57a0977b5d950 Mon Sep 17 00:00:00 2001 From: nHackel Date: Sat, 10 Aug 2024 19:18:40 +0200 Subject: [PATCH 3/7] Working HandFreeKaczmarz port --- src/Algorithms/HandsFreeAlgorithms/HandsFreeKaczmarz.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Algorithms/HandsFreeAlgorithms/HandsFreeKaczmarz.jl b/src/Algorithms/HandsFreeAlgorithms/HandsFreeKaczmarz.jl index eb70a19..5fd0966 100644 --- a/src/Algorithms/HandsFreeAlgorithms/HandsFreeKaczmarz.jl +++ b/src/Algorithms/HandsFreeAlgorithms/HandsFreeKaczmarz.jl @@ -126,7 +126,6 @@ function RegularizedLeastSquares.init!(solver::Kaczmarz, state::HandsFreeKaczmar end function iterate(solver::Kaczmarz{matT, Nothing}, state::HandsFreeKaczmarzState{T, Tc, vecTc, vecT} = solver.state) where {matT, T, Tc, vecTc, vecT} - if RegularizedLeastSquares.done(solver,state) return nothing end iteration = state.iteration usedIndices = filterFrequencies(solver, state) @@ -154,6 +153,7 @@ function iterate(solver::Kaczmarz{matT, Nothing}, state::HandsFreeKaczmarzState{ state.curvatures[iteration] = ( (dcdr - dcdr_old) / (state.resid[iteration]-state.resid[iteration-1]) ) / ((1 + dcdr^2)^(3/2)) end + if RegularizedLeastSquares.done(solver,state) return nothing end state.iteration += 1 return state.x, state end From 3cb3402f9a4cad73358d0ad2de87ef8da6dc256c Mon Sep 17 00:00:00 2001 From: nHackel Date: Mon, 19 Aug 2024 13:07:51 +0200 Subject: [PATCH 4/7] Add getNoiseLevel from HandsFreeReco Co-authored-by: Konrad Scheffler <83588313+KonradScheffler@users.noreply.github.com> --- src/PreProcessing/FrequencyFilterParameter.jl | 19 +++++++++++++ src/SystemMatrix/SystemMatrix.jl | 27 +++++++++++++++++++ 2 files changed, 46 insertions(+) diff --git a/src/PreProcessing/FrequencyFilterParameter.jl b/src/PreProcessing/FrequencyFilterParameter.jl index a164b04..3b73b10 100644 --- a/src/PreProcessing/FrequencyFilterParameter.jl +++ b/src/PreProcessing/FrequencyFilterParameter.jl @@ -85,4 +85,23 @@ Base.@kwdef struct CompositeFrequencyFilterParameters <: AbstractFrequencyFilter end function process(algoT::Type{<:AbstractMPIRecoAlgorithm}, params::CompositeFrequencyFilterParameters, file::MPIFile) return reduce(intersect, filter(!isnothing, map(p -> process(algoT, p, file), params.filters))) +end + +export NoiseLevelFrequencyFilterParameter +Base.@kwdef struct NoiseLevelFrequencyFilterParameter <: AbstractFrequencyFilterParameter + noiseMeas::MPIFile + frames::Union{Vector{Int64}, UnitRange{Int64}} = measBGFrameIdx(noiseMeas) + minFreq::Float64 = 0.0 + maxFreq::Union{Float64, Nothing} = nothing + recChannels::Union{UnitRange{Int64}, Nothing} = nothing + SNRThresh::Float64=-1.0 + numPeriodAverages::Int64 = 1 + numPeriodGrouping::Int64 = 1 + maxMixingOrder::Int64 = -1 + numSidebandFreqs::Int64 = -1 +end +function process(::Type{<:AbstractMPIRecoAlgorithm})(params::NoiseLevelFrequencyFilterParameter, file::MPIFile) + noiseLevel = getNoiseLevel(params.noiseMeas, params.frames, params.recChannels) + kwargs = toKwargs(params, ignore = [:noiseMeas, :frames], default = Dict{Symbol, Any}(:maxFreq => rxBandwidth(file), :recChannels => 1:rxNumChannels(file))) + return filterFrequencies(file; kwargs..., SNRThresh = params.SNRThresh * noiseLevel) end \ No newline at end of file diff --git a/src/SystemMatrix/SystemMatrix.jl b/src/SystemMatrix/SystemMatrix.jl index 111935a..ab567a3 100644 --- a/src/SystemMatrix/SystemMatrix.jl +++ b/src/SystemMatrix/SystemMatrix.jl @@ -306,3 +306,30 @@ function getSF(bSF::MPIFile, mx::Int, my::Int, mz::Int, recChan::Int; bgcorrecti end include("Sparse.jl") + +""" +Calculates a noise level from empty measurement `bEmpty`. + + NoiseLevel = getNoiseLevel(bEmpty,bgframes,channels) +""" +function getNoiseLevel(bEmpty, bgframes, channels) + tmp = getMeasurementsFD(bEmpty, true; frames=bgframes)[:, channels, :, :] + if length(channels) > 1 + numfreq, numchannels, numpatches, numframes = size(tmp) + measBG = reshape(permutedims(tmp, [1, 3, 2, 4]), (numfreq * numpatches, numchannels, numframes)) + else + numchannels = 1 + numfreq, numpatches, numframes = size(tmp) + measBG = reshape(tmp, (numfreq * numpatches, numchannels, numframes)) + end + noise = zeros(numfreq * numpatches, numchannels) + for r = 1:numchannels + for k = 1:numfreq*numpatches + tmp = view(measBG, k, r, :) + #maxBG = mapreduce(abs, max, tmp) # maximum(abs.(measBG[k, r, :])) + meanBG = mean(tmp) + noise[k, r] = mean(abs.((tmp .- meanBG)))#./maxBG)) + end + end + return mean(noise) +end From c9a872f3672c4623fdb2181f0a2a2c4d21441b99 Mon Sep 17 00:00:00 2001 From: nHackel Date: Mon, 19 Aug 2024 16:15:44 +0200 Subject: [PATCH 5/7] Implement high-level handsfree single patch reco --- .../HandsFreeAlgorithms.jl | 4 +- .../HandsFreeLeastSquares.jl | 81 +++++++++++++++++++ .../HandsFreeSinglePatch.jl | 43 ++++++++++ src/PreProcessing/FrequencyFilterParameter.jl | 11 +-- 4 files changed, 133 insertions(+), 6 deletions(-) create mode 100644 src/Algorithms/HandsFreeAlgorithms/HandsFreeLeastSquares.jl create mode 100644 src/Algorithms/HandsFreeAlgorithms/HandsFreeSinglePatch.jl diff --git a/src/Algorithms/HandsFreeAlgorithms/HandsFreeAlgorithms.jl b/src/Algorithms/HandsFreeAlgorithms/HandsFreeAlgorithms.jl index 436e7be..77846c9 100644 --- a/src/Algorithms/HandsFreeAlgorithms/HandsFreeAlgorithms.jl +++ b/src/Algorithms/HandsFreeAlgorithms/HandsFreeAlgorithms.jl @@ -1 +1,3 @@ -include("HandsFreeKaczmarz.jl") \ No newline at end of file +include("HandsFreeKaczmarz.jl") +include("HandsFreeLeastSquares.jl") +include("HandsFreeSinglePatch.jl") \ No newline at end of file diff --git a/src/Algorithms/HandsFreeAlgorithms/HandsFreeLeastSquares.jl b/src/Algorithms/HandsFreeAlgorithms/HandsFreeLeastSquares.jl new file mode 100644 index 0000000..f18d9c5 --- /dev/null +++ b/src/Algorithms/HandsFreeAlgorithms/HandsFreeLeastSquares.jl @@ -0,0 +1,81 @@ +export HandsFreeSolverParameters +Base.@kwdef struct HandsFreeSolverParameters <: AbstractSolverParameters{Kaczmarz} + iterbounds::Tuple{Int64, Int64}=(1, 25) + enforceReal::Bool=true + enforcePositive::Bool=true + normalizeReg::AbstractRegularizationNormalization = SystemMatrixBasedNormalization() + startλ::Float64=5.0 + SNRbounds::Tuple{Float64, Float64}=(0, 1) + stoppingParas::Tuple{Float64, Float64}=(0.25, 2.0) + equalizeIters::Bool=false + flattenIters::Bool=false +end + +function process(t::Type{<:AbstractMPIRecoAlgorithm}, params::LeastSquaresParameters{Kaczmarz, O, SF, R, SL, W}, u::AbstractArray, snr::AbstractVector) where {O, SF, R, SL <: HandsFreeSolverParameters, W} + + solverParams = params.solverParams + N = size(params.S, 2) + M = div(length(params.S), N) + L = size(u)[end] + u = reshape(u, M, L) + c = zeros(N, L) + + reg, _ = prepareRegularization([L2Regularization(real(eltype(u))(solverParams.startλ))], params) + + S = params.S + if !isnothing(params.weights) + S = ProdOp(WeightingOp(params.weights), S) + for l = 1:L + u[:, l] = params.weights.*u[:, l] + end + end + + it = zeros(L) + w_it = zeros(L) + curv = [zeros(solverParams.iterbounds[2]) for i=1:L] + + solv = Kaczmarz(S, snr; reg = reg, iterbounds = solverParams.iterbounds, stoppingParas = solverParams.stoppingParas, SNRbounds = solverParams.SNRbounds, normalizeReg = solverParams.normalizeReg) + + for l=1:L + d = solve!(solv, u[:, l]) + + # Update curve meta data + curv[l] = solv.state.curvatures + it[l] = solv.state.iteration + w_it[l] = solv.state.wanted_iters + + # Store frame + if !isnothing(params.op) + d[:] = params.op*d + end + c[:, l] = Array(real(d)) + + if solverParams.equalizeIters && l >= 3 + solv.state.expected_iters = round(Int, sum(it[l-2:l].*[0.1, 0.3, 0.6])) + end + end + + if solverParams.flattenIters + for l=1:L + iter = round(Int, mean(it[maximum((1, l-5)):minimum((l+5, L))])) + solv.state.iterbounds = (iter, iter) + solv.state.expected_iters = iter + + d = solve!(solv, u[:, l]) + + # Update curve meta data + curv[l] = solv.state.curvatures + it[l] = solv.state.iteration + w_it[l] = solv.state.iteration + + # Store frame + if !isnothing(params.op) + d[:] = params.op*d + end + c[:, l] = Array(real(d)) + end + end + + return c + +end \ No newline at end of file diff --git a/src/Algorithms/HandsFreeAlgorithms/HandsFreeSinglePatch.jl b/src/Algorithms/HandsFreeAlgorithms/HandsFreeSinglePatch.jl new file mode 100644 index 0000000..e0bd076 --- /dev/null +++ b/src/Algorithms/HandsFreeAlgorithms/HandsFreeSinglePatch.jl @@ -0,0 +1,43 @@ +export SinglePatchHandsFreeReconstructionParameter +Base.@kwdef struct SinglePatchHandsFreeReconstructionParameter{L<:AbstractSystemMatrixLoadingParameter, + arrT <: AbstractArray, SP<:HandsFreeSolverParameters, W<:AbstractWeightingParameters} <: AbstractSinglePatchReconstructionParameters + # File + sf::MPIFile + sfLoad::Union{L, ProcessResultCache{L}} + arrayType::Type{arrT} = Array + # Solver + solverParams::SP + #reg::Vector{R} = AbstractRegularization[] + weightingParams::Union{W, ProcessResultCache{W}} = NoWeightingParameters() +end + +function prepareSystemMatrix(reco::SinglePatchHandsFreeReconstructionParameter{L}) where {L<:AbstractSystemMatrixLoadingParameter} + freqs, sf, grid = process(AbstractMPIRecoAlgorithm, reco.sfLoad, reco.sf, Kaczmarz, reco.arrayType) + return freqs, sf, grid, reco.arrayType +end + +function prepareWeights(reco::SinglePatchHandsFreeReconstructionParameter{L,arrT,SP,W}, freqs, sf) where {L, arrT, SP, W<:AbstractWeightingParameters} + return process(AbstractMPIRecoAlgorithm, reco.weightingParams, freqs, sf, nothing, reco.arrayType) +end + +function process(algo::SinglePatchReconstructionAlgorithm, params::SinglePatchHandsFreeReconstructionParameter, u) + weights = process(algo, params.weightingParams, u, WeightingType(params.weightingParams)) + + B = getLinearOperator(algo, params) + + solver = LeastSquaresParameters(op = B, S = algo.S, reg = L2Regularization[], solverParams = params.solverParams, weights = weights) + + snr = real(eltype(algo.S)).(vec(MPIFiles.getCalibSNR(algo.sf)[algo.freqs, 1])) + + result = process(algo, solver, u, snr) + + return gridresult(result, algo.grid, algo.sf) +end + +function getLinearOperator(algo::SinglePatchReconstructionAlgorithm, params::SinglePatchHandsFreeReconstructionParameter{<:DenseSystemMatixLoadingParameter, S}) where {S} + return nothing +end + +function getLinearOperator(algo::SinglePatchReconstructionAlgorithm, params::SinglePatchHandsFreeReconstructionParameter{<:SparseSystemMatrixLoadingParameter, S}) where {S} + return process(algo, params.sfLoad, eltype(algo.S), algo.arrayType, tuple(shape(algo.grid)...)) +end \ No newline at end of file diff --git a/src/PreProcessing/FrequencyFilterParameter.jl b/src/PreProcessing/FrequencyFilterParameter.jl index 3b73b10..cda6347 100644 --- a/src/PreProcessing/FrequencyFilterParameter.jl +++ b/src/PreProcessing/FrequencyFilterParameter.jl @@ -87,21 +87,22 @@ function process(algoT::Type{<:AbstractMPIRecoAlgorithm}, params::CompositeFrequ return reduce(intersect, filter(!isnothing, map(p -> process(algoT, p, file), params.filters))) end +#= export NoiseLevelFrequencyFilterParameter Base.@kwdef struct NoiseLevelFrequencyFilterParameter <: AbstractFrequencyFilterParameter noiseMeas::MPIFile - frames::Union{Vector{Int64}, UnitRange{Int64}} = measBGFrameIdx(noiseMeas) + levelFactor::Float64 + noiseFrames::Union{Vector{Int64}, UnitRange{Int64}} = measBGFrameIdx(noiseMeas) minFreq::Float64 = 0.0 maxFreq::Union{Float64, Nothing} = nothing recChannels::Union{UnitRange{Int64}, Nothing} = nothing - SNRThresh::Float64=-1.0 numPeriodAverages::Int64 = 1 numPeriodGrouping::Int64 = 1 maxMixingOrder::Int64 = -1 numSidebandFreqs::Int64 = -1 end -function process(::Type{<:AbstractMPIRecoAlgorithm})(params::NoiseLevelFrequencyFilterParameter, file::MPIFile) +function process(::Type{<:AbstractMPIRecoAlgorithm}, params::NoiseLevelFrequencyFilterParameter, file::MPIFile) noiseLevel = getNoiseLevel(params.noiseMeas, params.frames, params.recChannels) kwargs = toKwargs(params, ignore = [:noiseMeas, :frames], default = Dict{Symbol, Any}(:maxFreq => rxBandwidth(file), :recChannels => 1:rxNumChannels(file))) - return filterFrequencies(file; kwargs..., SNRThresh = params.SNRThresh * noiseLevel) -end \ No newline at end of file + return filterFrequencies(file; kwargs..., SNRThresh = params.levelFactor * noiseLevel) +end=# \ No newline at end of file From ead4b6b2b1f881e962dac89d8d303b35c100b3f3 Mon Sep 17 00:00:00 2001 From: nHackel Date: Fri, 23 Aug 2024 11:33:10 +0200 Subject: [PATCH 6/7] Commented out getNoiseLevel (for now) --- src/SystemMatrix/SystemMatrix.jl | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/SystemMatrix/SystemMatrix.jl b/src/SystemMatrix/SystemMatrix.jl index ab567a3..ccd61f3 100644 --- a/src/SystemMatrix/SystemMatrix.jl +++ b/src/SystemMatrix/SystemMatrix.jl @@ -307,6 +307,7 @@ end include("Sparse.jl") +#= """ Calculates a noise level from empty measurement `bEmpty`. @@ -333,3 +334,4 @@ function getNoiseLevel(bEmpty, bgframes, channels) end return mean(noise) end +=# \ No newline at end of file From c345c04ea66044af8bad620cfb5808ef60ea84f6 Mon Sep 17 00:00:00 2001 From: nHackel Date: Fri, 23 Aug 2024 11:36:33 +0200 Subject: [PATCH 7/7] Add SinglePatchHandsFree plan --- config/SinglePatchHandsFree.toml | 98 ++++++++++++++++++++++++++++++++ 1 file changed, 98 insertions(+) create mode 100644 config/SinglePatchHandsFree.toml diff --git a/config/SinglePatchHandsFree.toml b/config/SinglePatchHandsFree.toml new file mode 100644 index 0000000..f595c4a --- /dev/null +++ b/config/SinglePatchHandsFree.toml @@ -0,0 +1,98 @@ +_module = "MPIReco" +_type = "RecoPlan{SinglePatchReconstructionAlgorithm}" + +[parameter] +_module = "MPIReco" +_type = "RecoPlan{SinglePatchParameters}" + + [parameter.reco] + _module = "MPIReco" + _type = "RecoPlan{SinglePatchHandsFreeReconstructionParameter}" + + [parameter.reco.sfLoad] + _module = "AbstractImageReconstruction" + _type = "RecoPlan{ProcessResultCache}" + + [parameter.reco.sfLoad.param] + _module = "MPIReco" + _type = "RecoPlan{DenseSystemMatixLoadingParameter}" + + [parameter.reco.sfLoad.param.freqFilter] + _module = "MPIReco" + _type = "RecoPlan{SNRThresholdFrequencyFilterParameter}" + + [parameter.reco.sfLoad.param.gridding] + _module = "MPIReco" + _type = "RecoPlan{SystemMatrixGriddingParameter}" + + + [[parameter.reco._listener.sf]] + field = "fov" + _module = "AbstractImageReconstruction" + _type = "LinkedFieldListener" + plan = ["parameter", "reco", "sfLoad", "param", "gridding"] + + [parameter.reco._listener.sf.fn] + _module = "MPIReco" + _type = "defaultParameterCalibFov" + [[parameter.reco._listener.sf]] + field = "center" + _module = "AbstractImageReconstruction" + _type = "LinkedFieldListener" + plan = ["parameter", "reco", "sfLoad", "param", "gridding"] + + [parameter.reco._listener.sf.fn] + _module = "MPIReco" + _type = "defaultParameterCalibCenter" + [[parameter.reco._listener.sf]] + field = "gridsize" + _module = "AbstractImageReconstruction" + _type = "LinkedFieldListener" + plan = ["parameter", "reco", "sfLoad", "param", "gridding"] + + [parameter.reco._listener.sf.fn] + _module = "MPIReco" + _type = "defaultParameterGridSize" + [[parameter.reco._listener.sf]] + field = "minFreq" + _module = "AbstractImageReconstruction" + _type = "LinkedFieldListener" + plan = ["parameter", "reco", "sfLoad", "param", "freqFilter"] + + [parameter.reco._listener.sf.fn] + _module = "MPIReco" + _type = "defaultParameterMinFreq" + [[parameter.reco._listener.sf]] + field = "maxFreq" + _module = "AbstractImageReconstruction" + _type = "LinkedFieldListener" + plan = ["parameter", "reco", "sfLoad", "param", "freqFilter"] + + [parameter.reco._listener.sf.fn] + _module = "MPIReco" + _type = "defaultParameterMaxFreq" + [[parameter.reco._listener.sf]] + field = "recChannels" + _module = "AbstractImageReconstruction" + _type = "LinkedFieldListener" + plan = ["parameter", "reco", "sfLoad", "param", "freqFilter"] + + [parameter.reco._listener.sf.fn] + _module = "MPIReco" + _type = "defaultParameterRecChannels" + + [parameter.reco.solverParams] + _module = "MPIReco" + _type = "RecoPlan{HandsFreeSolverParameters}" + + [parameter.pre] + _module = "AbstractImageReconstruction" + _type = "RecoPlan{ProcessResultCache}" + + [parameter.pre.param] + _module = "MPIReco" + _type = "RecoPlan{CommonPreProcessingParameters}" + + [parameter.post] + _module = "MPIReco" + _type = "RecoPlan{NoPostProcessing}"