diff --git a/.github/workflows/Semgrep.yml b/.github/workflows/Semgrep.yml
new file mode 100644
index 00000000..62f63cfe
--- /dev/null
+++ b/.github/workflows/Semgrep.yml
@@ -0,0 +1,66 @@
+# Based on:
+# https://semgrep.dev/docs/semgrep-ci/sample-ci-configs#github-actions
+# https://0xdbe.github.io/GitHub-HowToEnableCodeScanningWithSemgrep/
+# https://medium.com/@mostafa.elnakeb/supercharging-your-code-quality-with-semgrep-sast-in-github-actions-c8f30eb26655
+# Name of this GitHub Actions workflow.
+name: Semgrep OSS scan
+
+on:
+  # Scan on-demand through GitHub Actions interface:
+  workflow_dispatch:
+    branches:
+      - master
+  # Schedule the CI job (this method uses cron syntax):
+  schedule:
+    - cron: '0 0 * * 1' # Run at start of week 
+    
+jobs:
+  semgrep:
+    # User definable name of this GitHub Actions job.
+    name: semgrep-oss/scan
+    # If you are self-hosting, change the following `runs-on` value: 
+    runs-on: ubuntu-latest
+
+    steps:
+        # Checkout the repository.
+      - name: Clone source code
+        uses: actions/checkout@v4
+
+        # Checkout custom rules
+      - name: Checkout custom rules
+        uses: actions/checkout@v4
+        with:
+          repository: JuliaComputing/semgrep-rules-julia
+          ref: main
+          path: ./JuliaRules
+
+        # Prepare Python
+      - uses: actions/setup-python@v5
+        with:
+          python-version: '3.10'
+
+        # Install Semgrep
+      - name: Install Semgrep
+        run: python3 -m pip install semgrep
+
+        # Run Semgrep
+      - name: Scan with Semgrep
+        run: | 
+            semgrep scan \
+                --config ./JuliaRules/rules \
+                --metrics=off \
+                --sarif --output report.sarif \
+                --oss-only \
+                --exclude=JuliaRules
+    
+      - name: Save Semgrep report
+        uses: actions/upload-artifact@v4
+        with:
+          name: report.sarif
+          path: report.sarif
+      
+      - name: Upload Semgrep report
+        uses: github/codeql-action/upload-sarif@v3
+        with:
+          sarif_file: report.sarif
+          category: semgrep
diff --git a/.github/workflows/breakage.yml b/.github/workflows/breakage.yml
index 495e05b4..a96586d7 100644
--- a/.github/workflows/breakage.yml
+++ b/.github/workflows/breakage.yml
@@ -1,15 +1,15 @@
 name: Breakage 
 # Based on: https://github.com/JuliaSmoothOptimizers/LinearOperators.jl/blob/main/.github/workflows/Breakage.yml
 on:
-    pull_request:
-        branches:
-          - master
-    push:
-        branches:
-          - master
-        tags: '*'
+  pull_request:
+    branches:
+      - master
+    types: [opened, reopened, labeled, synchronize]
+
+  
 jobs:
   break:
+    if: ${{  contains(github.event.pull_request.labels.*.name, 'breakage')  }}
     runs-on: ubuntu-latest
     strategy:
       fail-fast: false
diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
index dfff1bf7..78d8128b 100644
--- a/.github/workflows/ci.yml
+++ b/.github/workflows/ci.yml
@@ -43,7 +43,7 @@ jobs:
       - uses: julia-actions/julia-buildpkg@v1
       - uses: julia-actions/julia-runtest@v1
       - uses: julia-actions/julia-processcoverage@v1
-      - uses: codecov/codecov-action@v4
+      - uses: codecov/codecov-action@v5
         env:
           CODECOV_TOKEN: ${{ secrets.CODECOV_TOKEN }}
         with:
diff --git a/.gitignore b/.gitignore
index f698500e..8a6861d6 100644
--- a/.gitignore
+++ b/.gitignore
@@ -25,5 +25,7 @@ test/data/*
 # committed for packages, but should be committed for applications that require a static
 # environment.
 Manifest.toml
+# also ignore julia-version specific Manifests (>1.11)
+Manifest-v*.toml
 
 .vscode
\ No newline at end of file
diff --git a/Project.toml b/Project.toml
index 3c8cd589..38ba9250 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.16.0"
+version = "0.16.1"
 
 [deps]
 AxisArrays = "39de3d68-74b9-583c-8d2d-e117c070f3a9"
@@ -35,25 +35,41 @@ UnitfulAngles = "6fb2a4bd-7999-5318-a3b2-8ad61056cd98"
 UnitfulParsableString = "06c00241-927a-4d5b-bb5e-6b5a2ada3567"
 
 [compat]
-Aqua = "0.6"
+Aqua = "0.8"
 AxisArrays = "0.3, 0.4"
 CodecZlib = "0.7"
+Dates = "1"
 DelimitedFiles = "1"
 DocStringExtensions = "0.8, 0.9"
 FFTW = "1.3"
 FileIO = "1.0"
 Graphics = "0.4, 1.0"
-HDF5 = "0.14, 0.15, 0.16"
+HDF5 = "0.14, 0.15, 0.16, 0.17"
+HTTP = "1"
 ImageAxes = "0.6"
 ImageMetadata = "0.9"
 Inflate = "0.1.2"
-Interpolations = "0.12, 0.13, 0.14"
-LinearOperatorCollection = "1.1"
-StableRNGs = "1.0.0"
+InteractiveUtils = "1"
+Interpolations = "0.12, 0.13, 0.14, 0.15"
+LazyArtifacts = "1"
+LinearAlgebra = "1"
+LinearOperatorCollection = "1, 2"
+Mmap = "1"
+Pkg = "1"
+REPL = "1"
+Random = "1"
+Scratch = "1"
+SHA = "0.7"
+StableRNGs = "1"
+Statistics = "1"
+Test = "1"
+TOML = "1"
+Tar = "1"
+UUIDs = "1"
 Unitful = "1.13, 1.14, 1.15, 1.16, 1.17"
-UnitfulAngles = "0.6.1"
+UnitfulAngles = "0.6.1, 0.7"
 UnitfulParsableString = "0.1.6"
-julia = "1.5"
+julia = "1.9"
 
 [extras]
 Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
diff --git a/src/Brukerfile.jl b/src/Brukerfile.jl
index 33f0fd6c..7cdb04f1 100644
--- a/src/Brukerfile.jl
+++ b/src/Brukerfile.jl
@@ -186,7 +186,7 @@ studyNameOld(b::BrukerFile) = string(latin1toutf8(b["VisuSubjectId"])*latin1tout
                                   b["VisuStudyNumber"])
 studyNumber(b::BrukerFile) = parse(Int64,b["VisuStudyNumber"])
 function studyUuid(b::BrukerFile)
-  rng = MersenneTwister(hash(b["VisuStudyUid"])) # use VisuStudyUid as seed to generate uuid4
+  rng = StableRNG(hash(b["VisuStudyUid"])) # use VisuStudyUid as seed to generate uuid4
   return uuid4(rng)
 end
 studyDescription(b::BrukerFile) = "n.a."
@@ -207,7 +207,7 @@ function experimentNumber(b::BrukerFile)
 end
 
 function experimentUuid(b::BrukerFile)
-  rng = MersenneTwister(hash(b["VisuUid"])) # use VisuUid as seed to generate uuid4
+  rng = StableRNG(hash(b["VisuUid"])) # use VisuUid as seed to generate uuid4
   return uuid4(rng)
 end
 experimentDescription(b::BrukerFile) = latin1toutf8(b["ACQ_scan_name"])
diff --git a/src/Derived.jl b/src/Derived.jl
index 051b1197..9bf8a950 100644
--- a/src/Derived.jl
+++ b/src/Derived.jl
@@ -4,19 +4,12 @@ export acqNumFGFrames, acqNumBGFrames, acqOffsetFieldShift, acqFramePeriod,
        rxNumFrequencies, rxFrequencies, rxTimePoints,
        measFGFrameIdx, measBGFrameIdx, measBGFrameBlockLengths
 
-rxNumFrequencies(f::MPIFile, numPeriodGrouping=1) = floor(Int,rxNumSamplingPoints(f)*numPeriodGrouping ./ 2 .+ 1)
+rxNumFrequencies(f::MPIFile, numPeriodGrouping=1) = length(rfftfreq(rxNumSamplingPoints(f)*numPeriodGrouping))
 
-function rxFrequencies(f::MPIFile)
-  numFreq = rxNumFrequencies(f)
-  a = collect(0:(numFreq-1))./(numFreq-1).*rxBandwidth(f)
-  return a
-end
+rxFrequencies(f::MPIFile, numPeriodGrouping=1) = rfftfreq(rxNumSamplingPoints(f)*numPeriodGrouping, 2rxBandwidth(f)) |> collect
+
+rxTimePoints(f::MPIFile, numPeriodGrouping=1) = range(0, step=1/2rxBandwidth(f), length=rxNumSamplingPoints(f)*numPeriodGrouping) |> collect
 
-function rxTimePoints(f::MPIFile)
-  numTP = rxNumSamplingPoints(f)
-  a = collect(0:(numTP-1))./(numTP).*dfCycle(f)
-  return a
-end
 
 function acqGradientDiag(f::MPIFile)
   g = acqGradient(f)
@@ -44,36 +37,8 @@ end
 acqNumFGFrames(f::MPIFile) = acqNumFrames(f) - acqNumBGFrames(f)
 acqNumBGFrames(f::MPIFile) = sum(measIsBGFrame(f))
 
-function measBGFrameIdx(f::MPIFile)
-  idx = zeros(Int64, acqNumBGFrames(f))
-  j = 1
-  mask = measIsBGFrame(f)
-  for i=1:acqNumFrames(f)
-    if mask[i]
-      idx[j] = i
-      j += 1
-    end
-  end
-  return idx
-end
-
-function measFGFrameIdx(f::MPIFile)
-  mask = measIsBGFrame(f)
-  if !any(mask)
-    #shortcut
-    return 1:acqNumFrames(f)
-  end
-  idx = zeros(Int64, acqNumFGFrames(f))
-  j = 1
-  for i=1:acqNumFrames(f)
-    if !mask[i]
-      idx[j] = i
-      j += 1
-    end
-  end
-  return idx
-end
-
+measBGFrameIdx(f::MPIFile) = findall(measIsBGFrame(f))
+measFGFrameIdx(f::MPIFile) = findall(.!measIsBGFrame(f))
 
 function measBGFrameBlockLengths(mask)
   len = Vector{Int}(undef,0)
diff --git a/src/FrequencyFilter.jl b/src/FrequencyFilter.jl
index f5f60a1b..df2f3d3d 100644
--- a/src/FrequencyFilter.jl
+++ b/src/FrequencyFilter.jl
@@ -28,10 +28,13 @@ Supported keyword arguments:
 * stepsize
 * maxMixingOrder
 * sortByMixFactors
+* numPeriodGrouping
+* numSideBandFreqs
+* stopBands
 """
 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, numSidebandFreqs = -1)
+                maxMixingOrder=-1, numPeriodAverages=1, numPeriodGrouping=1, numSidebandFreqs = -1, stopBands = nothing)
 
   nFreq = rxNumFrequencies(f, numPeriodGrouping)
   nReceivers = rxNumChannels(f)
@@ -39,7 +42,7 @@ function filterFrequencies(f::MPIFile; SNRThresh=-1, minFreq=0,
   freqs = measIsFrequencySelection(f) ? measFrequencySelection(f) : 1:nFreq
 
   if numPeriodGrouping == 1
-    freqIndices = vec([CartesianIndex{2}(i, j) for i in freqs, j in recChannels])
+    freqIndices = vec([CartesianIndex{2}(i, j) for i in freqs, j in intersect(1:nReceivers, recChannels)])
   else
     freqIndices = collect(vec(CartesianIndices((nFreq, nReceivers))))
     filterFrequenciesBySelection!(freqIndices, f) 
@@ -88,6 +91,10 @@ function filterFrequencies(f::MPIFile; SNRThresh=-1, minFreq=0,
     filterFrequenciesByStepsize!(freqIndices, stepsize)
   end
 
+  if !isnothing(stopBands)
+    filterFrequenciesByStopBands!(freqIndices, f, stopBands; numPeriodGrouping = numPeriodGrouping)
+  end
+
   return collect(vec(freqIndices))
 end
 
@@ -104,19 +111,17 @@ filterFrequenciesByChannel!(indices, channels, sorted = issorted(channels)) = fi
 
 export filterFrequenciesByMinFreq!
 function filterFrequenciesByMinFreq!(indices, f::MPIFile, minFreq; numPeriodGrouping = 1)
-  nFreq = rxNumFrequencies(f, numPeriodGrouping)
-  minIdx = floor(Int, minFreq / rxBandwidth(f) * (nFreq-1) ) + 1
+  minIdx = searchsortedfirst(rfftfreq(rxNumSamplingPoints(f)*numPeriodGrouping, 2rxBandwidth(f)), minFreq)
   return filterFrequenciesByMinIdx!(indices, minIdx)
 end
-filterFrequenciesByMinIdx!(indices, minIdx) = minIdx > 0 ?  filter!(x -> x[1] > minIdx, indices) : indices 
+filterFrequenciesByMinIdx!(indices, minIdx) = minIdx > 0 ?  filter!(x -> x[1] >= minIdx, indices) : indices 
 
 export filterFrequenciesByMaxFreq!
 function filterFrequenciesByMaxFreq!(indices, f::MPIFile, maxFreq; numPeriodGrouping = 1)
-  nFreq = rxNumFrequencies(f, numPeriodGrouping)
-  maxIdx = ceil(Int, maxFreq / rxBandwidth(f) * (nFreq-1) ) + 1
+  maxIdx = searchsortedlast(rfftfreq(rxNumSamplingPoints(f)*numPeriodGrouping, 2rxBandwidth(f)), maxFreq)
   return filterFrequenciesByMaxIdx!(indices, maxIdx)
 end
-filterFrequenciesByMaxIdx!(indices, maxIdx) = filter!(x-> x[1] < maxIdx, indices)
+filterFrequenciesByMaxIdx!(indices, maxIdx) = filter!(x-> x[1] <= maxIdx, indices)
 
 export filterFrequenciesByMaxMixingOrder!
 filterFrequenciesByMaxMixingOrder!(indices, maxMixingOrder, f::MPIFile) = filterFrequenciesByMaxMixingOrder!(indices, maxMixingOrder, mixingFactors(f))
@@ -177,6 +182,29 @@ function filterFrequenciesByStepsize!(indices, stepsize)
   filter!(x -> insorted(x[1], stepIndices), indices)
 end
 
+export filterFrequenciesByStopBands!, filterFrequenciesByStopBand!
+function filterFrequenciesByStopBands!(indices, f::MPIFile, stopBands::Vector; kwargs...)
+  for stopBand in stopBands
+    filterFrequenciesByStopBand!(indices, f, stopBand; kwargs...)
+  end
+end
+filterFrequenciesByStopBands!(indices, f::MPIFile, stopBand::Union{Vector{Int64}, UnitRange, NTuple}; kwargs...) = filterFrequenciesByStopBands!(indices, f, [stopBand]; kwargs...)
+function filterFrequenciesByStopBand!(indices, f::MPIFile, stopBand::Vector{Int64}; numPeriodGrouping = 1)
+  if length(stopBand) != 2
+    error("Stop band are only defined for a start and stop value. Found $(length(stopBand)) values")
+  end
+  nFreq = rxNumFrequencies(f, numPeriodGrouping)
+  minIdx = searchsortedlast(rfftfreq(rxNumSamplingPoints(f)*numPeriodGrouping, 2rxBandwidth(f)), first(stopBand))
+  maxIdx = searchsortedlast(rfftfreq(rxNumSamplingPoints(f)*numPeriodGrouping, 2rxBandwidth(f)), last(stopBand))
+  return filterFrequenciesByStopBand!(indices, minIdx, maxIdx)
+end
+function filterFrequenciesByStopBand!(indices, f::MPIFile, stopBand::Union{UnitRange, NTuple{2, Int64}}; numPeriodGrouping = 1)
+  minIdx = searchsortedlast(rfftfreq(rxNumSamplingPoints(f)*numPeriodGrouping, 2rxBandwidth(f)), first(stopBand))
+  maxIdx = searchsortedlast(rfftfreq(rxNumSamplingPoints(f)*numPeriodGrouping, 2rxBandwidth(f)), last(stopBand))
+  return filterFrequenciesByStopBand!(indices, minIdx, maxIdx)
+end
+filterFrequenciesByStopBand!(indices, minIdx, maxIdx) = filter!(x-> x[1] < minIdx || x[1] > maxIdx, indices)
+
 function sortFrequencies(indices, f::MPIFile; numPeriodGrouping = 1, stepsize = 1, sortBySNR = false, sortByMixFactors = false)
   if sortBySNR && !sortByMixFactors
     indices = sortFrequenciesBySNR(indices, f, numPeriodGrouping = numPeriodGrouping, stepsize = stepsize)
diff --git a/src/MDF.jl b/src/MDF.jl
index 20051f38..9f965d8b 100644
--- a/src/MDF.jl
+++ b/src/MDF.jl
@@ -226,7 +226,7 @@ dfBaseFrequency(f::MDFFile)::Union{Float64, Missing} = @keyrequired f["/acquisit
 dfCustomWaveform(f::MDFFileV2)::Union{String, Nothing} = @keyoptional f["/acquisition/drivefield/customWaveform"]
 dfDivider(f::MDFFileV1) = @keyrequired addTrailingSingleton(
                 f["/acquisition/drivefield/divider"],2)
-dfDivider(f::MDFFileV2) = f["/acquisition/drivefield/divider"]
+dfDivider(f::MDFFileV2) = @keyrequired f["/acquisition/drivefield/divider"]
 dfWaveform(f::MDFFileV1) = "sine"
 function dfWaveform(f::MDFFileV2)::Union{Array{String, 2}, Missing}
   value = @keyrequired f["/acquisition/drivefield/waveform"]
@@ -401,7 +401,7 @@ measFramePermutation(f::MDFFileV1)::Union{Vector{Int64}, Nothing} = nothing
 measFramePermutation(f::MDFFileV2)::Union{Vector{Int64}, Nothing} = @keyoptional f["/measurement/framePermutation"]
 measSparsityTransformation(f::MDFFileV1)::Union{String, Nothing} = nothing
 measSparsityTransformation(f::MDFFileV2)::Union{String, Nothing} = @keyoptional f["/measurement/sparsityTransformation"]
-measSubsamplingIndices(f::MDFFileV2)::Union{Array{Integer, 4}, Nothing} = @keyoptional f["/measurement/subsamplingIndices"]
+measSubsamplingIndices(f::MDFFileV2)::Union{Array{Int64, 4}, Nothing} = @keyoptional f["/measurement/subsamplingIndices"]
 
 fullFramePermutation(f::MDFFile) = fullFramePermutation(f, calibIsMeanderingGrid(f))
 
diff --git a/src/MDFInMemory.jl b/src/MDFInMemory.jl
index bde3346e..401de0cf 100644
--- a/src/MDFInMemory.jl
+++ b/src/MDFInMemory.jl
@@ -418,7 +418,7 @@ mutable struct MDFv2Measurement <: MDFv2InMemoryPart
   "Name of the applied sparsity transformation; optional if !isSparsityTransformed"
   sparsityTransformation::Union{String, Nothing}
   "Subsampling indices \\beta{j,c,k,b}; optional if !isSparsityTransformed"
-  subsamplingIndices::Union{Array{Integer, 4}, Nothing}
+  subsamplingIndices::Union{Array{Int64, 4}, Nothing}
 
   function MDFv2Measurement(;
     data = missing,
@@ -1303,6 +1303,57 @@ for (fieldname, fieldtype) in zip(fieldnames(MDFv2InMemory), fieldtypes(MDFv2InM
   end
 end
 
+# Extract individual parts of an MDFFile as in-memory parts
+for (partFieldnameStr, partPrefix) ∈ prefixes
+  partFieldnameSymbol = Symbol(partFieldnameStr)
+  partInMemorySymbol = Symbol(lowercase(partFieldnameStr[6:end]))
+  type_ = getfield(MPIFiles, partFieldnameSymbol)
+  fields = fieldnames(type_)
+  functionNames = [let fieldString=string(field_); partPrefix != "" ? Symbol(partPrefix*uppercase(fieldString[1])*fieldString[2:end]) : Symbol(fieldString) end for field_ ∈ fields if field_ ∉ [:drivefield, :receiver]]
+
+  @eval begin
+    @doc $"""
+      $partFieldnameStr(file::MDFFile)
+
+    Create a `$partFieldnameStr` from the respective section in the given `file`.
+    """
+    function $partFieldnameSymbol(file::MDFFile)::$type_
+      instance = $partFieldnameSymbol()
+      for (fieldname_, functionName_) in zip($fields, $functionNames)
+        setfield!(instance, fieldname_, eval(functionName_)(file))
+      end
+
+      return instance
+    end
+  end
+
+  if !(partFieldnameStr == "MDFv2Drivefield" || partFieldnameStr == "MDFv2Receiver")
+    @eval begin
+      @doc $"""
+        $partFieldnameStr(file::MDFv2InMemory)
+
+      Create a `$partFieldnameStr` from the respective section in the given `file`.
+      """
+      function $partFieldnameSymbol(file::MDFv2InMemory)::$type_
+        return file.$partInMemorySymbol
+      end
+    end
+  else
+    @eval begin
+      @doc $"""
+        $partFieldnameStr(file::MDFv2InMemory)
+
+      Create a `$partFieldnameStr` from the respective section in the given `file`.
+      """
+      function $partFieldnameSymbol(file::MDFv2InMemory)::$type_
+        return file.acquisition.$partInMemorySymbol
+      end
+    end
+  end
+end
+MDFv2Drivefield(part::MDFv2Acquisition) = part.drivefield
+MDFv2Receiver(part::MDFv2Acquisition) = part.receiver
+
 # And some utility functions
 measIsCalibProcessed(mdf::MDFv2InMemory)::Union{Bool, Missing} = measIsFramePermutation(mdf) && 
                                                                  measIsFourierTransformed(mdf) &&
diff --git a/src/Positions/Positions.jl b/src/Positions/Positions.jl
index 36def6f6..9f95ee31 100644
--- a/src/Positions/Positions.jl
+++ b/src/Positions/Positions.jl
@@ -489,8 +489,8 @@ function getindex(rpos::UniformRandomPositions{AxisAlignedBox}, i::Integer)
     throw(BoundsError(rpos,i))
   else
     # make sure Positions are randomly generated from given seed
-    mersenneTwister = MersenneTwister(seed(rpos))
-    rP = rand(mersenneTwister, 3, i)[:,i]
+    rng = StableRNG(seed(rpos))
+    rP = rand(rng, 3, i)[:,i]
     return (rP.-0.5).*fieldOfView(rpos)+fieldOfViewCenter(rpos)
   end
 end
@@ -500,9 +500,9 @@ function getindex(rpos::UniformRandomPositions{Ball}, i::Integer)
     throw(BoundsError(rpos,i))
   else
     # make sure Positions are randomly generated from given seed
-    mersenneTwister = MersenneTwister(seed(rpos))
-    D = rand(mersenneTwister, i)[i]
-    P = randn(mersenneTwister, 3, i)[:,i]
+    rng = StableRNG(seed(rpos))
+    D = rand(rng, i)[i]
+    P = randn(rng, 3, i)[:,i]
     return radius(rpos)*D^(1/3)*normalize(P)+fieldOfViewCenter(rpos)
   end
 end
diff --git a/src/RecoData.jl b/src/RecoData.jl
index d568df86..db7a52e9 100644
--- a/src/RecoData.jl
+++ b/src/RecoData.jl
@@ -21,8 +21,8 @@ function saveRecoData(filename, image::ImageMeta)
 
   params = properties(image)
   params[:recoData] = c
-  params[:recoFov] = collect(grid) .* collect(converttometer(pixelspacing(image)))
-  params[:recoFovCenter] = collect(converttometer(imcenter(image)))[1:3]
+  params[:recoFov] = Float64.(collect(grid) .* collect(converttometer(pixelspacing(image))))
+  params[:recoFovCenter] = Float64.(collect(converttometer(imcenter(image)))[1:3])
   params[:recoSize] = collect(grid)
   params[:recoOrder] = "xyz"
   if haskey(params,:recoParams)
diff --git a/test/FrequencyFilter.jl b/test/FrequencyFilter.jl
new file mode 100644
index 00000000..f814a568
--- /dev/null
+++ b/test/FrequencyFilter.jl
@@ -0,0 +1,98 @@
+@testset "Frequency Filter" begin
+  mdf = MDFv2InMemory()
+  numChannels = 3
+  bandwidth = 256
+  numSamplingPoints = 2 * bandwidth
+  recv = MDFv2Receiver(;bandwidth = bandwidth, numSamplingPoints = numSamplingPoints, numChannels = numChannels)
+  mdf.acquisition = MDFv2Acquisition(;receiver = recv)
+  mdf.measurement = MDFv2Measurement(;isFrequencySelection = false)
+
+  unfiltered = filterFrequencies(mdf)
+  @test length(unfiltered) == (bandwidth + 1) * numChannels
+  @test eltype(unfiltered) == CartesianIndex{2}
+  nFreq = bandwidth + 1
+  frequencies = rxFrequencies(mdf)
+
+  @testset "Filter Receive Channels" begin
+    @test length(filterFrequencies(mdf, recChannels = 1:numChannels)) == length(unfiltered)
+    @test length(filterFrequencies(mdf, recChannels = 1:numChannels + 1)) == length(unfiltered)
+
+    freqs = filterFrequencies(mdf, recChannels = 1:numChannels - 1)
+    @test length(freqs) == length(unfiltered) - (nFreq)
+    @test all(i -> in(i[2], 1:(numChannels - 1)), freqs)
+
+
+    freqs = filterFrequencies(mdf, recChannels = 2:2)
+    @test length(freqs) == nFreq
+    @test all(i -> in(i[2], 2:2), freqs)
+
+    chSelection = [1, 3]
+    freqs = filterFrequencies(mdf, recChannels = chSelection)
+    @test length(freqs) == length(unfiltered) - (nFreq)
+    @test all(i -> in(i[2], chSelection), freqs)
+
+    @test length(filterFrequencies(mdf, recChannels = 0:0)) == 0
+    @test length(filterFrequencies(mdf, recChannels = numChannels+1:numChannels+2)) == 0
+  end
+
+  @testset "Min and Max Frequencies" begin
+    # minFreq = 0 -> Offset Frequency with index 1
+    @test length(filterFrequencies(mdf, minFreq = frequencies[1])) == length(unfiltered)
+    @test length(filterFrequencies(mdf, minFreq = 1.0)) == length(unfiltered) - numChannels
+    @test length(filterFrequencies(mdf, minFreq = frequencies[end])) == numChannels
+    @test length(filterFrequencies(mdf, minFreq = frequencies[end] + 1.0)) == 0
+
+    @test length(filterFrequencies(mdf, maxFreq = frequencies[end])) == length(unfiltered)
+    @test length(filterFrequencies(mdf, maxFreq = frequencies[end - 1])) == length(unfiltered) - numChannels
+    @test length(filterFrequencies(mdf, maxFreq = 0.0)) == numChannels
+    @test isempty(filterFrequencies(mdf, maxFreq = -1.0))
+  end
+
+  @testset "Stop Bands" begin
+    @test length(filterFrequencies(mdf, stopBands = [(0, 0)])) == length(unfiltered) - numChannels
+    @test length(filterFrequencies(mdf, stopBands = [[0, 0]])) == length(unfiltered) - numChannels
+    @test length(filterFrequencies(mdf, stopBands = [0:0])) == length(unfiltered) - numChannels
+
+    @test length(filterFrequencies(mdf, stopBands = (0, 0))) == length(unfiltered) - numChannels
+    @test length(filterFrequencies(mdf, stopBands = [0, 0])) == length(unfiltered) - numChannels
+    @test length(filterFrequencies(mdf, stopBands = 0:0)) == length(unfiltered) - numChannels
+
+    @test length(filterFrequencies(mdf, stopBands = 0:bandwidth)) == 0
+
+    freqs = filterFrequencies(mdf, stopBands = 5:10)
+    @test length(freqs) == length(unfiltered) - length(5:10) * numChannels
+    @test !any(i -> in(i[1], (5:10) .+ 1), freqs)
+
+    stopBands = [(5, 10), [15, 20], 18:22]
+    freqs = filterFrequencies(mdf, stopBands = stopBands)
+    @test length(freqs) == length(unfiltered) - (length(5:10) + length(15:22)) * numChannels
+    @test !any(i -> in(i[1], (5:10) .+ 1) || in(i[1], (15:22) .+ 1), freqs)
+  end
+
+
+  snr = fill(0.0, nFreq, numChannels, 1)
+  # Each freq has an SNR equal to its index
+  snr[:, :, 1] .= 1:nFreq
+  mdf.calibration = MDFv2Calibration(;snr = snr)
+
+  @testset "SNR Threshold" begin
+    @test length(filterFrequencies(mdf, SNRThresh = 0)) == length(unfiltered)
+    @test length(filterFrequencies(mdf, SNRThresh = -1)) == length(unfiltered)
+    @test length(filterFrequencies(mdf, SNRThresh = nFreq)) == numChannels
+    @test length(filterFrequencies(mdf, SNRThresh = nFreq + 1)) == 0
+
+    freqs = filterFrequencies(mdf, SNRThresh = floor(nFreq / 2))
+    @test all(i -> i[1] >= floor(nFreq / 2), freqs)
+  end
+
+  @testset "Num Used Frequencies" begin
+    @test length(filterFrequencies(mdf, numUsedFreqs = 1)) == numChannels skip = true
+    @test length(filterFrequencies(mdf, numUsedFreqs = nFreq)) == length(unfiltered) skip = true
+    @test length(filterFrequencies(mdf, numUsedFreqs = nFreq + 1)) == length(unfiltered) skip = true
+    @test length(filterFrequencies(mdf, numUsedFreqs = nFreq - 1)) == length(unfiltered) - numChannels skip = true
+
+    #freqs = filterFrequencies(mdf, numUsedFreqs = 1)
+    #@test all(i -> i[1] == nFreq, freqs) skip = true
+  end
+
+end
\ No newline at end of file
diff --git a/test/General.jl b/test/General.jl
index f6c72293..c059177d 100644
--- a/test/General.jl
+++ b/test/General.jl
@@ -328,3 +328,25 @@ end
 	# test new measurement parameters
 	@test [any(measIsBGFrame(mdf))] == measIsBGFrame(mdfavg)
 end
+
+@testset "Testing rxFrequencies" begin
+  f = MDFv2InMemory(Dict{Symbol,Any}(:rxNumSamplingPoints=>625, :rxBandwidth=>125e6/8/2))
+  @test rxFrequencies(f)[1] ≈ 0.0
+  @test rxFrequencies(f)[2] ≈ 25e3
+  @test rxNumFrequencies(f) == 313
+  @test rxNumFrequencies(f, 10) == 3126
+  @test rxFrequencies(f,10)[1] ≈ 0
+  @test rxFrequencies(f,10)[2] ≈ 2500
+  @test rxFrequencies(f,10)[11] ≈ 25e3
+  @test rxTimePoints(f,10)[1:625] ≈ rxTimePoints(f)[1:625]
+
+  f = MDFv2InMemory(Dict{Symbol,Any}(:rxNumSamplingPoints=>96, :rxBandwidth=>125e6/50/2))
+  @test rxFrequencies(f)[1] ≈ 0.0
+  @test rxFrequencies(f)[2] ≈ 26041.666666666666
+  @test rxNumFrequencies(f) == 49
+  @test rxNumFrequencies(f, 10) == 481
+  @test rxFrequencies(f,10)[1] ≈ 0
+  @test rxFrequencies(f,10)[2] ≈ 26041.666666666666/10
+  @test rxFrequencies(f,10)[11] ≈ 26041.666666666666
+  @test rxTimePoints(f,10)[1:96] ≈ rxTimePoints(f)[1:96]
+end
\ No newline at end of file
diff --git a/test/MDFInMemory.jl b/test/MDFInMemory.jl
index c2b4ec47..27769499 100644
--- a/test/MDFInMemory.jl
+++ b/test/MDFInMemory.jl
@@ -41,6 +41,7 @@ end
 
   fnMeasV2_converted = joinpath(tmpdir,"mdfim","measurement_V2_converted.mdf")
   fnSMV2_converted = joinpath(tmpdir,"mdfim","systemMatrix_V2_converted.mdf")
+  fnPart_converted = joinpath(tmpdir,"mdfim","imparts.mdf")
 
   @testset "Conversion" begin
     @testset "Fields" begin
@@ -396,7 +397,8 @@ end
       order = "xyz",
       positions = fill(0.0, (3, O)),
       size = [1, 1, 1],
-      snr = fill(0.0, (K, C, J))
+      snr = fill(0.0, (K, C, J)),
+      isMeanderingGrid = false # Not in specs for MDF v2!
     )
     mdf.reconstruction = MDFv2Reconstruction(
       data = fill(0, (S, P, Q)),
@@ -582,4 +584,148 @@ end
       test_mdf_replacement(mdf, :recoSize, [1, 1, 2], "The product of `size` with `[1, 1, 2]` must equal P.")
     end
   end
+
+  @testset "Part retrieval" begin
+    A = 1 # tracer materials/injections for multi-color MPI
+    N = 2 # acquired frames (N = O + E), same as a spatial position for calibration
+    E = 1 # acquired background frames (E = N − O)
+    O = 1 # acquired foreground frames (O = N − E)
+    B = 1 # coefficients stored after sparsity transformation (B \\le O)
+    J = 1 # periods within one frame
+    Y = 1 # partitions of each patch position
+    C = 1 # receive channels
+    D = 1 # drive-field channels
+    F = 1 # frequencies describing the drive-field waveform
+    V = 800 # points sampled at receiver during one drive-field period
+    W = V # sampling points containing processed data (W = V if no frequency selection or bandwidth reduction has been applied)
+    K = Int(V/2 + 1)# frequencies describing the processed data (K = V/2 + 1 if no frequency selection or bandwidth reduction has been applied)
+    Q = 1 # frames in the reconstructed MPI data set
+    P = 1 # voxels in the reconstructed MPI data set
+    S = 1 # channels in the reconstructed MPI data set
+
+    mdf = MDFv2InMemory()
+    mdf.root = MDFv2Root(
+      time=DateTime("2021-04-26T17:12:21.686"),
+      uuid=UUID("946a039e-48de-47ee-957e-a15af437e0be"),
+      version=VersionNumber("2.1.0")
+    )
+    mdf.study = MDFv2Study(
+      description = "n.a.",
+      name = "n.a.",
+      number = 1,
+      time = DateTime("2021-04-26T17:12:21.686"),
+      uuid = UUID("f2f3ae66-2fc3-49f7-91f9-71c1c2dc15e7")
+    )
+    mdf.experiment = MDFv2Experiment(
+      description = "n.a.",
+      isSimulation = true,
+      name = "n.a.",
+      number = 1,
+      subject = "n.a.",
+      uuid = UUID("3076d4bc-a2ef-46ef-8a99-dcd047379b8a")
+    )
+    mdf.tracer = MDFv2Tracer(
+      batch = fill("n.a.", A),
+      concentration = fill(1.0, A),
+      injectionTime = fill(DateTime("2021-04-26T17:12:21.686"), A),
+      name = fill("MyAwesomeTracer", A),
+      solute = fill("n.a.", A),
+      vendor = fill("Me", A),
+      volume = fill(1.0, A),
+    )
+    mdf.scanner = MDFv2Scanner(
+      boreSize = 0.3,
+      facility = "MyAwesomeInstitute",
+      manufacturer = "MeMyselfAndI",
+      name = "MyAwesomeScanner",
+      operator = "JustMe",
+      topology = "FFL" # Of course ;)
+    )
+    drivefield = MDFv2Drivefield(
+      baseFrequency = 40e3,
+      cycle = 1600/40e3,
+      divider = fill(1600, (F, D)),
+      numChannels = D,
+      phase = fill(0.0, (F, D, J)),
+      strength = fill(1.0, (F, D, J)),
+      waveform = fill("sine", (F, D)),
+    )
+    receiver = MDFv2Receiver(
+      bandwidth = 20e3,
+      dataConversionFactor = repeat([1/2^16 0]', outer=(1, C)),
+      inductionFactor = fill(1.0, C),
+      numChannels = C,
+      numSamplingPoints = V,
+      transferFunction = fill(1+0.5im, (K, C)),
+      unit = "V"
+    )
+    mdf.acquisition = MDFv2Acquisition(
+      gradient = fill(0.0, (3, 3, Y, J)),
+      numAverages = 1,
+      numFrames = N,
+      numPeriodsPerFrame = J,
+      offsetField = fill(0.0, (3, Y, J)),
+      startTime = DateTime("2021-04-26T17:12:21.686"),
+      drivefield = drivefield,
+      receiver = receiver
+    )
+    mdf.measurement = MDFv2Measurement(;
+      data = fill(0, (K, C, J, N)),
+      # framePermutation = fill(0, N),
+      # frequencySelection = collect(1:K),
+      isBackgroundCorrected = false,
+      isBackgroundFrame = vcat(fill(true, E), fill(false, O)),
+      isFastFrameAxis = false,
+      isFourierTransformed = false,
+      isFramePermutation = false,
+      isFrequencySelection = false,
+      isSparsityTransformed = false,
+      isSpectralLeakageCorrected = false,
+      isTransferFunctionCorrected = false,
+      # sparsityTransformation = "DCT-I",
+      # subsamplingIndices = fill(0, (B, K, C, J))
+    )
+    mdf.calibration = MDFv2Calibration(
+      deltaSampleSize = fill(0.001, 3),
+      fieldOfView = fill(0.2, 3),
+      fieldOfViewCenter = fill(0.0, 3),
+      method = "robot",
+      offsetFields = fill(0.0, (3, O)),
+      order = "xyz",
+      positions = fill(0.0, (3, O)),
+      size = [1, 1, 1],
+      snr = fill(0.0, (K, C, J)),
+      isMeanderingGrid = false
+    )
+    mdf.reconstruction = MDFv2Reconstruction(
+      data = fill(0, (S, P, Q)),
+      fieldOfView = fill(0.2, 3),
+      fieldOfViewCenter = fill(0.0, 3),
+      isOverscanRegion = fill(false, P),
+      order = "xyz",
+      positions = fill(0.0, (3, P)),
+      size = [1, 1, 1]
+    )
+
+    saveasMDF(fnPart_converted, mdf)
+    file = MPIFile(fnPart_converted)
+
+    for partFieldnameStr in keys(MPIFiles.prefixes)
+      partInMemorySymbol = Symbol(lowercase(partFieldnameStr[6:end]))
+
+      partim = eval(Symbol(partFieldnameStr))(mdf)
+      part = eval(Symbol(partFieldnameStr))(file)
+      
+      if !(partFieldnameStr == "MDFv2Drivefield" || partFieldnameStr == "MDFv2Receiver")
+        @test partim == getfield(mdf, partInMemorySymbol)
+        @test all([getfield(part, fieldname_) == getfield(getfield(mdf, partInMemorySymbol), fieldname_) for fieldname_ ∈ fieldnames(typeof(part)) if fieldname_ ∉ [:drivefield, :receiver]])
+      else
+        @test partim == getfield(mdf.acquisition, partInMemorySymbol)
+        @test all([getfield(part, fieldname_) == getfield(getfield(mdf.acquisition, partInMemorySymbol), fieldname_) for fieldname_ ∈ fieldnames(typeof(part))])
+      end 
+    end
+
+    @test MDFv2Drivefield(MDFv2Acquisition(mdf)) == drivefield
+    @test MDFv2Receiver(MDFv2Acquisition(mdf)) == receiver
+  end
 end
\ No newline at end of file
diff --git a/test/Positions.jl b/test/Positions.jl
index 275a2bd4..d559df6b 100644
--- a/test/Positions.jl
+++ b/test/Positions.jl
@@ -228,9 +228,9 @@ pospath = joinpath(tmpdir,"positions","Positions.h5")
   @test domain.fov == fov
   @test domain.center == ctr
   rP1 = UniformRandomPositions(N,seed,domain)
-  @test rP1[1] == [0.09954904813158394,-0.13791259323857274,-1.446939519855107]Unitful.mm
-  @test rP1[2] == [-0.9812009131891462,1.3767776289892044,1.4206979394110573]Unitful.mm
-  @test rP1[3] == [-0.5883911667396526,-0.9692742011014337,1.3707474722677764]Unitful.mm
+  @test rP1[1] == [0.24154458805558643, -0.9262755616254505, 1.413400404619357]Unitful.mm
+  @test rP1[2] == [0.7303493695629726, -0.9870943521080697, 0.6143278667437428]Unitful.mm
+  @test rP1[3] == [-0.17686922698492702, 0.911915746834733, 0.817151846349423]Unitful.mm
   @test_throws BoundsError rP1[0]
   @test_throws BoundsError rP1[4]
   h5open(pospath, "w") do file
@@ -249,7 +249,7 @@ pospath = joinpath(tmpdir,"positions","Positions.h5")
   @test domain.radius == radius
   @test domain.center == ctr
   rP3 = UniformRandomPositions(N,seed,domain)
-  @test rP3[1] == [-6.715713750009747,0.4103832286623821,-4.525933276650638]Unitful.mm
+  @test rP3[1] == [1.9129764588101597, 5.876973430472611, 5.6027641441895994]Unitful.mm
   @test_throws BoundsError rP3[0]
   @test_throws BoundsError rP3[4]
   h5open(pospath, "w") do file
diff --git a/test/runtests.jl b/test/runtests.jl
index 6080e3d6..8263ae33 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -34,6 +34,7 @@ include("MultiMPIFile.jl")
 include("Reco.jl")
 include("IMT.jl")
 include("TransferFunction.jl")
+include("FrequencyFilter.jl")
 include("CustomSFMeas.jl")
 include("MDFInMemory.jl")
 include("MagneticFieldMeasurement.jl")