diff --git a/docs/src/transferfunction.md b/docs/src/transferfunction.md index 0c9dc5b2..24fa27c5 100644 --- a/docs/src/transferfunction.md +++ b/docs/src/transferfunction.md @@ -7,6 +7,7 @@ The `TransferFunction` object is constructed from samples in frequency space and ```@setup tf using MPIFiles +using MPIFiles.Unitful ``` ```@repl tf f = collect(range(0,1e6,step=1e3)); @@ -33,7 +34,9 @@ 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`. +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. Alternatively `data` can just be a Unitful.Quantity. Then `units` is ignored. + +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 @@ -42,6 +45,14 @@ RL = TransferFunction(f, R .+ im*2pi*f*L, units=["V/A"]) RL([0,100e3]) ``` +```@repl tf +f_unitful = collect(range(0u"Hz",1u"MHz",step=1u"kHz")); +R = 1u"Ω"; +L = 10u"µH"; +RL = TransferFunction(f_unitful, R .+ im*2pi*f_unitful*L .|> u"V/A") +RL([0,100e3]) +``` + ## Saving and loading A `TransferFunction` object can be saved to and loaded from a .h5 file. @@ -50,20 +61,20 @@ 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. +## Constructors +The `TransferFunction` constructor can take either a complex data array or two arrays representing the amplitude and phase of the transfer function. Unitful conversion is automatically done for all parameters. 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 + TransferFunction MPIFiles.TransferFunction(::MPIFile) ``` ## Other interesting functions ```@docs - TransferFunction combine(::TransferFunction, ::TransferFunction) + MPIFiles.load_tf_fromVNA + MPIFiles.processRxTransferFunction ``` diff --git a/src/MDF.jl b/src/MDF.jl index 9f965d8b..419fc8a8 100644 --- a/src/MDF.jl +++ b/src/MDF.jl @@ -252,17 +252,9 @@ function rxTransferFunction(f::MDFFile) return nothing end end -function rxTransferFunctionFileName(f::MDFFile) - parameter = "/acquisition/receiver/transferFunctionFileName" - if haskey(f.file, parameter) - return f[parameter] - else - return nothing - end -end -function rxHasTransferFunction(f::MDFFile) - haskey(f.file, "/acquisition/receiver/transferFunction") -end +rxTransferFunctionFileName(f::MDFFile) = @keyoptional f["/acquisition/receiver/transferFunctionFileName"] +rxHasTransferFunction(f::MDFFile) = haskey(f.file, "/acquisition/receiver/transferFunction") + rxInductionFactor(f::MDFFileV1) = nothing rxInductionFactor(f::MDFFileV2) = @keyoptional f["/acquisition/receiver/inductionFactor"] diff --git a/src/Measurements.jl b/src/Measurements.jl index 19fcc84e..e6aeda14 100644 --- a/src/Measurements.jl +++ b/src/Measurements.jl @@ -271,6 +271,9 @@ function getMeasurements(f::MPIFile, neglectBGFrames=true; if tfCorrection && !measIsTFCorrected(f) tf = rxTransferFunction(f) inductionFactor = rxInductionFactor(f) + if isnothing(tf) + error("No transfer function available in file, please use tfCorrection=false") + end if isnothing(inductionFactor) @warn "The file is missing the induction factor. The induction factor will be set to 1." inductionFactor = ones(Float64, rxNumChannels(f)) @@ -278,11 +281,8 @@ function getMeasurements(f::MPIFile, neglectBGFrames=true; J = size(data,1) dataF = rfft(data, 1) - dataF[2:end,:,:,:] ./= tf[2:end,:,:,:] - - if all(tf[1,:,:,:] .!= 0) && !any(isnan.(tf[1,:,:,:])) - dataF[1,:,:,:] ./= tf[1,:,:,:] - end + dataF ./= tf + map!(x -> isnan(x) ? zero(eltype(dataF)) : x, dataF, dataF) @warn "This measurement has been corrected with a Transfer Function. Name of TF: $(rxTransferFunctionFileName(f))" if !isnothing(inductionFactor) @@ -330,16 +330,15 @@ function getMeasurementsFD(f::MPIFile, args...; if tfCorrection && !measIsTFCorrected(f) tf = rxTransferFunction(f) inductionFactor = rxInductionFactor(f) + if isnothing(tf) + error("No transfer function available in file, please use tfCorrection=false") + end if isnothing(inductionFactor) @warn "The file is missing the induction factor. The induction factor will be set to 1." inductionFactor = ones(Float64, rxNumChannels(f)) end - data[2:end,:,:,:] ./= tf[2:end,:,:,:] - - if all(tf[1,:,:,:] .!= 0) && !any(isnan.(tf[1,:,:,:])) - data[1,:,:,:] ./= tf[1,:,:,:] - end + data ./= tf @warn "This measurement has been corrected with a Transfer Function. Name of TF: $(rxTransferFunctionFileName(f))" if !isnothing(inductionFactor) diff --git a/src/TransferFunction.jl b/src/TransferFunction.jl index 1961206e..45577a90 100644 --- a/src/TransferFunction.jl +++ b/src/TransferFunction.jl @@ -5,13 +5,15 @@ $(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)]) + TransferFunction(freq::Vector, data::Array, [phasedata]; [inductionFactor], [units]) -Create a `TransferFunction` from a complex data array at frequencies `freq_`. +Create a `TransferFunction` from a data array `data` at frequencies `freq`. +`freq` is given in Hz or should have a Unitful frequency unit attached +`data` should have the shape [frequencies, channels]. `data` can be either complex or real, if `data` is real a second array `phasedata` can be passed representing the phase in radians. Both the amplitude and the phase can have Unitful units. # 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 +- `units::Vector`: units for each channel, can be either Unitful.FreeUnits or a string that can be parsed as a Unitful unit. Instead of using this keyword, `data` can also have Unitful units attached, then the units keyword-argument is ignored. """ mutable struct TransferFunction @@ -22,7 +24,7 @@ mutable struct TransferFunction 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)]) + function 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)]) parsed_units = Unitful.FreeUnits[] for tmp in units if isa(tmp, String); tmp = uparse(tmp) end # get correct unit from unit strings @@ -42,25 +44,38 @@ end 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 +function TransferFunction(freq::Vector{<:Real}, ampdata::Array{<:Union{Real,Unitful.AbstractQuantity{<: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; kwargs...) end +TransferFunction(freq::Vector{<:Unitful.Frequency}, args...; kwargs...) = TransferFunction(ustrip.(u"Hz", freq), args...; kwargs...) +TransferFunction(freq::Vector{<:Real}, ampdata::Array{<:Union{Real,Unitful.AbstractQuantity{<:Real}},N}, phasedata::Array{<:Unitful.DimensionlessQuantity{<:Real},N}; kwargs...) where N = TransferFunction(freq, ampdata, ustrip.(u"rad", phasedata); kwargs...) +TransferFunction(freq::Vector{<:Real}, ampdata::Array{<:Union{Real,Unitful.AbstractQuantity{<:Real}},N}; kwargs...) where N = TransferFunction(freq, ampdata, zeros(size(ampdata)); kwargs...) + +function TransferFunction(freq::Vector{<:Real}, datain::Array{<:Unitful.AbstractQuantity{<:Complex},N}; units=nothing, kwargs...) where N + if !isnothing(units) + @warn "You passed explicit units combined with a Unitful data array. The explicit units will be ignored!" + end + units = Unitful.FreeUnits[] + for ch in eachcol(datain) + ch_units = Unitful.unit.(ch) + if all(y->y==ch_units[1], ch_units) + push!(units, ch_units[1]) + else + error("One TransferFunction channel must have identical units!") + end + end + return TransferFunction(freq, ustrip.(datain); units=units, 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` +The file can be either a h5-File created with this package or a file that is written by a VNA. Keyword arguments will be passed to `load_tf_fromVNA` """ function TransferFunction(filename::String; kargs...) filenamebase, ext = splitext(filename) @@ -80,8 +95,12 @@ Create a `TransferFunction` from the tf data saved in a MPIFile (see `rxTransfer 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=inductionFactor) + f = rxFrequencies(file) + if isnothing(inductionFactor) + return TransferFunction(f, abs.(tf_file), angle.(tf_file)) + else + return TransferFunction(f, abs.(tf_file), angle.(tf_file), inductionFactor=inductionFactor) + end end """ @@ -132,7 +151,6 @@ end """ $(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) @@ -142,16 +160,18 @@ function combine(tf1::TransferFunction, tf2::TransferFunction; interpolate=false 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) + data = cat(ustrip.(tf1(freq, :)), ustrip.(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 +combine(tf::TransferFunction, ::Nothing; interpolate=false) = combine(tf, TransferFunction(tf.freq, zeros(length(tf.freq)))) +combine(::Nothing, tf::TransferFunction; interpolate=false) = combine(TransferFunction(tf.freq, zeros(length(tf.freq))), tf) + """ $(TYPEDSIGNATURES) - Save `tf` as a h5 file to `filename` """ function save(filename::String, tf::TransferFunction) @@ -162,13 +182,11 @@ function save(filename::String, tf::TransferFunction) return nothing end -"""TODO: fix function and write docs""" -function load_tf_fromVNA(filename::String; - frequencyWeighting=false, - R = 50.0, #Ω - N = 10, #5# Turns - A = 7.4894*10.0^-4) #m^2 #1.3e-3^2*pi;) - +""" +$(TYPEDSIGNATURES) +Read the data recorded with the VNA from file `filename` into three Julia arrays `freq`, `ampdata`, `phasedata` +""" +function readVNAdata(filename::String) file = open(filename) lines = readlines(file) apdata = Float64[] @@ -231,11 +249,62 @@ function load_tf_fromVNA(filename::String; end end close(file) + return freq, apdata, aϕdata +end + +""" +$(TYPEDSIGNATURES) +Load data receive calibration from file recorded with the VNA. Data will be processed by `processRxTransferFunction`, see there for keyword arguments +""" +function load_tf_fromVNA(filename::String; kwargs...) + freq, ampdata, phasedata = readVNAdata(filename) + compdata = ampdata.*exp.(im.*phasedata) + if any([:R,:N,:A] .∈ [keys(kwargs)]) + @warn "In v0.16.1 and below load_tf_fromVNA mistakenly ignored the keyword parameters R, N and A. The current version includes these parametes if set, resulting in the correct scaling in magnetic moment domain." + end + return processRxTransferFunction(freq, compdata; kwargs...) +end + +""" +$(TYPEDSIGNATURES) +Process the data from a receive calibration measurement using a calibration coil into a TransferFunction + +Keyword parameters: +- `frequencyWeighting`: if true corrects for the frequency term in the TF, which results in a TF that does not integrate on application, but instead shows derivative of magnetic moment +- `R`: value of the resistance of the calibration coil in Ω +- `N`: number of turns of the calibration coil +- `A`, `r`, `d`: Area in m² or radius/diameter in m of the calibration coil, define only one +""" +function processRxTransferFunction(freq, compdata; frequencyWeighting::Bool=false, + R::Union{Real,Nothing} = nothing, # Ω + N::Union{Real,Nothing} = nothing, # Turns + A::Union{Real,Nothing} = nothing, # m² + r::Union{Real,Nothing} = nothing, # m + d::Union{Real,Nothing} = nothing) # m + + if sum(.!isnothing.([A,r,d])) > 1 + error("You can only define one of the keyword parameters A, r or d defining the geometry of the calibration coil") + elseif !isnothing(d) + A = pi * (d/2)^2 + elseif !isnothing(r) + A = pi * r^2 + end + + if all(isnothing.([R,N,A])) + unit = u"V/V" + @warn "No calibration coil parameters set in `processRxTransferFunction`. Make sure this is intended" + elseif any(isnothing.([R,N,A])) + error("To process the rxTransferFunction all three parameters describing the calibration coil (R, N and A) need to be defined.") + else + unit = u"V/A*m^2" + compdata .*= R/(N*A) + end + if frequencyWeighting - apdata ./= (freq.*2*pi) # As TF is defined as u_ADC = u_coil *TF the derivative from magnetic moment is applied as the division of the TF by w + compdata ./= (im*freq.*2*pi) # As TF is defined as u_ADC = u_coil *TF the derivative from magnetic moment is applied as the division of the TF by jw end - return TransferFunction(freq, apdata, aϕdata) + return TransferFunction(freq, compdata, units=[unit]) end @@ -290,9 +359,9 @@ function setTF(b::BrukerFile, filenameTF::AbstractString) return end +setTF(f::MDFFile, filenameTF::AbstractString) = setTF(f, TransferFunction(filenameTF), filenameTF) -function setTF(f::MDFFile, filenameTF::AbstractString) - tmf = TransferFunction(filenameTF) +function setTF(f::MDFFile, tmf::TransferFunction, filenameTF::Union{AbstractString,Nothing}=nothing) tf = sampleTF(tmf, f) # We need to close the HDF5 file handle before we can write to it @@ -302,7 +371,9 @@ function setTF(f::MDFFile, filenameTF::AbstractString) if haskey(file, "/acquisition/receiver/transferFunctionFileName") delete_object(file, "/acquisition/receiver/transferFunctionFileName") end - write(file, "/acquisition/receiver/transferFunctionFileName", filenameTF) + if !isnothing(filenameTF) + write(file, "/acquisition/receiver/transferFunctionFileName", filenameTF) + end if haskey(file, "/acquisition/receiver/transferFunction") delete_object(file, "/acquisition/receiver/transferFunction") end @@ -312,10 +383,13 @@ function setTF(f::MDFFile, filenameTF::AbstractString) end write(file, "/acquisition/receiver/inductionFactor", tmf.inductionFactor) end + # reopen filehandler so f stays usable + f.file = h5open(filepath(f), "r") return end function setTF(f::MDFv2InMemory, tf::TransferFunction) rxTransferFunction(f, sampleTF(tf, f)) rxInductionFactor(f, tf.inductionFactor) + return end \ No newline at end of file diff --git a/test/TransferFunction.jl b/test/TransferFunction.jl index 7d2b8a16..095c3c71 100644 --- a/test/TransferFunction.jl +++ b/test/TransferFunction.jl @@ -2,24 +2,28 @@ using MPIFiles using Test fnMeasBruker = joinpath(datadir,"BrukerStore","20150915_102110_Wuerfelphantom_1_1","18") +fnMeasMDFv1 = joinpath(datadir, "mdf", "measurement_V1.mdf") +fnTmpMDFv1 = joinpath(tmpdir, "mdf", "measurement_V1_tmp.mdf") +fnMeasMDFv2 = joinpath(tmpdir,"mdfim","measurement_V2.mdf") tfpath = joinpath(datadir,"transferFunction","example.s1p") -tfh5path = joinpath(tmpdir,"transferFunction","example.h5") +tfh5pathTmp = joinpath(tmpdir,"transferFunction","example.h5") +tfh5path = joinpath(datadir,"transferFunction","tf.h5") @testset "Testing TransferFunction submodule" begin a = TransferFunction(tfpath, frequencyWeighting=true) - if isfile(tfh5path) - rm(tfh5path) + if isfile(tfh5pathTmp) + rm(tfh5pathTmp) end - MPIFiles.save(tfh5path, a) - b = TransferFunction(tfh5path) + MPIFiles.save(tfh5pathTmp, a) + b = TransferFunction(tfh5pathTmp) @test a.freq == b.freq @test a.data == b.data @test a[1,1] == a(0.0,1) c = MPIFiles.combine(MPIFiles.combine(a,a),a) - rm(tfh5path) - MPIFiles.save(tfh5path, c) + rm(tfh5pathTmp) + MPIFiles.save(tfh5pathTmp, c) @test a[1,1] == c[1,2] @test a[1,1] == c[1,3] @@ -48,6 +52,19 @@ tfh5path = joinpath(tmpdir,"transferFunction","example.h5") @test tf(0) == 1u"V/V" @test tf(0,2) == 1u"A/V" + tf = TransferFunction(f, data) + tf_kHz = TransferFunction(f*u"kHz", data) + @test tf(100) ≈ tf_kHz(100e3) + + data_units = [1u"V/T" ./(1 .+im*f/1e4) 1u"V/(A*m^2)" ./(1 .+im*f/1e3)] + tf = TransferFunction(f, data_units) + @test_throws ErrorException TransferFunction([0,1], [1.0u"V/T" 1.0u"V/A"; 1.0u"V/T" 1.0u"m"]) + + tf_comp = TransferFunction([0,1,2], [1,im*2,-im*3]) + tf_deg = TransferFunction([0,1,2], [1,2,3], [0,90,270]u"°") + tf_rad = TransferFunction([0,1,2], [1,2,3], [0,π/2,-π/2]) + @test all(tf_comp([0,1,2]) .≈ tf_deg([0,1,2])) + @test all(tf_deg([0,1,2]) .≈ tf_rad([0,1,2])) f1 = collect(range(0,1e6,step=1e3)); data1 = 1 ./ (1 .+ im*f1/1e4 ); @@ -62,4 +79,46 @@ tfh5path = joinpath(tmpdir,"transferFunction","example.h5") tf_c = combine(tf1,tf2,interpolate=true) @test tf_c(101.5e3, :) ≈ [tf1(101.5e3) tf2(101.5e3)] atol=1e-12 + # prepare files + cp(fnMeasMDFv1, fnTmpMDFv1, force=true) + chmod(fnTmpMDFv1, 0o777) + f_mdfv1 = MPIFile(fnMeasMDFv1) + f_mdfv1_tmp = MPIFile(fnTmpMDFv1) + f_mdfv2 = MPIFile(fnMeasMDFv2) + + + # create TF from file + tf = TransferFunction(f_mdfv2) + @test sampleTF(tf, f_mdfv2) ≈ rxTransferFunction(f_mdfv2) + + # set TF to file that has no TF + @test_throws ErrorException getMeasurementsFD(f_mdfv1, tfCorrection=true) + newtf = TransferFunction([0,10e6],[1.0+0.0*im 1.0+0.0*im 1.0+0.0*im; 1.0+0.0*im 1.0+0.0*im 1.0+0.0*im]) + setTF(f_mdfv1_tmp, newtf) + @test rxHasTransferFunction(f_mdfv1_tmp) + @test getMeasurementsFD(f_mdfv1_tmp, tfCorrection=true) == getMeasurementsFD(f_mdfv1_tmp, tfCorrection=false) + + # read TF from h5 file + tf = TransferFunction(tfh5path) + # set TF using file path + setTF(f_mdfv1_tmp, tfh5path) + @test rxTransferFunction(f_mdfv1_tmp) == sampleTF(tf, f_mdfv1_tmp) + + setTF(MDFv2InMemory(f_mdfv2), newtf) + + # test printing + out = sprint() do io + show(io, MIME"text/plain"(),newtf) + end + @test out == "MPIFiles.TransferFunction: \n\t3 channel(s), units of [\"NoUnits\", \"NoUnits\", \"NoUnits\"]\n\t2 frequency samples from 0.0 Hz to 1.0e7 Hz" + + + # loading and processing tf + tf1 = MPIFiles.load_tf_fromVNA(tfpath, R=50, N=8, A=1e-3^2*pi) + tf2 = MPIFiles.load_tf_fromVNA(tfpath, R=50, N=8, d=2e-3) + tf3 = MPIFiles.load_tf_fromVNA(tfpath, R=50, N=8, r=1e-3) + @test tf1.data == tf2.data == tf3.data + @test_throws ErrorException MPIFiles.load_tf_fromVNA(tfpath, R=50, d=2e-3) + @test_throws ErrorException MPIFiles.load_tf_fromVNA(tfpath, R=50, N=8, d=2e-3, r=2e-3) + MPIFiles.load_tf_fromVNA(tfpath, R=50, N=8, d=2e-3, frequencyWeighting=true) end diff --git a/test/runtests.jl b/test/runtests.jl index 8263ae33..b07c75d1 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -33,10 +33,10 @@ include("MDFv1.jl") include("MultiMPIFile.jl") include("Reco.jl") include("IMT.jl") -include("TransferFunction.jl") -include("FrequencyFilter.jl") include("CustomSFMeas.jl") include("MDFInMemory.jl") +include("TransferFunction.jl") +include("FrequencyFilter.jl") include("MagneticFieldMeasurement.jl") @info "The unit tests are done!"