Updated interface of TransferFunction

Project.toml

@@ -29,16 +29,16 @@ StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 TOML = "fa267f1f-6049-4f14-aa54-33bafae1ed76"
 Tar = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e"
-Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 UUIDs = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"
 Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 UnitfulAngles = "6fb2a4bd-7999-5318-a3b2-8ad61056cd98"
+UnitfulParsableString = "06c00241-927a-4d5b-bb5e-6b5a2ada3567"
 
 [compat]
 Aqua = "0.6"
 AxisArrays = "0.3, 0.4"
 CodecZlib = "0.7"
-DelimitedFiles = "1.9.0"
+DelimitedFiles = "1"
 DocStringExtensions = "0.8, 0.9"
 FFTW = "1.3"
 FileIO = "1.0"
@@ -52,6 +52,7 @@ LinearOperatorCollection = "1.1"
 StableRNGs = "1.0.0"
 Unitful = "1.13, 1.14, 1.15, 1.16, 1.17"
 UnitfulAngles = "0.6.1"
+UnitfulParsableString = "0.1.6"
 julia = "1.5"
 
 [extras]

docs/make.jl

@@ -13,7 +13,8 @@ makedocs(
  "Measurements" => "measurements.md",
  "System Matrices" => "systemmatrix.md",
  "Frequency Filter" => "frequencyFilter.md",
- "Reconstruction Results" => "reconstruction.md"
+ "Reconstruction Results" => "reconstruction.md",
+ "Transfer Functions" => "transferfunction.md"
  # "Positions" => "positions.md"
  ],
  warnonly = [:missing_docs]

docs/src/transferfunction.md

@@ -0,0 +1,69 @@
+# Transfer functions
+
+MPIFiles defines the type `TransferFunction` which represents the system-properties of a linear time-invariant system in frequency space ([see also](https://en.wikipedia.org/wiki/Transfer_function)). In the context of MPIFiles this usually includes the properties of the receive chain of an MPI system, but can also hold information about other parts of the signal chain. 
+
+## Basic construction and data access
+The `TransferFunction` object is constructed from samples in frequency space and offers two ways to access the underlying data. The first uses indexing to directly access the underlying data:
+
+```@setup tf
+using MPIFiles
+```
+```@repl tf
+ f = collect(range(0,1e6,step=1e3));
+ # TransferFunction of a simple lowpass filter
+ tf = TransferFunction(f, 1 ./ (1 .+ im*f/1e4 )) 
+ tf[1] # Value of tf at first frequency sample (f = 0 Hz)
+```
+
+The second method has to be called on the object itself and uses linear interpolation to access the tf at any frequency:
+```@repl tf
+ tf(0) # Value of tf at f = 0 Hz
+ tf(1e4) # Value of tf at f = 10000 Hz
+```
+
+A `TransferFunction` can have multiple channels, which adds a second index to both data access functions. Directly accessing multiple channels is also possible. The complex array used to construct the `TransferFunction` needs to have the shape [number of frequency samples, channels].
+
+```@repl tf
+tf = TransferFunction(f, [1 ./(1 .+im*f/1e4) 1 ./(1 .+im*f/1e3)])
+tf[11,1]
+tf[11,2]
+tf(1e4,1)
+tf(1e4,2)
+tf(1e4, [1,2])
+tf(1e4,:)
+```
+## Units 
+To attach units to the `TransferFunction` the keyword-parameter `units` can to be used to give a `Unitful` unit to every channel of the tf. This can be useful if the transfer function is not dimensionless but relates two physical quantities, e.g. voltage and current in the form of an impedance. All **interpolated** accesses to tf data then return a `Unitful.Quantity`.
+
+```@repl tf
+R = 1; # Ohm
+L = 10e-6; # Henry
+RL = TransferFunction(f, R .+ im*2pi*f*L, units=["V/A"])
+RL([0,100e3])
+```
+
+## Saving and loading
+A `TransferFunction` object can be saved to and loaded from a .h5 file.
+
+```@docs
+MPIFiles.save(filename::String, tf::TransferFunction)
+MPIFiles.TransferFunction(::String)
+```
+
+## Additional constructors
+
+In addition to the constructor taking a single (complex) array it is also possible to give two arrays representing amplitude and phase.
+
+It is also possible to construct a `TransferFunction` from the transfer function data included in an `MPIFile`.
+
+```@docs
+ MPIFiles.TransferFunction(freq::Vector{<:Real}, ampdata::Array{<:Real,N}, phasedata::Array{<:Real,N}) where N
+ MPIFiles.TransferFunction(::MPIFile)
+```
+
+## Other interesting functions
+```@docs
+ TransferFunction
+ combine(::TransferFunction, ::TransferFunction)
+```
+

src/MDF.jl

@@ -144,7 +144,7 @@ tracerSolute(f::MDFFileV1)::Union{Vector{String}, Missing} = ["Fe"]
 function tracerInjectionTime(f::MDFFile)::Union{Vector{DateTime}, Nothing}
  p = typeof(f) <: MDFFileV1 ? "/tracer/time" : "/tracer/injectionTime"
  time = @keyoptional f[p]
- if isnothing(nothing)
+ if isnothing(time)
  return nothing
  end
 

src/MPIFiles.jl

@@ -26,6 +26,7 @@ using Pkg.Artifacts
 using Unitful
 using InteractiveUtils
 using UnitfulAngles
+using UnitfulParsableString
 using Inflate, SHA
 using StableRNGs
 using REPL: fielddoc

src/TransferFunction.jl

@@ -1,84 +1,168 @@
-export TransferFunction, sampleTF, setTF
+export TransferFunction, sampleTF, setTF, combine
 
+"""
+$(TYPEDEF)
+$(TYPEDFIELDS)
+
+
+ TransferFunction(freq_::Vector{<:Real}, datain::Array{<:Complex}; inductionFactor::Vector{<:Real}=ones(size(datain, 2)), units::Vector=Unitful.FreeUnits[Unitful.NoUnits for i in 1:size(datain, 2)])
+
+Create a `TransferFunction` from a complex data array at frequencies `freq_`.
+
+# Optional Keyword-Arguments:
+- `inductionFactor::Vector{<:Real}`: induction factor for each channel
+- `units::Vector`: units for each channel, can be either Unitful.FreeUnits or a string that can be parsed as a Unitful unit 
+
+"""
 mutable struct TransferFunction
  freq::Vector{Float64}
  data::Matrix{ComplexF64}
+ interpolator::Vector{AbstractInterpolation}
  inductionFactor::Vector{Float64}
+ units::Vector{Unitful.FreeUnits}
+
+
+ function TransferFunction(freq_::Vector{<:Real}, datain::Array{<:Number}; inductionFactor::Vector{<:Real}=ones(size(datain, 2)), units::Vector=Unitful.FreeUnits[Unitful.NoUnits for i in 1:size(datain, 2)])
+ parsed_units = Unitful.FreeUnits[]
+ for tmp in units
+ if isa(tmp, String); tmp = uparse(tmp) end # get correct unit from unit strings
+ if !isa(tmp, Unitful.FreeUnits); tmp = Unitful.unit(tmp) end # get correct unit for numbers and quantities, e.g. unit(1) or unit(1u"V/T"), includes the case that the string parsing returned one of these types
+ push!(parsed_units, tmp)
+ end
+
+ if length(freq_) != size(datain,1); error("Length of frequency axis ($(length(freq_))) does not match the number of samples in the data ($(size(datain,1)))!") end
+ if size(datain, 2) != length(inductionFactor); error("Number of channels in data ($(size(datain, 2))) does not match the number of channels in the given inductionFactor ($(size(inductionFactor,1)))!") end
+ if size(datain, 2) != length(units); error("Number of channels in data ($(size(datain, 2))) does not match the number of channels in the given units ($(size(units,1)))!") end
 
- function TransferFunction(freq_, datain::Array{T}, inductionFactor=ones(size(datain, 2))) where {T<:Complex}
- freq = freq_[1]:(freq_[2]-freq_[1]):freq_[end]
  data = reshape(deepcopy(datain), size(datain,1), size(datain, 2))
- return new(freq_, data, inductionFactor)
+ interpolator = [extrapolate(interpolate((freq_,), channel, Gridded(Linear())), Interpolations.Flat()) for channel in eachcol(data)]
+ return new(freq_, data, interpolator, inductionFactor, parsed_units)
  end
 end
 
-function TransferFunction(freq_, ampdata, phasedata, args...)
+Base.show(io::IO, ::MIME"text/plain", tf::TransferFunction) = print(io, "MPIFiles.TransferFunction: \n\t$(size(tf.data,2)) channel(s), units of $(string.(tf.units))\n\t$(size(tf.data,1)) frequency samples from $(tf.freq[1]) Hz to $(tf.freq[end]) Hz")
+
+"""
+$(TYPEDSIGNATURES)
+
+Create a `TransferFunction` from separate amplitude and phase arrays at frequencies `freq`.
+
+`ampdata` and `phasedata` should have the following shape: [frequencies, channels]
+"""
+function TransferFunction(freq::Vector{<:Real}, ampdata::Array{<:Real,N}, phasedata::Array{<:Real,N}; kwargs...) where N
+ if size(ampdata) != size(phasedata); error("The size of ampdata and phasedata must match!") end
  data = ampdata.*exp.(im.*phasedata)
- return TransferFunction(freq_, data, args...)
+ return TransferFunction(freq, data; kwargs...)
 end
 
+"""
+$(TYPEDSIGNATURES)
+
+Create a `TransferFunction` from a data file at `filename`.
+
+The file can be either a h5-File created with this package. Keyword arguments will be passed to `load_tf_fromVNA` 
+"""
 function TransferFunction(filename::String; kargs...)
  filenamebase, ext = splitext(filename)
- if  ext == ".h5"
+ if ext == ".h5"
  tf = load_tf(filename)
  else #if ext == "s1p" || ext == "s2p"
  tf = load_tf_fromVNA(filename; kargs...)
  end
  return tf
 end
 
+"""
+$(TYPEDSIGNATURES)
+
+Create a `TransferFunction` from the tf data saved in a MPIFile (see `rxTransferFunction`)
+"""
 function TransferFunction(file::MPIFile)
  tf_file = rxTransferFunction(file)
  inductionFactor = rxInductionFactor(file)
  f = collect(rfftfreq(rxNumSamplingPoints(file), rxBandwidth(file)*2))
  return TransferFunction(f, abs.(tf_file), angle.(tf_file), inductionFactor)
 end
 
-function getindex(tmf::TransferFunction, x::UnitRange, chan::Integer=1)
- a = tmf.data[x,chan]
- return a
-end
+"""
+ tf[i,j]
 
-#function getindex(tmf::TransferFunction, x::Real, chan::Integer=1)
-# a = tmf.interp[chan](x)
-# return a
-#end
+Directly access the underlying data of a `TransferFunction`
+"""
+function getindex(tf::TransferFunction, args...)
+ try 
+ return getindex(tf.data, args...) 
+ catch e
+ @warn "The indexing using square brackets on TransferFunction objects now always operates on the integer indizes of the underlying transfer function data. To use frequency interpolation, use tf(freq, channel) instead of tf[[freq],channel]."
+ rethrow(e)
+ end
+end
 
-function getindex(tmf::TransferFunction, X::Union{Vector,AbstractRange}, chan::Integer=1)
- I = extrapolate(interpolate((tmf.freq,), tmf.data[:,chan], Gridded(Linear())), Interpolations.Flat())
+"""
+  tf(f, chan::Integer)
 
- return [I(x) for x in X]
+Interpolated access to a `TransferFunction` at frequencies `f` and single channel `chan`
+"""
+function (tf::TransferFunction)(f, chan::Integer=1)
+ if chan>length(tf.interpolator); error("The TransferFunction only has $(length(tf.interpolator)) channel(s), unable to access channel $(chan)") end
+ return tf.interpolator[chan](f) .* tf.units[chan]
 end
 
-function getindex(tmf::TransferFunction, X::Union{Vector,AbstractRange}, chan::AbstractRange)
- out = zeros(ComplexF64, length(X), length(chan))
- for d=1:length(chan)
- out[:,d] = tmf[X,d]
- end
- return out
-end
+"""
+ tf(f, chan::AbstractVector{<:Integer})
+
+Interpolated access to a `TransferFunction` at frequencies `f` and channels `chan`
+"""
+(tf::TransferFunction)(f, chan::AbstractVector{<:Integer}) = hcat([tf(f,c) for c in chan]...)
+
+(tf::TransferFunction)(f, ::Colon) = tf(f, axes(tf.data,2))
 
 function load_tf(filename::String)
  tf = h5read(filename,"/transferFunction")
  freq = h5read(filename,"/frequencies")
  inductionFactor = h5read(filename,"/inductionFactor")
- return TransferFunction(freq,tf,inductionFactor)
+ unit = []
+ try 
+ unit = h5read(filename, "/unit")
+ catch # if h5read fails, it should mean that there is no units, maybe do something better here
+ return TransferFunction(freq, tf, inductionFactor=inductionFactor)
+ end
+ return TransferFunction(freq, tf, inductionFactor=inductionFactor, units=uparse.(unit))
 end
 
-function combine(tf1, tf2)
- freq = tf1.freq
- data = cat(tf1.data,tf2.data, dims=2)
- inductionFactor = cat(tf1.inductionFactor, tf2.inductionFactor, dims=1)
- return TransferFunction(freq, data, inductionFactor)
+"""
+$(TYPEDSIGNATURES)
+
+Combine two `TransferFunctions` along their channel dimension. If interpolate=false, will only work if the frequency samples are identical.
+"""
+function combine(tf1::TransferFunction, tf2::TransferFunction; interpolate=false)
+ if !interpolate
+ if tf1.freq != tf2.freq; error("The frequency axes of the transfer functions do not match. Can not combine!") end
+ freq = tf1.freq
+ data = cat(tf1.data,tf2.data, dims=2)
+ else
+ freq = unique(sort(cat(tf1.freq, tf2.freq, dims=1)))
+ data = cat(tf1(freq, :), tf2(freq, :), dims=2)
+ end
+ inductionFactor = cat(tf1.inductionFactor, tf2.inductionFactor, dims=1)
+ units = cat(tf1.units, tf2.units, dims=1)
+ return TransferFunction(freq, data, inductionFactor=inductionFactor, units=units)
 end
 
+"""
+$(TYPEDSIGNATURES)
+
+Save `tf` as a h5 file to `filename`
+"""
 function save(filename::String, tf::TransferFunction)
  h5write(filename, "/transferFunction", tf.data)
  h5write(filename, "/frequencies", tf.freq)
  h5write(filename, "/inductionFactor", tf.inductionFactor)
+ h5write(filename, "/unit", string.(tf.units))
  return nothing
 end
 
+"""TODO: fix function and write docs"""
 function load_tf_fromVNA(filename::String;
  frequencyWeighting=false,
  R = 50.0, #Ω
@@ -130,11 +214,10 @@
 end
 
 
-function sampleTF(tmf::TransferFunction, f::MPIFile)
+function sampleTF(tf::TransferFunction, f::MPIFile)
  freq = rxFrequencies(f)
  numChan = rxNumChannels(f)
- numFreq = length(freq)
- return tmf[freq,1:numChan]
+ return tf(freq,1:numChan)
 end
 
 

test/TransferFunction.jl

@@ -15,15 +15,51 @@ tfh5path = joinpath(tmpdir,"transferFunction","example.h5")
 
  @test a.freq == b.freq
  @test a.data == b.data
- @test a[[1],1] == a[[0.0],1]
+ @test a[1,1] == a(0.0,1)
 
  c = MPIFiles.combine(MPIFiles.combine(a,a),a)
  rm(tfh5path)
  MPIFiles.save(tfh5path, c)
- @test a[[1],1] == c[[1],2]
- @test a[[1],1] == c[[1],3]
+ @test a[1,1] == c[1,2]
+ @test a[1,1] == c[1,3]
 
  measBruker = MPIFile(fnMeasBruker)
  tf = sampleTF(c, measBruker)
  @test size(tf) == (rxNumFrequencies(measBruker), rxNumChannels(measBruker))
+
+ f = collect(range(0,1e6,step=1e3));
+ data = 1 ./ (1 .+ im*f/1e4 );
+ tf = TransferFunction(f, data) # TransferFunction of a simple lowpass filter
+ @test tf[1] == 1.0
+ @test tf[1:6,1] == data[1:6]
+ @test tf(0) == tf(0,1)
+
+ data = [1 ./(1 .+im*f/1e4) 1 ./(1 .+im*f/1e3)]
+ tf = TransferFunction(f, data)
+ @test tf[11,1] == data[11,1]
+ @test tf[11,2] == data[11,2]
+ @test tf[11,:] == data[11,:]
+ @test tf(1e4, [1,2]) == [data[11,1] 1 ./(1 .+im*1e4/1e3)]
+ @test tf(1e4,:) == tf(1e4, [1,2])
+ @test tf(1e4,:) == tf(1e4, 1:2)
+ @test tf(0) == tf(0,1)
+
+ tf = TransferFunction(f, data, units=["V/V", "A/V"])
+ @test tf(0) == 1u"V/V"
+ @test tf(0,2) == 1u"A/V"
+
+
+ f1 = collect(range(0,1e6,step=1e3));
+ data1 = 1 ./ (1 .+ im*f1/1e4 );
+
+ f2 = collect(range(0,1e7,step=1e4));
+ data2 = 1 ./ (1 .+ im*f2/1e4 );
+
+ tf1 = TransferFunction(f1, data1)
+ tf2 = TransferFunction(f2, data2)
+
+ @test_throws ErrorException tf_c = combine(tf1,tf2)
+ tf_c = combine(tf1,tf2,interpolate=true)
+ @test tf_c(101.5e3, :) ≈ [tf1(101.5e3) tf2(101.5e3)] atol=1e-12
+
 end