Merge pull request #96 from MagneticParticleImaging/stopBands

This PR adds a stopBands kwarg to frequency filtering + Tests

src/FrequencyFilter.jl

@@ -28,18 +28,21 @@ 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)
   nPeriods = 1 #acqNumPeriodsPerFrame(f)
   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
 
@@ -108,15 +115,15 @@ function filterFrequenciesByMinFreq!(indices, f::MPIFile, minFreq; numPeriodGrou
   minIdx = floor(Int, minFreq / rxBandwidth(f) * (nFreq-1) ) + 1
   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
   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 +184,30 @@ 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 = floor(Int, first(stopBand) / rxBandwidth(f) * (nFreq-1) ) + 1
+  maxIdx = ceil(Int, last(stopBand) / rxBandwidth(f) * (nFreq-1) ) + 1
+  return filterFrequenciesByStopBand!(indices, minIdx, maxIdx)
+end
+function filterFrequenciesByStopBand!(indices, f::MPIFile, stopBand::Union{UnitRange, NTuple{2, Int64}}; numPeriodGrouping = 1)
+  nFreq = rxNumFrequencies(f, numPeriodGrouping)
+  minIdx = floor(Int, first(stopBand) / rxBandwidth(f) * (nFreq-1) ) + 1
+  maxIdx = ceil(Int, last(stopBand) / rxBandwidth(f) * (nFreq-1) ) + 1
+  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)

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

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")