From 83acb49a114644d5cd3878af1f4cb93a010fa6af Mon Sep 17 00:00:00 2001
From: nHackel <nHackel@users.noreply.github.com>
Date: Wed, 30 Aug 2023 10:29:16 +0000
Subject: [PATCH 01/22] Init updating freqs

---
 src/FrequencyFilter.jl | 59 ++++++++++++++++++++++++++++--------------
 1 file changed, 39 insertions(+), 20 deletions(-)

diff --git a/src/FrequencyFilter.jl b/src/FrequencyFilter.jl
index ad715281..a2a84170 100644
--- a/src/FrequencyFilter.jl
+++ b/src/FrequencyFilter.jl
@@ -22,30 +22,17 @@ function filterFrequencies(f::MPIFile; SNRThresh=-1, minFreq=0,
   nFreq = rxNumFrequencies(f, numPeriodGrouping)
   nReceivers = rxNumChannels(f)
   nPeriods = 1 #acqNumPeriodsPerFrame(f)
+  freqIndices = collect(vec(CartesianIndices((nFreq, nReceivers))))
 
   minIdx = floor(Int, minFreq / rxBandwidth(f) * (nFreq-1) ) + 1
   maxIdx = ceil(Int, maxFreq / rxBandwidth(f) * (nFreq-1) ) + 1
 
-  freqMask = zeros(Bool,nFreq,nReceivers,nPeriods)
-
-  freqMask[:,recChannels,:] .= true
-
-  if measIsFrequencySelection(f)
-    freqMask[:,recChannels,:] .= false
-    idx = measFrequencySelection(f)
-    freqMask[idx,recChannels,:] .= true
-  else
-    freqMask[:,recChannels,:] .= true
-  end
-
-  if minFreq > 0
-    freqMask[1:(minIdx),:,:] .= false
-  end
-
-  if maxIdx < nFreq
-    freqMask[(maxIdx):end,:,:] .= false
-  end
+  filterFrequenciesBySelection!(freqIndices, f)
+  filterFrequenciesByChannel!(freqIndices, recChannels)
+  filterFrequenciesByMinIdx!(freqIndices, minIdx)
+  filterFrequenciesByMaxIdx!(freqIndices, maxIdx)
 
+  #=
   if maxMixingOrder > 0
     if numPeriodGrouping == 1
       mf = mixingFactors(f)
@@ -71,7 +58,7 @@ function filterFrequencies(f::MPIFile; SNRThresh=-1, minFreq=0,
       freqMask .&= freqMaskMO
     end
   end
-
+  =#
   SNRAll = nothing
 
   if SNRThresh > 0 || numUsedFreqs > 0 || sortBySNR
@@ -137,6 +124,38 @@ function filterFrequencies(f::MPIFile; SNRThresh=-1, minFreq=0,
   freq
 end
 
+function filterFrequenciesBySelection!(indices::Vector{CartesianIndex}, f::MPIFile)
+  if measIsFrequencySelection(f)
+    return filterFrequenciesBySelection!(indices, measFrequencySelection(f))
+  end
+end
+filterFrequenciesBySelection!(indices::Vector{CartesianIndex}, selection::Vector{Int64}) = filter!(x -> in(x[1], selection), indices)
+
+filterFrequenciesByChannel!(indices::Vector{CartesianIndex}, channels) = filter!(x-> in(x[2], channels), indices)
+
+function filterFrequenciesByMinFreq!(indices::Vector{CartesianIndex}, f::MPIFile, minFreq)
+  nFreq = rxNumFrequencies(f, numPeriodGrouping)
+  minIdx = floor(Int, minFreq / rxBandwidth(f) * (nFreq-1) ) + 1
+  return filterFrequenciesByMinIdx!(indices, minIdx)
+end
+filterFrequenciesByMinIdx!(indices::Vector{CartesianIndex}, minIdx) = minIdx > 0 ?  filter!(x -> x[1] > minFreq, freqIndices) : indices 
+
+
+function filterFrequenciesByMaxFreq!(indices::Vector{CartesianIndex}, f::MPIFile, maxFreq)
+  nFreq = rxNumFrequencies(f, numPeriodGrouping)
+  maxIdx = ceil(Int, maxFreq / rxBandwidth(f) * (nFreq-1) ) + 1
+  return filterFrequenciesByMaxIdx!(indices, maxIdx)
+end
+filterFrequenciesByMaxIdx!(indices::Vector{CartesianIndex}, maxIdx) = filter!(x-> x[1] <= maxIdx, indices)
+
+function sortFrequenciesBySNR!()
+end
+
+
+
+function sortFrequenciesByMixFactors()
+end
+
 
 function rowsToSubsampledRows(f::MPIFile, rows)
   if measIsFrequencySelection(f)

From 3208ea9466a8e1f2b119417e25165681439170c2 Mon Sep 17 00:00:00 2001
From: nHackel <nHackel@users.noreply.github.com>
Date: Wed, 30 Aug 2023 16:01:04 +0200
Subject: [PATCH 02/22] Update filterFrequencies and seperate filtering and
 sorting

Atm features duplicate SNR computation, could be further optimized later
---
 src/FrequencyFilter.jl | 156 +++++++++++++++++++++++++----------------
 1 file changed, 97 insertions(+), 59 deletions(-)

diff --git a/src/FrequencyFilter.jl b/src/FrequencyFilter.jl
index a2a84170..280e1c0c 100644
--- a/src/FrequencyFilter.jl
+++ b/src/FrequencyFilter.jl
@@ -1,4 +1,19 @@
-export filterFrequencies
+export filterFrequencies, sortFrequencies
+
+function getCalibSNR(f::MPIFile; numPeriodGrouping = 1, stepsize = 1)
+  nFreq = rxNumFrequencies(f, numPeriodGrouping)
+  nReceivers = rxNumChannels(f)
+  SNR = zeros(nFreq, nReceivers)
+  idx = measIsFrequencySelection(f) ? measFrequencySelection(f) : idx = 1:nFreq
+
+  SNRAll = calibSNR(f)
+  if SNRAll != nothing
+    SNR[idx,:] = SNRAll[:,:,1]
+  end
+
+  SNR = SNR[1:stepsize:nFreq,:,:]
+  return SNR
+end
 
 """
   filterFrequencies(f; kargs...) => Vector{Int64}
@@ -15,9 +30,8 @@ Supported keyword arguments:
 * sortByMixFactors
 """
 function filterFrequencies(f::MPIFile; SNRThresh=-1, minFreq=0,
-               maxFreq=rxBandwidth(f), recChannels=1:rxNumChannels(f),
-               sortBySNR=false, numUsedFreqs=-1, stepsize=1, maxMixingOrder=-1,
-               sortByMixFactors=false, numPeriodAverages=1, numPeriodGrouping=1)
+               maxFreq=rxBandwidth(f), recChannels=1:rxNumChannels(f),numUsedFreqs=-1, stepsize=1,
+                maxMixingOrder=-1, numPeriodAverages=1, numPeriodGrouping=1)
 
   nFreq = rxNumFrequencies(f, numPeriodGrouping)
   nReceivers = rxNumChannels(f)
@@ -61,7 +75,7 @@ function filterFrequencies(f::MPIFile; SNRThresh=-1, minFreq=0,
   =#
   SNRAll = nothing
 
-  if SNRThresh > 0 || numUsedFreqs > 0 || sortBySNR
+  if SNRThresh > 0 || numUsedFreqs > 0
     SNR = zeros(nFreq, nReceivers)
     idx = measIsFrequencySelection(f) ? measFrequencySelection(f) : idx = 1:nFreq
 
@@ -71,90 +85,114 @@ function filterFrequencies(f::MPIFile; SNRThresh=-1, minFreq=0,
     end
   end
 
-  if SNRThresh > 0 && numUsedFreqs > 0
+  if SNRThresh > 0 && numUsedFreqs <= 0
+    filterFrequenciesBySNRThresh!(freqIndices, SNRThresh, SNR)
+  elseif  numUsedFreqs > 0 && SNRThresh <= 0
+    filterFrequenciesByNumUsedFrequencies!(freqIndices, numUsedFreqs)
+  else
     error("It is not possible to use SNRThresh and SNRFactorUsedFreq similtaneously")
   end
 
-  if SNRThresh > 0 && SNRAll != nothing
-    freqMask[ findall(vec(SNR) .< SNRThresh) ] .= false
-  end
-
-  if numUsedFreqs > 0 && SNRAll != nothing
-    numFreqsAlreadyFalse = sum(!freqMask)
-    numFreqsFalse = round(Int,length(freqMask)* (1-numUsedFreqs))
-    S = sortperm(vec(SNR))
-
-    l = 1
-    j = 1
-    while j<  (numFreqsFalse-numFreqsAlreadyFalse)
-      if freqMask[S[l]] == true
-        freqMask[S[l]] = false
-        j += 1
-      end
-      l += 1
-    end
-
-  end
-
 
   if stepsize > 1
-    freqStepsizeMask = zeros(Bool,nFreq, nReceivers, nPatches)
-    freqStepsizeMask[1:stepsize:nFreq,:,:] = freqMask[1:stepsize:nFreq,:,:]
-    freqMask = freqStepsizeMask
-  end
-
-  freq = findall( vec(freqMask) )
-
-  if sortBySNR && !sortByMixFactors && SNRAll != nothing
-    SNR = vec(SNR[1:stepsize:nFreq,:,:])
-
-    freq = freq[reverse(sortperm(SNR[freq]),dims=1)]
+    filterFrequenciesByStepsize!(freqIndices, stepsize)
   end
 
-  if !sortBySNR && sortByMixFactors
-    mfNorm = zeros(nFreq,nReceivers,nPeriods)
-    mf = mixingFactors(f)
-    for k=1:nFreq
-      mfNorm[k,:,:] = norm(mf[k,1:3])
-    end
-
-    freq = freq[sortperm(mfNorm[freq])]
-  end
-
-  freq
+  return freqIndices
 end
 
-function filterFrequenciesBySelection!(indices::Vector{CartesianIndex}, f::MPIFile)
+export filterFrequenciesBySelection!
+function filterFrequenciesBySelection!(indices::Vector{CartesianIndex{2}}, f::MPIFile)
   if measIsFrequencySelection(f)
     return filterFrequenciesBySelection!(indices, measFrequencySelection(f))
   end
 end
-filterFrequenciesBySelection!(indices::Vector{CartesianIndex}, selection::Vector{Int64}) = filter!(x -> in(x[1], selection), indices)
+filterFrequenciesBySelection!(indices::Vector{CartesianIndex{2}}, selection::Vector{Int64}) = filter!(x -> in(x[1], selection), indices)
 
-filterFrequenciesByChannel!(indices::Vector{CartesianIndex}, channels) = filter!(x-> in(x[2], channels), indices)
+export filterFrequenciesByChannel!
+filterFrequenciesByChannel!(indices::Vector{CartesianIndex{2}}, channels) = filter!(x-> in(x[2], channels), indices)
 
-function filterFrequenciesByMinFreq!(indices::Vector{CartesianIndex}, f::MPIFile, minFreq)
+export filterFrequenciesByMinFreq!
+function filterFrequenciesByMinFreq!(indices::Vector{CartesianIndex{2}}, f::MPIFile, minFreq; numPeriodGrouping = 1)
   nFreq = rxNumFrequencies(f, numPeriodGrouping)
   minIdx = floor(Int, minFreq / rxBandwidth(f) * (nFreq-1) ) + 1
   return filterFrequenciesByMinIdx!(indices, minIdx)
 end
-filterFrequenciesByMinIdx!(indices::Vector{CartesianIndex}, minIdx) = minIdx > 0 ?  filter!(x -> x[1] > minFreq, freqIndices) : indices 
+filterFrequenciesByMinIdx!(indices::Vector{CartesianIndex{2}}, minIdx) = minIdx > 0 ?  filter!(x -> x[1] > minIdx, indices) : indices 
 
-
-function filterFrequenciesByMaxFreq!(indices::Vector{CartesianIndex}, f::MPIFile, maxFreq)
+export filterFrequenciesByMaxFreq!
+function filterFrequenciesByMaxFreq!(indices::Vector{CartesianIndex{2}}, f::MPIFile, maxFreq; numPeriodGrouping = 1)
   nFreq = rxNumFrequencies(f, numPeriodGrouping)
   maxIdx = ceil(Int, maxFreq / rxBandwidth(f) * (nFreq-1) ) + 1
   return filterFrequenciesByMaxIdx!(indices, maxIdx)
 end
-filterFrequenciesByMaxIdx!(indices::Vector{CartesianIndex}, maxIdx) = filter!(x-> x[1] <= maxIdx, indices)
+filterFrequenciesByMaxIdx!(indices::Vector{CartesianIndex{2}}, maxIdx) = filter!(x-> x[1] <= maxIdx, indices)
+
+
+export filterFrequenciesBySNRThresh!
+function filterFrequenciesBySNRThresh!(indices::Vector{CartesianIndex{2}}, f::MPIFile, snrthresh; numPeriodGrouping = 1)
+  SNR = getCalibSNR(f, numPeriodGrouping = numPeriodGrouping)
+  return filterFrequenciesBySNRThresh!(indices, snrthresh, SNR)
+end
+filterFrequenciesBySNRThresh!(indices::Vector{CartesianIndex{2}}, SNRThresh, SNR) = filter!(x-> SNR[x] >= SNRThresh , indices)
 
-function sortFrequenciesBySNR!()
+export filterFrequenciesByNumUsedFrequencies!
+function filterFrequenciesByNumUsedFrequencies!(indices::Vector{CartesianIndex{2}}, f::MPIFile, maxFreq)
+  error("Not implemented")
 end
+filterFrequenciesByNumUsedFrequencies!(indices::Vector{CartesianIndex{2}}, maxIdx) = error("not implemented")
+#=
+    numFreqsAlreadyFalse = sum(!freqMask)
+    numFreqsFalse = round(Int,length(freqMask)* (1-numUsedFreqs))
+    S = sortperm(vec(SNR))
+
+    l = 1
+    j = 1
+    while j<  (numFreqsFalse-numFreqsAlreadyFalse)
+      if freqMask[S[l]] == true
+        freqMask[S[l]] = false
+        j += 1
+      end
+      l += 1
+    end
+=#
 
+export filterFrequenciesByStepsize!
+filterFrequenciesByStepsize!(indices::Vector{CartesianIndex{2}}, stepsize) = error("Not implemented")
+#=
+    freqStepsizeMask = zeros(Bool,nFreq, nReceivers, nPatches)
+    freqStepsizeMask[1:stepsize:nFreq,:,:] = freqMask[1:stepsize:nFreq,:,:]
+    freqMask = freqStepsizeMask
+=#
 
+function sortFrequencies!(indices::Vector{CartesianIndex{2}}, f::MPIFile; numPeriodGrouping = 1, stepsize = 1, sortBySNR = false, sortByMixFactors = false)
+  if sortBySNR && !sortByMixFactors
+    indices = sortFrequenciesBySNR(indices, f, numPeriodGrouping = numPeriodGrouping, stepsize = stepsize)
+  elseif !sortBySNR && sortByMixFactors
+    indices = sortFrequenciesByMixFactors(indices, f, numPeriodGrouping = numPeriodGrouping)
+  else
+    error("Can not apply multiple sorting algorithms to frequencies")
+  end
+  return indices
+end
+function sortFrequenciesBySNR(indices::Vector{CartesianIndex{2}}, f::MPIFile; numPeriodGrouping = 1, stepsize = 1)
+  SNR = getCalibSNR(f, numPeriodGrouping = numPeriodGrouping, stepsize = stepsize)
+  sortFrequenciesBySNR(indices, SNR)
+end
+sortFrequenciesBySNR(indices::Vector{CartesianIndex{2}}, SNR::AbstractArray) = indices[reverse(sortperm(SNR[indices]),dims=1)]
 
-function sortFrequenciesByMixFactors()
+function sortFrequenciesByMixFactors(indices, f::MPIFile; numPeriodGrouping = 1)
+  nFreq = rxNumFrequencies(f, numPeriodGrouping)
+  nReceivers = rxNumChannels(f)
+  nPeriods = 1 #acqNumPeriodsPerFrame(f)
+  mfNorm = zeros(nFreq,nReceivers,nPeriods)
+  mf = mixingFactors(f)
+  for k=1:nFreq
+    mfNorm[k,:,:] = norm(mf[k,1:3])
+  end
+  sortFrequenciesByMixFactors(indices, mfNorm)
 end
+sortFrequenciesByMixFactors(indices::Vector{CartesianIndex{2}}, mfNorm::AbstractArray) = indices[sortperm(mfNorm[indices])]
 
 
 function rowsToSubsampledRows(f::MPIFile, rows)

From 177b7a1ca0c4bd74c625f78f63bb2253c9158f93 Mon Sep 17 00:00:00 2001
From: nHackel <nHackel@users.noreply.github.com>
Date: Wed, 30 Aug 2023 17:06:14 +0200
Subject: [PATCH 03/22] Fix frequencyfilter throwing error if both snr filters
 are set to 0

---
 src/FrequencyFilter.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/FrequencyFilter.jl b/src/FrequencyFilter.jl
index 280e1c0c..11346d35 100644
--- a/src/FrequencyFilter.jl
+++ b/src/FrequencyFilter.jl
@@ -89,7 +89,7 @@ function filterFrequencies(f::MPIFile; SNRThresh=-1, minFreq=0,
     filterFrequenciesBySNRThresh!(freqIndices, SNRThresh, SNR)
   elseif  numUsedFreqs > 0 && SNRThresh <= 0
     filterFrequenciesByNumUsedFrequencies!(freqIndices, numUsedFreqs)
-  else
+  elseif numUsedFreqs > 0 && SNRThresh > 0
     error("It is not possible to use SNRThresh and SNRFactorUsedFreq similtaneously")
   end
 

From 955c7643176f189942eee002927d5de967689b3b Mon Sep 17 00:00:00 2001
From: nHackel <nHackel@users.noreply.github.com>
Date: Wed, 30 Aug 2023 18:49:06 +0200
Subject: [PATCH 04/22] Access FD data with cartesianindex instead of reshape

---
 src/MDFCommon.jl    | 8 ++------
 src/Measurements.jl | 1 -
 2 files changed, 2 insertions(+), 7 deletions(-)

diff --git a/src/MDFCommon.jl b/src/MDFCommon.jl
index 9777e07f..8083a320 100644
--- a/src/MDFCommon.jl
+++ b/src/MDFCommon.jl
@@ -34,9 +34,7 @@ function systemMatrix(f::Union{MDFFileV2, MDFv2InMemory}, rows, bgCorrection=tru
 
   data_ = measDataRaw(f)
 
-  data_ = reshape(data_, size(data_,1),
-          size(data_,2)*size(data_,3),
-          size(data_,4))[:, rows_, :]
+  data_ = data_[:, rows_, :]
   data = reshape(data_, Val(2))
 
   fgdata = data[measFGFrameIdx(f),:]
@@ -46,9 +44,7 @@ function systemMatrix(f::Union{MDFFileV2, MDFv2InMemory}, rows, bgCorrection=tru
     B = linearOperator(measSparsityTransformation(f), calibSize(f))
 
     tmp = measSubsamplingIndices(f)
-    subsamplingIndices_ = reshape(tmp, size(tmp,1),
-                          size(tmp,2)*size(tmp,3),
-                          size(tmp,4))[:, rows_, :]
+    subsamplingIndices_ = tmp[:, rows_, :]
     subsamplingIndices = reshape(subsamplingIndices_, Val(2))
 
     for l=1:size(fgdata,2)
diff --git a/src/Measurements.jl b/src/Measurements.jl
index be3b7e41..150c4e4f 100644
--- a/src/Measurements.jl
+++ b/src/Measurements.jl
@@ -344,7 +344,6 @@ function getMeasurementsFD(f::MPIFile, args...;
 
   if frequencies != nothing
     # here we merge frequencies and channels
-    data = reshape(data, size(data,1)*size(data,2), size(data,3), size(data,4))
     data = data[frequencies, :, :]
   end
 

From 17cdee6393211b93f6a5a0e1ccafb5938a1f37fe Mon Sep 17 00:00:00 2001
From: nHackel <nHackel@users.noreply.github.com>
Date: Wed, 30 Aug 2023 18:50:25 +0200
Subject: [PATCH 05/22] Attempt to update rowToSubsampledRows

---
 src/FrequencyFilter.jl | 13 ++++++++-----
 1 file changed, 8 insertions(+), 5 deletions(-)

diff --git a/src/FrequencyFilter.jl b/src/FrequencyFilter.jl
index 11346d35..df1a3730 100644
--- a/src/FrequencyFilter.jl
+++ b/src/FrequencyFilter.jl
@@ -198,15 +198,18 @@ sortFrequenciesByMixFactors(indices::Vector{CartesianIndex{2}}, mfNorm::Abstract
 function rowsToSubsampledRows(f::MPIFile, rows)
   if measIsFrequencySelection(f)
     # In this case we need to convert indices
-    tmp = zeros(Int64, rxNumFrequencies(f), rxNumChannels(f) )
+    tmp = Array{Union{Nothing, CartesianIndex{2}}}(undef, rxNumFrequencies(f), rxNumChannels(f))
     idxAvailable = measFrequencySelection(f)
-    for d=1:rxNumChannels(f)
-      tmp[idxAvailable, d] = (1:length(idxAvailable)) .+ (d-1)*length(idxAvailable)
+    for (i, idx) in enumerate(idxAvailable)
+      for d=1:rxNumChannels(f)
+        tmp[idx, d] = CartesianIndex(i, d)
+      end
     end
-    rows_ = vec(tmp)[rows]
-    if findfirst(x -> x == 0, rows_) != nothing
+    rows_ = tmp[rows]
+    if any(isnothing, rows_)
       @error "Indices applied to systemMatrix are not available in the file"
     end
+    return identity.(rows_) # Go from Vector{Union{Nothing, Index}} to Vector{Index}
   else
     rows_ = rows
   end

From 66751d4d94b65601cbbffcada1eee9f759d6ac75 Mon Sep 17 00:00:00 2001
From: nHackel <nHackel@users.noreply.github.com>
Date: Thu, 31 Aug 2023 11:34:36 +0200
Subject: [PATCH 06/22] Fix error where lowered bound for MinFreq was not >=

---
 src/FrequencyFilter.jl | 5 ++++-
 1 file changed, 4 insertions(+), 1 deletion(-)

diff --git a/src/FrequencyFilter.jl b/src/FrequencyFilter.jl
index df1a3730..79fb6f4e 100644
--- a/src/FrequencyFilter.jl
+++ b/src/FrequencyFilter.jl
@@ -118,7 +118,7 @@ function filterFrequenciesByMinFreq!(indices::Vector{CartesianIndex{2}}, f::MPIF
   minIdx = floor(Int, minFreq / rxBandwidth(f) * (nFreq-1) ) + 1
   return filterFrequenciesByMinIdx!(indices, minIdx)
 end
-filterFrequenciesByMinIdx!(indices::Vector{CartesianIndex{2}}, minIdx) = minIdx > 0 ?  filter!(x -> x[1] > minIdx, indices) : indices 
+filterFrequenciesByMinIdx!(indices::Vector{CartesianIndex{2}}, minIdx) = minIdx > 0 ?  filter!(x -> x[1] >= minIdx, indices) : indices 
 
 export filterFrequenciesByMaxFreq!
 function filterFrequenciesByMaxFreq!(indices::Vector{CartesianIndex{2}}, f::MPIFile, maxFreq; numPeriodGrouping = 1)
@@ -175,12 +175,15 @@ function sortFrequencies!(indices::Vector{CartesianIndex{2}}, f::MPIFile; numPer
   end
   return indices
 end
+
+export sortFrequenciesBySNR
 function sortFrequenciesBySNR(indices::Vector{CartesianIndex{2}}, f::MPIFile; numPeriodGrouping = 1, stepsize = 1)
   SNR = getCalibSNR(f, numPeriodGrouping = numPeriodGrouping, stepsize = stepsize)
   sortFrequenciesBySNR(indices, SNR)
 end
 sortFrequenciesBySNR(indices::Vector{CartesianIndex{2}}, SNR::AbstractArray) = indices[reverse(sortperm(SNR[indices]),dims=1)]
 
+export sortFrequenciesByMixFactors
 function sortFrequenciesByMixFactors(indices, f::MPIFile; numPeriodGrouping = 1)
   nFreq = rxNumFrequencies(f, numPeriodGrouping)
   nReceivers = rxNumChannels(f)

From c6b7e9336691692bda397ba04ccb463156c81b78 Mon Sep 17 00:00:00 2001
From: Jonas Schumacher <github735@jonasschumacher.de>
Date: Wed, 6 Sep 2023 11:22:10 +0200
Subject: [PATCH 07/22] Move some stuff to Base and add Base methods to Regular
 GridPositions

---
 src/Positions/Positions.jl | 16 +++++++++-------
 1 file changed, 9 insertions(+), 7 deletions(-)

diff --git a/src/Positions/Positions.jl b/src/Positions/Positions.jl
index 8c5af8f5..0fa98702 100644
--- a/src/Positions/Positions.jl
+++ b/src/Positions/Positions.jl
@@ -80,6 +80,8 @@ function range(grid::RegularGridPositions, dim::Int)
     return 1:1
   end
 end
+Base.ndims(grid::RegularGridPositions) = length(grid.shape)
+Base.axes(grid::RegularGridPositions) = tuple([range(grid, i) for i in 1:ndims(grid)]...)
 
 RegularGridPositions(shape, fov, center) = RegularGridPositions(shape, fov, center, ones(Int,length(shape)))
 
@@ -670,17 +672,17 @@ function Base.:(==)(val1::Positions, val2::Positions)
 end
 
 # fuction related to looping
-length(tdes::SphericalTDesign) = size(tdes.positions,2)
-length(apos::ArbitraryPositions) = size(apos.positions,2)
-length(grid::GridPositions) = prod(grid.shape)
-length(rpos::UniformRandomPositions) = rpos.N
-length(mgrid::MeanderingGridPositions) = length(mgrid.grid)
-length(bgrid::BreakpointGridPositions) = length(bgrid.grid)+length(bgrid.breakpointIndices)
+Base.length(tdes::SphericalTDesign) = size(tdes.positions,2)
+Base.length(apos::ArbitraryPositions) = size(apos.positions,2)
+Base.length(grid::GridPositions) = prod(grid.shape)
+Base.length(rpos::UniformRandomPositions) = rpos.N
+Base.length(mgrid::MeanderingGridPositions) = length(mgrid.grid)
+Base.length(bgrid::BreakpointGridPositions) = length(bgrid.grid)+length(bgrid.breakpointIndices)
 
 start_(grid::Positions) = 1
 next_(grid::Positions,state) = (grid[state],state+1)
 done_(grid::Positions,state) = state > length(grid)
-iterate(grid::Positions, s=start_(grid)) = done_(grid, s) ? nothing : next_(grid, s)
+Base.iterate(grid::Positions, s=start_(grid)) = done_(grid, s) ? nothing : next_(grid, s)
 
 
 include("Interpolation.jl")

From 1d8e8806d0172e8361490f9d30ade1152b63ccb2 Mon Sep 17 00:00:00 2001
From: Jonas Schumacher <github735@jonasschumacher.de>
Date: Mon, 11 Sep 2023 14:09:40 +0200
Subject: [PATCH 08/22] Make RegularGridPositions work with CartesianIndices

---
 src/Positions/Positions.jl | 4 ++++
 1 file changed, 4 insertions(+)

diff --git a/src/Positions/Positions.jl b/src/Positions/Positions.jl
index 0fa98702..74603e27 100644
--- a/src/Positions/Positions.jl
+++ b/src/Positions/Positions.jl
@@ -211,6 +211,10 @@ function getindex(grid::RegularGridPositions, idx::Vector{T}) where T<:Number
   return 0.5.*fieldOfView(grid) .* (-1 .+ (2 .* idx .- 1) ./ shape(grid)) .+ fieldOfViewCenter(grid)
 end
 
+function getindex(grid::RegularGridPositions, idx::CartesianIndex)
+  return getindex(grid, LinearIndices(tuple(grid.shape...))[idx])
+end
+
 function posToIdxFloat(grid::RegularGridPositions,pos::Vector)
   idx = 0.5 .* (shape(grid) .* ((pos .- fieldOfViewCenter(grid)) ./
               ( 0.5 .* fieldOfView(grid) ) .+ 1) .+ 1)

From ee93455f2b581726d74fd2fe42635dd2990853df Mon Sep 17 00:00:00 2001
From: Jonas Schumacher <github735@jonasschumacher.de>
Date: Thu, 14 Sep 2023 08:15:02 +0200
Subject: [PATCH 09/22] Fix plural

---
 docs/src/lowlevel.md | 2 +-
 src/Brukerfile.jl    | 2 +-
 src/Conversion.jl    | 4 ++--
 src/IMT.jl           | 2 +-
 src/MDF.jl           | 2 +-
 src/MDFInMemory.jl   | 4 ----
 src/MPIFiles.jl      | 2 +-
 src/MultiMPIFile.jl  | 2 +-
 test/IMT.jl          | 2 +-
 test/MDFv1.jl        | 2 +-
 10 files changed, 10 insertions(+), 14 deletions(-)

diff --git a/docs/src/lowlevel.md b/docs/src/lowlevel.md
index 5399dee8..9af862f4 100644
--- a/docs/src/lowlevel.md
+++ b/docs/src/lowlevel.md
@@ -40,7 +40,7 @@ measIsBGFrame, measIsSpectralLeakageCorrected, measFramePermutation
 
 # calibrations
 calibSNR, calibFov, calibFovCenter, calibSize, calibOrder, calibPositions,
-calibOffsetField, calibDeltaSampleSize, calibMethod, calibIsMeanderingGrid
+calibOffsetFields, calibDeltaSampleSize, calibMethod, calibIsMeanderingGrid
 
 # reconstruction results
 recoData, recoFov, recoFovCenter, recoSize, recoOrder, recoPositions
diff --git a/src/Brukerfile.jl b/src/Brukerfile.jl
index e7cc1438..80785249 100644
--- a/src/Brukerfile.jl
+++ b/src/Brukerfile.jl
@@ -659,7 +659,7 @@ calibFovCenter(b::BrukerFile) =
 calibSize(b::BrukerFile) = [parse(Int64,s) for s in b["PVM_Matrix"]]
 calibOrder(b::BrukerFile) = "xyz"
 calibPositions(b::BrukerFile) = nothing
-calibOffsetField(b::BrukerFile) = nothing
+calibOffsetFields(b::BrukerFile) = nothing
 calibDeltaSampleSize(b::BrukerFile) = [0.0, 0.0, 0.0] #TODO
 calibMethod(b::BrukerFile) = "robot"
 calibIsMeanderingGrid(b::BrukerFile) = true
diff --git a/src/Conversion.jl b/src/Conversion.jl
index eb565a0c..bee5fee3 100644
--- a/src/Conversion.jl
+++ b/src/Conversion.jl
@@ -115,7 +115,7 @@ end
 function loadCalibParams(f, params = Dict{Symbol,Any}())
   if experimentIsCalibration(f)
     for op in [:calibFov, :calibFovCenter,
-               :calibSize, :calibOrder, :calibPositions, :calibOffsetField,
+               :calibSize, :calibOrder, :calibPositions, :calibOffsetFields,
                :calibDeltaSampleSize, :calibMethod]
       setparam!(params, op, eval(op)(f))
     end
@@ -825,7 +825,7 @@ function saveasMDF(file::HDF5.File, params::Dict{Symbol,Any})
   writeIfAvailable(file, "/calibration/size",  params, :calibSize)
   writeIfAvailable(file, "/calibration/order",  params, :calibOrder)
   writeIfAvailable(file, "/calibration/positions",  params, :calibPositions)
-  writeIfAvailable(file, "/calibration/offsetField",  params, :calibOffsetField)
+  writeIfAvailable(file, "/calibration/offsetFields",  params, :calibOffsetFields)
   writeIfAvailable(file, "/calibration/deltaSampleSize",  params, :calibDeltaSampleSize)
   writeIfAvailable(file, "/calibration/method",  params, :calibMethod)
   if hasKeyAndValue(params, :calibIsMeanderingGrid)
diff --git a/src/IMT.jl b/src/IMT.jl
index 2bacaabf..adb7e150 100644
--- a/src/IMT.jl
+++ b/src/IMT.jl
@@ -201,7 +201,7 @@ calibSize(f::IMTFile) = [div(size(f["/systemResponseFrequencies"],1),2),
                          size(f["/systemResponseFrequencies"],2),
                          size(f["/systemResponseFrequencies"],3)]
 calibOrder(f::IMTFile) = "xyz"
-calibOffsetField(f::IMTFileCalib) = nothing
+calibOffsetFields(f::IMTFileCalib) = nothing
 calibDeltaSampleSize(f::IMTFile) = [0.0,0.0,0.0]
 calibMethod(f::IMTFile) = "simulation"
 #calibIsMeanderingGrid(f::IMTFile) = Bool(f["/calibration/isMeanderingGrid", 0])
diff --git a/src/MDF.jl b/src/MDF.jl
index 38a19e64..f4a22856 100644
--- a/src/MDF.jl
+++ b/src/MDF.jl
@@ -425,7 +425,7 @@ calibFov(f::MDFFile) = @keyoptional f["/calibration/fieldOfView"]
 calibFovCenter(f::MDFFile) = @keyoptional f["/calibration/fieldOfViewCenter"]
 calibSize(f::MDFFile) = @keyoptional f["/calibration/size"]
 calibOrder(f::MDFFile) = @keyoptional f["/calibration/order"]
-calibOffsetField(f::MDFFile) = @keyoptional f["/calibration/offsetField"]
+calibOffsetFields(f::MDFFile) = @keyoptional f["/calibration/offsetFields"]
 calibDeltaSampleSize(f::MDFFile) = @keyoptional f["/calibration/deltaSampleSize",[0.0,0.0,0.0]]
 calibMethod(f::MDFFile) = @keyrequired f["/calibration/method"]
 calibIsMeanderingGrid(f::MDFFile) = @keyoptional Bool(f["/calibration/isMeanderingGrid", 0])
diff --git a/src/MDFInMemory.jl b/src/MDFInMemory.jl
index c8c1b3a3..8eb7adb0 100644
--- a/src/MDFInMemory.jl
+++ b/src/MDFInMemory.jl
@@ -1303,10 +1303,6 @@ for (fieldname, fieldtype) in zip(fieldnames(MDFv2InMemory), fieldtypes(MDFv2InM
   end
 end
 
-
-# The MDF specification uses the plural
-calibOffsetFields(f::MDFFile) = calibOffsetField(f)
-
 # And some utility functions
 measIsCalibProcessed(mdf::MDFv2InMemory)::Union{Bool, Missing} = measIsFramePermutation(mdf) && 
                                                                  measIsFourierTransformed(mdf) &&
diff --git a/src/MPIFiles.jl b/src/MPIFiles.jl
index 3405d97c..ed80113a 100644
--- a/src/MPIFiles.jl
+++ b/src/MPIFiles.jl
@@ -199,7 +199,7 @@ abstract type MPIFile end
 @mustimplement calibSize(f::MPIFile)
 @mustimplement calibOrder(f::MPIFile)
 @mustimplement calibPositions(f::MPIFile)
-@mustimplement calibOffsetField(f::MPIFile)
+@mustimplement calibOffsetFields(f::MPIFile)
 @mustimplement calibDeltaSampleSize(f::MPIFile)
 @mustimplement calibMethod(f::MPIFile)
 @mustimplement calibIsMeanderingGrid(f::MPIFile)
diff --git a/src/MultiMPIFile.jl b/src/MultiMPIFile.jl
index e682e6ca..acde4ab2 100644
--- a/src/MultiMPIFile.jl
+++ b/src/MultiMPIFile.jl
@@ -144,7 +144,7 @@ end
 # TODO: define functions for multi calibration data
 #  if experimentIsCalibration(f)
 #    for op in [:calibSNR, :calibFov, :calibFovCenter,
-#               :calibSize, :calibOrder, :calibPositions, :calibOffsetField,
+#               :calibSize, :calibOrder, :calibPositions, :calibOffsetFields,
 #               :calibDeltaSampleSize, :calibMethod]
 #      setparam!(params, string(op), eval(op)(f))
 #    end
diff --git a/test/IMT.jl b/test/IMT.jl
index 1395b129..2a150681 100644
--- a/test/IMT.jl
+++ b/test/IMT.jl
@@ -93,7 +93,7 @@ end
   @test calibSize(calibIMT) == [20; 20; 1]
   @test calibOrder(calibIMT) == "xyz"
   @test calibPositions(calibIMT) == nothing
-  @test calibOffsetField(calibIMT) == nothing
+  @test calibOffsetFields(calibIMT) == nothing
   @test calibDeltaSampleSize(calibIMT) == [0.0; 0.0; 0.0]
   @test calibMethod(calibIMT) == "simulation"
 
diff --git a/test/MDFv1.jl b/test/MDFv1.jl
index 11f178ca..a5e275a2 100644
--- a/test/MDFv1.jl
+++ b/test/MDFv1.jl
@@ -125,7 +125,7 @@ for sm in (smv1,smv2)
   @test calibSize(sm) == [44; 44; 1]
   @test calibOrder(sm) == "xyz"
   @test calibPositions(smv1) == nothing
-  @test calibOffsetField(smv1) == nothing
+  @test calibOffsetFields(smv1) == nothing
   @test calibDeltaSampleSize(sm) == [0.001; 0.001; 0.001]
   @test calibMethod(sm) == "robot"
 

From d6b110eac259fc5a3207e974248e96ead55ad6d0 Mon Sep 17 00:00:00 2001
From: Jonas Schumacher <github735@jonasschumacher.de>
Date: Thu, 14 Sep 2023 08:21:07 +0200
Subject: [PATCH 10/22] Missed export

---
 src/MPIFiles.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/MPIFiles.jl b/src/MPIFiles.jl
index ed80113a..9b0df446 100644
--- a/src/MPIFiles.jl
+++ b/src/MPIFiles.jl
@@ -90,7 +90,7 @@ export measData, measDataTDPeriods, measIsFourierTransformed, measIsTFCorrected,
 # calibrations
 export calibSNR, calibSnr, calibFov, calibFieldOfView, calibFovCenter,
   calibFieldOfViewCenter, calibSize, calibOrder, calibPositions,
-  calibOffsetField, calibDeltaSampleSize,
+  calibOffsetFields, calibDeltaSampleSize,
   calibMethod, calibIsMeanderingGrid
 
 # reconstruction results

From 895edb506fdc44ac5c1c2a1e52b97f2cd6ee0b52 Mon Sep 17 00:00:00 2001
From: Jonas Schumacher <github735@jonasschumacher.de>
Date: Thu, 14 Sep 2023 08:49:31 +0200
Subject: [PATCH 11/22] Remove Base prefix since functions are already imported

---
 src/Positions/Positions.jl | 18 +++++++++---------
 1 file changed, 9 insertions(+), 9 deletions(-)

diff --git a/src/Positions/Positions.jl b/src/Positions/Positions.jl
index 74603e27..0d0f2c57 100644
--- a/src/Positions/Positions.jl
+++ b/src/Positions/Positions.jl
@@ -80,8 +80,8 @@ function range(grid::RegularGridPositions, dim::Int)
     return 1:1
   end
 end
-Base.ndims(grid::RegularGridPositions) = length(grid.shape)
-Base.axes(grid::RegularGridPositions) = tuple([range(grid, i) for i in 1:ndims(grid)]...)
+ndims(grid::RegularGridPositions) = length(grid.shape)
+axes(grid::RegularGridPositions) = tuple([range(grid, i) for i in 1:ndims(grid)]...)
 
 RegularGridPositions(shape, fov, center) = RegularGridPositions(shape, fov, center, ones(Int,length(shape)))
 
@@ -676,17 +676,17 @@ function Base.:(==)(val1::Positions, val2::Positions)
 end
 
 # fuction related to looping
-Base.length(tdes::SphericalTDesign) = size(tdes.positions,2)
-Base.length(apos::ArbitraryPositions) = size(apos.positions,2)
-Base.length(grid::GridPositions) = prod(grid.shape)
-Base.length(rpos::UniformRandomPositions) = rpos.N
-Base.length(mgrid::MeanderingGridPositions) = length(mgrid.grid)
-Base.length(bgrid::BreakpointGridPositions) = length(bgrid.grid)+length(bgrid.breakpointIndices)
+length(tdes::SphericalTDesign) = size(tdes.positions,2)
+length(apos::ArbitraryPositions) = size(apos.positions,2)
+length(grid::GridPositions) = prod(grid.shape)
+length(rpos::UniformRandomPositions) = rpos.N
+length(mgrid::MeanderingGridPositions) = length(mgrid.grid)
+length(bgrid::BreakpointGridPositions) = length(bgrid.grid)+length(bgrid.breakpointIndices)
 
 start_(grid::Positions) = 1
 next_(grid::Positions,state) = (grid[state],state+1)
 done_(grid::Positions,state) = state > length(grid)
-Base.iterate(grid::Positions, s=start_(grid)) = done_(grid, s) ? nothing : next_(grid, s)
+iterate(grid::Positions, s=start_(grid)) = done_(grid, s) ? nothing : next_(grid, s)
 
 
 include("Interpolation.jl")

From 37ad2cdd43c69e8bb6b7e3a501762923e9ae8901 Mon Sep 17 00:00:00 2001
From: Jonas Schumacher <github735@jonasschumacher.de>
Date: Thu, 14 Sep 2023 09:05:04 +0200
Subject: [PATCH 12/22] Add tests

---
 test/Positions.jl | 4 ++++
 1 file changed, 4 insertions(+)

diff --git a/test/Positions.jl b/test/Positions.jl
index 53f97dcd..7a351a3b 100644
--- a/test/Positions.jl
+++ b/test/Positions.jl
@@ -6,6 +6,8 @@ pospath = joinpath(tmpdir,"positions","Positions.h5")
   ctr = [0.0,0.0,0.0]Unitful.mm
   caG = RegularGridPositions(shp,fov,ctr)
   @test shape(caG) == shp
+  @test ndims(caG) == 3
+  @test axes(cag) == ((-1.0:1.0:1.0)u"mm", (-1.0:1.0:1.0)u"mm", (-1.0:1.0:1.0)u"mm")
   @test fieldOfView(caG) == fov
   @test fieldOfViewCenter(caG) == ctr
   @test_throws BoundsError caG[0]
@@ -15,6 +17,8 @@ pospath = joinpath(tmpdir,"positions","Positions.h5")
   @test caG[3] == [1,-1,-1]Unitful.mm
   @test caG[4] == [-1,0,-1]Unitful.mm
   @test caG[27] == [1,1,1]Unitful.mm
+  @test getindex(caG, CartesianIndex(1,1,1)) == [-1,-1,-1]Unitful.mm
+  @test getindex(caG, CartesianIndex(3,3,3)) == [1,1,1]Unitful.mm
   h5open(pospath, "w") do file
     write(file, caG)
   end

From f7fc1df46764b98cf67780bf43fb6ca1df68d7d2 Mon Sep 17 00:00:00 2001
From: Jonas Schumacher <github735@jonasschumacher.de>
Date: Thu, 14 Sep 2023 09:12:59 +0200
Subject: [PATCH 13/22] Fix tests

---
 test/Positions.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/test/Positions.jl b/test/Positions.jl
index 7a351a3b..638a60cb 100644
--- a/test/Positions.jl
+++ b/test/Positions.jl
@@ -7,7 +7,7 @@ pospath = joinpath(tmpdir,"positions","Positions.h5")
   caG = RegularGridPositions(shp,fov,ctr)
   @test shape(caG) == shp
   @test ndims(caG) == 3
-  @test axes(cag) == ((-1.0:1.0:1.0)u"mm", (-1.0:1.0:1.0)u"mm", (-1.0:1.0:1.0)u"mm")
+  @test axes(caG) == ((-1.0:1.0:1.0)u"mm", (-1.0:1.0:1.0)u"mm", (-1.0:1.0:1.0)u"mm")
   @test fieldOfView(caG) == fov
   @test fieldOfViewCenter(caG) == ctr
   @test_throws BoundsError caG[0]

From 2c4c74048c54c027a5446adb2002669fce0308bd Mon Sep 17 00:00:00 2001
From: Jonas Schumacher <github735@jonasschumacher.de>
Date: Thu, 14 Sep 2023 09:26:07 +0200
Subject: [PATCH 14/22] Bump version

---
 Project.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Project.toml b/Project.toml
index b49a0607..937482c3 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "MPIFiles"
 uuid = "371237a9-e6c1-5201-9adb-3d8cfa78fa9f"
 authors = ["Tobias Knopp <tobias@knoppweb.de>"]
-version = "0.12.8"
+version = "0.12.9"
 
 [deps]
 AxisArrays = "39de3d68-74b9-583c-8d2d-e117c070f3a9"

From e8b026b4adfce3478276d4d5a34da521abc6ca41 Mon Sep 17 00:00:00 2001
From: nHackel <nHackel@users.noreply.github.com>
Date: Thu, 14 Sep 2023 14:37:05 +0200
Subject: [PATCH 15/22] Implement frequency filter for mixing order and num
 sidebands

---
 src/FrequencyFilter.jl | 61 +++++++++++++++++++++++-------------------
 1 file changed, 34 insertions(+), 27 deletions(-)

diff --git a/src/FrequencyFilter.jl b/src/FrequencyFilter.jl
index 79fb6f4e..266d8668 100644
--- a/src/FrequencyFilter.jl
+++ b/src/FrequencyFilter.jl
@@ -31,7 +31,7 @@ Supported keyword arguments:
 """
 function filterFrequencies(f::MPIFile; SNRThresh=-1, minFreq=0,
                maxFreq=rxBandwidth(f), recChannels=1:rxNumChannels(f),numUsedFreqs=-1, stepsize=1,
-                maxMixingOrder=-1, numPeriodAverages=1, numPeriodGrouping=1)
+                maxMixingOrder=-1, numPeriodAverages=1, numPeriodGrouping=1, numSidebandFreqs = -1)
 
   nFreq = rxNumFrequencies(f, numPeriodGrouping)
   nReceivers = rxNumChannels(f)
@@ -46,33 +46,18 @@ function filterFrequencies(f::MPIFile; SNRThresh=-1, minFreq=0,
   filterFrequenciesByMinIdx!(freqIndices, minIdx)
   filterFrequenciesByMaxIdx!(freqIndices, maxIdx)
 
-  #=
   if maxMixingOrder > 0
     if numPeriodGrouping == 1
-      mf = mixingFactors(f)
-      for l=1:size(mf,1)
-        if mf[l,4] > maxMixingOrder || mf[l,4] > maxMixingOrder
-          freqMask[(l-1)+1,recChannels] .= false
-        end
-      end
-    else # This is a hack until we have a general solution
-
-      numPatches = div(acqNumPeriodsPerFrame(f), numPeriodAverages)
-      fBands = (collect(0:(rxNumFrequencies(f)-1)).-1) .* numPatches
-      freqMaskMO = zeros(Bool,nFreq,nReceivers,nPeriods)
-
-      for i=1:length(fBands)
-        for y=-maxMixingOrder:maxMixingOrder
-          index = fBands[i]+y
-          if 1 <= index <= size(freqMask,1)
-            freqMaskMO[index,recChannels] .= true
-          end
-        end
-      end
-      freqMask .&= freqMaskMO
+      filterFrequenciesByMaxMixingOrder!(freqIndices, maxMixingOrder, f)
+    else
+      error("Can not apply max mixing order with a period grouping larger than one, found: $numPeriodGrouping")
     end
   end
-  =#
+
+  if numSidebandFreqs > 0 && numPeriodGrouping > 1
+    filterFrequenciesByNumSidebandFreqs!(freqIndices, numSidebandFreqs, f, numPeriodGrouping = numPeriodGrouping)
+  end
+
   SNRAll = nothing
 
   if SNRThresh > 0 || numUsedFreqs > 0
@@ -118,7 +103,7 @@ function filterFrequenciesByMinFreq!(indices::Vector{CartesianIndex{2}}, f::MPIF
   minIdx = floor(Int, minFreq / rxBandwidth(f) * (nFreq-1) ) + 1
   return filterFrequenciesByMinIdx!(indices, minIdx)
 end
-filterFrequenciesByMinIdx!(indices::Vector{CartesianIndex{2}}, minIdx) = minIdx > 0 ?  filter!(x -> x[1] >= minIdx, indices) : indices 
+filterFrequenciesByMinIdx!(indices::Vector{CartesianIndex{2}}, minIdx) = minIdx > 0 ?  filter!(x -> x[1] > minIdx, indices) : indices 
 
 export filterFrequenciesByMaxFreq!
 function filterFrequenciesByMaxFreq!(indices::Vector{CartesianIndex{2}}, f::MPIFile, maxFreq; numPeriodGrouping = 1)
@@ -126,8 +111,30 @@ function filterFrequenciesByMaxFreq!(indices::Vector{CartesianIndex{2}}, f::MPIF
   maxIdx = ceil(Int, maxFreq / rxBandwidth(f) * (nFreq-1) ) + 1
   return filterFrequenciesByMaxIdx!(indices, maxIdx)
 end
-filterFrequenciesByMaxIdx!(indices::Vector{CartesianIndex{2}}, maxIdx) = filter!(x-> x[1] <= maxIdx, indices)
-
+filterFrequenciesByMaxIdx!(indices::Vector{CartesianIndex{2}}, maxIdx) = filter!(x-> x[1] < maxIdx, indices)
+
+export filterFrequenciesByMaxMixingOrder!
+filterFrequenciesByMaxMixingOrder!(indices::Vector{CartesianIndex{2}}, maxMixingOrder, f::MPIFile) = filterFrequenciesByMaxMixingOrder!(indices, maxMixingOrder, mixingFactors(f))
+filterFrequenciesByMaxMixingOrder!(indices::Vector{CartesianIndex{2}}, maxMixingOrder, mf::Matrix) = filter!(x-> mf[x[1], 4] <= maxMixingOrder, indices)
+
+export filterFrequenciesByNumSidebandFreqs!
+function filterFrequenciesByNumSidebandFreqs!(indices::Vector{CartesianIndex{2}}, numSidebandFreqs::Int64, f::MPIFile; numPeriodGrouping = 1)
+  # https://en.wikipedia.org/wiki/Sideband
+  # Because of period grouping we have more frequencies than in original data
+
+  # Find "virtual" frequency indices of original frequencies
+  fBands = (collect(0:(rxNumFrequencies(f))).-1) .* numPeriodGrouping
+
+  delete = Int64[]
+  for (i, cart) in enumerate(indices)
+    freq = cart[1]
+    # Check if there is no frequency that is at most numSideBandFreqs away from our original frequency
+    if !any(fBand -> abs(fBand - freq) <= numSidebandFreqs, fBands)
+      push!(delete, i)
+    end
+  end
+  deleteat!(indices, delete)
+end
 
 export filterFrequenciesBySNRThresh!
 function filterFrequenciesBySNRThresh!(indices::Vector{CartesianIndex{2}}, f::MPIFile, snrthresh; numPeriodGrouping = 1)

From 6a9ec4106684d47d801e940d2b3ada9e20ce8bed Mon Sep 17 00:00:00 2001
From: nHackel <nHackel@users.noreply.github.com>
Date: Thu, 14 Sep 2023 15:42:07 +0200
Subject: [PATCH 16/22] Add stepwise freq filtering

---
 src/FrequencyFilter.jl | 17 ++++++-----------
 1 file changed, 6 insertions(+), 11 deletions(-)

diff --git a/src/FrequencyFilter.jl b/src/FrequencyFilter.jl
index 266d8668..66cee64e 100644
--- a/src/FrequencyFilter.jl
+++ b/src/FrequencyFilter.jl
@@ -38,13 +38,10 @@ function filterFrequencies(f::MPIFile; SNRThresh=-1, minFreq=0,
   nPeriods = 1 #acqNumPeriodsPerFrame(f)
   freqIndices = collect(vec(CartesianIndices((nFreq, nReceivers))))
 
-  minIdx = floor(Int, minFreq / rxBandwidth(f) * (nFreq-1) ) + 1
-  maxIdx = ceil(Int, maxFreq / rxBandwidth(f) * (nFreq-1) ) + 1
-
   filterFrequenciesBySelection!(freqIndices, f)
   filterFrequenciesByChannel!(freqIndices, recChannels)
-  filterFrequenciesByMinIdx!(freqIndices, minIdx)
-  filterFrequenciesByMaxIdx!(freqIndices, maxIdx)
+  filterFrequenciesByMinFreq!(freqIndices, f, minFreq)
+  filterFrequenciesByMaxFreq!(freqIndices, f, maxFreq)
 
   if maxMixingOrder > 0
     if numPeriodGrouping == 1
@@ -165,12 +162,10 @@ filterFrequenciesByNumUsedFrequencies!(indices::Vector{CartesianIndex{2}}, maxId
 =#
 
 export filterFrequenciesByStepsize!
-filterFrequenciesByStepsize!(indices::Vector{CartesianIndex{2}}, stepsize) = error("Not implemented")
-#=
-    freqStepsizeMask = zeros(Bool,nFreq, nReceivers, nPatches)
-    freqStepsizeMask[1:stepsize:nFreq,:,:] = freqMask[1:stepsize:nFreq,:,:]
-    freqMask = freqStepsizeMask
-=#
+function filterFrequenciesByStepsize!(indices::Vector{CartesianIndex{2}}, stepsize)
+  stepIndices = 1:stepsize:maximum(map(x -> x[1], indices))
+  filter!(x -> in(x[1], stepIndices), indices)
+end
 
 function sortFrequencies!(indices::Vector{CartesianIndex{2}}, f::MPIFile; numPeriodGrouping = 1, stepsize = 1, sortBySNR = false, sortByMixFactors = false)
   if sortBySNR && !sortByMixFactors

From e2c9118da9dd2f2007aa0c136eb3fc0f2c6d64e9 Mon Sep 17 00:00:00 2001
From: nHackel <nHackel@users.noreply.github.com>
Date: Thu, 14 Sep 2023 17:19:23 +0200
Subject: [PATCH 17/22] Increase version to 0.13

---
 Project.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Project.toml b/Project.toml
index 937482c3..163dab2d 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "MPIFiles"
 uuid = "371237a9-e6c1-5201-9adb-3d8cfa78fa9f"
 authors = ["Tobias Knopp <tobias@knoppweb.de>"]
-version = "0.12.9"
+version = "0.13.0"
 
 [deps]
 AxisArrays = "39de3d68-74b9-583c-8d2d-e117c070f3a9"

From ee44bd2ff6968314b462ba2283d2f519b7ab68ea Mon Sep 17 00:00:00 2001
From: Jonas Schumacher <jonas.schumacher@imte.fraunhofer.de>
Date: Mon, 18 Sep 2023 13:03:04 +0200
Subject: [PATCH 18/22] Prevent `NaN` errors with 2D grids

---
 src/Positions/Positions.jl | 1 +
 1 file changed, 1 insertion(+)

diff --git a/src/Positions/Positions.jl b/src/Positions/Positions.jl
index 0d0f2c57..1b30c439 100644
--- a/src/Positions/Positions.jl
+++ b/src/Positions/Positions.jl
@@ -218,6 +218,7 @@ end
 function posToIdxFloat(grid::RegularGridPositions,pos::Vector)
   idx = 0.5 .* (shape(grid) .* ((pos .- fieldOfViewCenter(grid)) ./
               ( 0.5 .* fieldOfView(grid) ) .+ 1) .+ 1)
+  idx = [isnan(val) ? one(eltype(idx)) : val for val in idx]
   return idx
 end
 

From 157cf588266533a4c274b03de10442976b9a989c Mon Sep 17 00:00:00 2001
From: Jonas Schumacher <jonas.schumacher@imte.fraunhofer.de>
Date: Mon, 18 Sep 2023 15:51:05 +0200
Subject: [PATCH 19/22] Remove Vector requirement for position functions

---
 src/Positions/Positions.jl | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/src/Positions/Positions.jl b/src/Positions/Positions.jl
index 1b30c439..57b3ab6f 100644
--- a/src/Positions/Positions.jl
+++ b/src/Positions/Positions.jl
@@ -215,14 +215,14 @@ function getindex(grid::RegularGridPositions, idx::CartesianIndex)
   return getindex(grid, LinearIndices(tuple(grid.shape...))[idx])
 end
 
-function posToIdxFloat(grid::RegularGridPositions,pos::Vector)
+function posToIdxFloat(grid::RegularGridPositions, pos)
   idx = 0.5 .* (shape(grid) .* ((pos .- fieldOfViewCenter(grid)) ./
               ( 0.5 .* fieldOfView(grid) ) .+ 1) .+ 1)
   idx = [isnan(val) ? one(eltype(idx)) : val for val in idx]
   return idx
 end
 
-function posToIdx(grid::RegularGridPositions,pos::Vector)
+function posToIdx(grid::RegularGridPositions, pos)
   idx = round.(Int64, posToIdxFloat(grid,pos))
   for d=1:length(idx)
     if grid.sign[d] == -1
@@ -232,7 +232,7 @@ function posToIdx(grid::RegularGridPositions,pos::Vector)
   return idx
 end
 
-function posToLinIdx(grid::RegularGridPositions,pos::Vector)
+function posToLinIdx(grid::RegularGridPositions, pos)
   return (LinearIndices(tuple(shape(grid)...)))[posToIdx(grid,pos)...]
 end
 

From 7b575b0e7c2d7a009a47a1862ce4d56b9e95cbf2 Mon Sep 17 00:00:00 2001
From: nHackel <nHackel@users.noreply.github.com>
Date: Tue, 19 Sep 2023 17:21:12 +0200
Subject: [PATCH 20/22] Fix error with unwanted filtering of min and max freq

---
 src/FrequencyFilter.jl | 10 ++++++++--
 1 file changed, 8 insertions(+), 2 deletions(-)

diff --git a/src/FrequencyFilter.jl b/src/FrequencyFilter.jl
index 66cee64e..0bd1b4f8 100644
--- a/src/FrequencyFilter.jl
+++ b/src/FrequencyFilter.jl
@@ -40,8 +40,14 @@ function filterFrequencies(f::MPIFile; SNRThresh=-1, minFreq=0,
 
   filterFrequenciesBySelection!(freqIndices, f)
   filterFrequenciesByChannel!(freqIndices, recChannels)
-  filterFrequenciesByMinFreq!(freqIndices, f, minFreq)
-  filterFrequenciesByMaxFreq!(freqIndices, f, maxFreq)
+  
+  if minFreq > 0
+    filterFrequenciesByMinFreq!(freqIndices, f, minFreq)
+  end
+  
+  if maxFreq < nFreq
+    filterFrequenciesByMaxFreq!(freqIndices, f, maxFreq)
+  end
 
   if maxMixingOrder > 0
     if numPeriodGrouping == 1

From 9b512bb8c6f39fb029f94b2ca4b866c6a21eefc4 Mon Sep 17 00:00:00 2001
From: nHackel <nHackel@users.noreply.github.com>
Date: Tue, 19 Sep 2023 17:21:34 +0200
Subject: [PATCH 21/22] Fix brukerfile handling of cartesian index frequencies

---
 src/Brukerfile.jl | 15 +++++----------
 1 file changed, 5 insertions(+), 10 deletions(-)

diff --git a/src/Brukerfile.jl b/src/Brukerfile.jl
index 80785249..d01eda03 100644
--- a/src/Brukerfile.jl
+++ b/src/Brukerfile.jl
@@ -540,22 +540,17 @@ function systemMatrix(b::BrukerFileCalib, rows, bgCorrection=true)
   nFreq = div(rxNumSamplingPoints(b)*numSubPeriods(b),2)+1
 
   if numSubPeriods(b) > 1
-    rows_ = collect(rows)
-    NFreq = rxNumFrequencies(b)
-    NRx = rxNumChannels(b)
     stepsize = numSubPeriods(b)
     for k=1:length(rows)
-      freq = mod1(rows[k],NFreq)
-      rec = div(rows[k]-1,NFreq)
-      rows_[k] = (freq-1)*stepsize+1 + rec*nFreq
+      freq = rows[k][1]
+      rec = rows[k][2]
+      rows[k] = CartesianIndex((freq-1)*stepsize+1, rec)
     end
-  else
-    rows_ = rows
   end
 
   s = open(sfFilename)
-  data = Mmap.mmap(s, Array{ComplexF64,2}, (prod(calibSize(b)),nFreq*rxNumChannels(b)))
-  S = data[:,rows_]
+  data = Mmap.mmap(s, Array{ComplexF64,3}, (prod(calibSize(b)), nFreq, rxNumChannels(b)))
+  S = data[:,rows]
   close(s)
   rmul!(S, 1.0/acqNumAverages(b))
   return S

From 0807542f0cd630c9ad92061261b2dbac76c82934 Mon Sep 17 00:00:00 2001
From: nHackel <nHackel@users.noreply.github.com>
Date: Tue, 19 Sep 2023 17:22:09 +0200
Subject: [PATCH 22/22] Fix tests to work with cartesian index frequencies

---
 test/Cartesian.jl   | 10 +++++-----
 test/General.jl     | 14 +++++++-------
 test/MDFInMemory.jl | 10 +++++-----
 test/MDFv1.jl       |  6 +++---
 4 files changed, 20 insertions(+), 20 deletions(-)

diff --git a/test/Cartesian.jl b/test/Cartesian.jl
index 05e26cb9..95636ee4 100644
--- a/test/Cartesian.jl
+++ b/test/Cartesian.jl
@@ -7,27 +7,27 @@ fnSMFDSP = joinpath(tmpdir,"mdf","systemMatrixCartesianSP.mdf")
 smTD = MPIFile(fnSMTD)
 @test typeof(smTD) <: MDFFileV2
 
-
+freqs = collect(vec(CartesianIndices((10, 1))))
 @test acqNumPeriodsPerFrame(smTD) == 6500
 @test rxNumSamplingPoints(smTD) == 76
-@test size(getSystemMatrix(smTD,1:10)) == (81,6500*10)
+@test size(getSystemMatrix(smTD,freqs)) == (81,6500*10)
 @test size(getMeasurements(smTD)) == (76, 2, 6500, 81)
 @test size(getMeasurements(smTD, numPeriodAverages=65)) == (76, 2, 100, 81)
 @test size(getMeasurements(smTD, numPeriodAverages=65, numPeriodGrouping=100)) == (7600, 2, 1, 81)
-@test size(getMeasurementsFD(smTD, frequencies=1:10)) == (10, 6500, 81)
+@test size(getMeasurementsFD(smTD, frequencies=freqs)) == (10, 6500, 81)
 
 saveasMDF(fnSMFDMP, smTD, numPeriodAverages=65, applyCalibPostprocessing=true)
 smFDMP = MPIFile(fnSMFDMP)
 
 @test acqNumPeriodsPerFrame(smFDMP) == 100
 @test rxNumSamplingPoints(smFDMP) == 76
-@test size(getSystemMatrix(smFDMP,1:10)) == (81,10*100)
+@test size(getSystemMatrix(smFDMP,freqs)) == (81,10*100)
 
 saveasMDF(fnSMFDSP, smTD, numPeriodAverages=65, applyCalibPostprocessing=true, numPeriodGrouping=100)
 smFDSP = MPIFile(fnSMFDSP)
 
 @test acqNumPeriodsPerFrame(smFDSP) == 1
 @test rxNumSamplingPoints(smFDSP) == 7600
-@test size(getSystemMatrix(smFDSP,1:10)) == (81,10)
+@test size(getSystemMatrix(smFDSP,freqs)) == (81,10)
 
 end
\ No newline at end of file
diff --git a/test/General.jl b/test/General.jl
index 156adee5..9b24da47 100644
--- a/test/General.jl
+++ b/test/General.jl
@@ -110,7 +110,7 @@
 
       @test size(getMeasurementsFD(mdf, numAverages=10, frames=1:100, loadasreal=true)) == (1634,3,1,10)
 
-      @test size(getMeasurementsFD(mdf,frequencies=1:10, numAverages=10)) == (10,1,50)
+      @test size(getMeasurementsFD(mdf,frequencies=collect(vec(CartesianIndices((10, 1)))), numAverages=10)) == (10,1,50)
 
     end
   end
@@ -141,7 +141,7 @@
       @info "Test $sm"
 
       @test size( systemMatrixWithBG(sm) ) == (1959,817,3,1)
-      @test size( systemMatrix(sm,1:10) ) == (1936,10)
+      @test size( systemMatrix(sm,collect(vec(CartesianIndices((10, 1))))) ) == (1936,10)
 
       @test size(rxTransferFunction(sm)) == (817, 3)
       @test rxHasTransferFunction(sm) == true
@@ -164,9 +164,9 @@
       @test size(filterFrequencies(sm, SNRThresh = 5)) == (147,)
       #@test size(filterFrequencies(sm, numUsedFreqs = 100)) == (100,) # not working
 
-      @test size(getSystemMatrix(sm,1:10)) == (1936,10)
-      @test size(getSystemMatrix(sm,1:10,loadasreal=true)) == (1936,20)
-      @test size(getSystemMatrix(sm,1:10,bgCorrection=true)) == (1936,10)
+      @test size(getSystemMatrix(sm,collect(vec(CartesianIndices((10, 1)))))) == (1936,10)
+      @test size(getSystemMatrix(sm,collect(vec(CartesianIndices((10, 1)))),loadasreal=true)) == (1936,20)
+      @test size(getSystemMatrix(sm,collect(vec(CartesianIndices((10, 1)))),bgCorrection=true)) == (1936,10)
       # test on the data level if the conversion was successful
       SNRThresh = 2
       freq = filterFrequencies(smBruker,SNRThresh=SNRThresh)
@@ -174,7 +174,7 @@
       S = getSystemMatrix(sm,frequencies=freq)
 
       relativeDeviation = zeros(Float32,length(freq))
-      for f in 1:length(freq)
+      for f in 1:length(map(f->f[1], freq))
         relativeDeviation[f] = norm(SBruker[:,f]-S[:,f])/norm(SBruker[:,f])
       end
       # test if relative deviation for most of the frequency components is below 0.003
@@ -253,7 +253,7 @@
       @info "Test $sm"
 
       @test size( systemMatrixWithBG(sm) ) == (67,52,3,1)
-      @test size( systemMatrix(sm,1:10) ) == (60,10)
+      @test size( systemMatrix(sm,collect(vec(CartesianIndices((10, 1))))) ) == (60,10)
       @test size( systemMatrix(sm) ) == (60,52,3,1)
 
       @test size(calibSNR(sm)) == (52, 3, 1)
diff --git a/test/MDFInMemory.jl b/test/MDFInMemory.jl
index 8698cce8..e94f8d80 100644
--- a/test/MDFInMemory.jl
+++ b/test/MDFInMemory.jl
@@ -147,7 +147,7 @@ end
 
       @test size(getMeasurementsFD(mdf, numAverages=10, frames=1:100, loadasreal=true)) == (1634,3,1,10)
 
-      @test size(getMeasurementsFD(mdf,frequencies=1:10, numAverages=10)) == (10,1,50)
+      @test size(getMeasurementsFD(mdf,frequencies=collect(vec(CartesianIndices((10, 1)))), numAverages=10)) == (10,1,50)
     end
 
     @testset "Calibration" begin
@@ -168,7 +168,7 @@ end
       @test typeof(sm) <: MDFFileV2
   
       @test size( systemMatrixWithBG(sm) ) == (1959,817,3,1)
-      @test size( systemMatrix(sm,1:10) ) == (1936,10)
+      @test size( systemMatrix(sm,collect(vec(CartesianIndices((10, 1))))) ) == (1936,10)
 
       @test size(rxTransferFunction(sm)) == (817, 3)
       @test rxHasTransferFunction(sm) == true
@@ -198,9 +198,9 @@ end
       @test size(filterFrequencies(sm, SNRThresh = 5)) == (147,)
       #@test size(filterFrequencies(sm, numUsedFreqs = 100)) == (100,) # not working
 
-      @test size(getSystemMatrix(sm,1:10)) == (1936,10)
-      @test size(getSystemMatrix(sm,1:10,loadasreal=true)) == (1936,20)
-      @test size(getSystemMatrix(sm,1:10,bgCorrection=true)) == (1936,10)
+      @test size(getSystemMatrix(sm,collect(vec(CartesianIndices((10, 1)))))) == (1936,10)
+      @test size(getSystemMatrix(sm,collect(vec(CartesianIndices((10, 1)))),loadasreal=true)) == (1936,20)
+      @test size(getSystemMatrix(sm,collect(vec(CartesianIndices((10, 1)))),bgCorrection=true)) == (1936,10)
       # test on the data level if the conversion was successfull
       SNRThresh = 2
       freq = filterFrequencies(smBruker,SNRThresh=SNRThresh)
diff --git a/test/MDFv1.jl b/test/MDFv1.jl
index a5e275a2..9afcc348 100644
--- a/test/MDFv1.jl
+++ b/test/MDFv1.jl
@@ -96,7 +96,7 @@ for mdf in (mdfv1,mdfv2)
 
   @test size(getMeasurementsFD(mdf, numAverages=10, frames=1:500, loadasreal=true)) == (1634,3,1,50)
 
-  @test size(getMeasurementsFD(mdf,frequencies=1:10, numAverages=10)) == (10,1,50)
+  @test size(getMeasurementsFD(mdf,frequencies=collect(vec(CartesianIndices((10, 1)))), numAverages=10)) == (10,1,50)
 
 end
 
@@ -132,8 +132,8 @@ for sm in (smv1,smv2)
   @test size(filterFrequencies(sm, SNRThresh = 5)) == (147,)
   #@test size(filterFrequencies(sm, numUsedFreqs = 100)) == (100,) # not working
 
-  @test size(getSystemMatrix(sm,1:10)) == (1936,10)
-  @test size(getSystemMatrix(sm,1:10,loadasreal=true)) == (1936,20)
+  @test size(getSystemMatrix(sm,collect(vec(CartesianIndices((10, 1)))))) == (1936,10)
+  @test size(getSystemMatrix(sm,collect(vec(CartesianIndices((10, 1)))),loadasreal=true)) == (1936,20)
 end
 
 end