kramers kronig, fixed phase unwrapping

Kernel/Material.wl

@@ -59,20 +59,40 @@ MaterialParameters[n_Association]["Transmission"] := QuantityArray[Transpose[{n[
 
 MaterialParameters[n_Association]["Phase"] := QuantityArray[Transpose[{n["Frequencies"] // Normal, Normal[n["Phase"]]}], {1/"Centimeters", 1}]
 
+MaterialParameters[a_Association]["Phase Features"] := With[{
+  offset = offsetPhase[ a ],
+  raw = Normal[ a["Phase"] ]
+},
+  QuantityArray[Transpose[{Normal @ a["Frequencies"], raw - offset }], {1/"Centimeters", 1}]
+]
+
 
 MaterialParameters[n_Association]["Best Transmission"] := SelectBestRange[n] @ QuantityArray[Transpose[{n["Frequencies"] // Normal, Normal[n["T"]]^2}], {1/"Centimeters", 1}]
 
 MaterialParameters[n_Association]["Best Phase"] := SelectBestRange[n] @ QuantityArray[Transpose[{n["Frequencies"] // Normal, Normal[n["Phase"]]}], {1/"Centimeters", 1}]
 
 
+MaterialParameters[a_Association]["Best Phase Features"] := With[{
+  offset = offsetPhase[ a ],
+  raw = Normal[ a["Phase"] ]
+},
+  SelectBestRange[n] @ QuantityArray[Transpose[{Normal @ a["Frequencies"], raw - offset }], {1/"Centimeters", 1}]
+]
+
+
 MaterialParameters[n_Association]["n"] := QuantityArray[Transpose[{n["Frequencies"] // Normal, n["n"] // Normal}], {1/"Centimeters", 1}]
 MaterialParameters[n_Association]["k"] := QuantityArray[Transpose[{n["Frequencies"] // Normal, n["k"] // Normal}], {1/"Centimeters", 1}]
 
 MaterialParameters[n_Association]["Raw k"] := QuantityArray[Transpose[{n["Frequencies"] // Normal, n["k debug"] // Normal}], {1/"Centimeters", 1}]
 MaterialParameters[n_Association]["Best Raw k"] := SelectBestRange[n] @ QuantityArray[Transpose[{n["Frequencies"] // Normal, n["k debug"] // Normal}], {1/"Centimeters", 1}]
 MaterialParameters[n_Association]["Raw Best k"] := SelectBestRange[n] @ QuantityArray[Transpose[{n["Frequencies"] // Normal, n["k debug"] // Normal}], {1/"Centimeters", 1}]
 
-
+offsetPhase[t_] := With[{
+  offset = 2 Pi (1/33.356) QuantityMagnitude[t["\[Delta]t"], "Picoseconds"],
+  freqs = Normal[t["Frequencies"]]
+},
+  offset freqs 
+]
 
 
 SelectBestRange[n_][list_] := With[{ranges = findFDCIRanges[MaterialParameters[n] ]},
@@ -96,7 +116,7 @@ MaterialParameters[n_Association]["Domain"] := Quantity[#, 1/"Centimeters"] &/@
 
 
 MaterialParameters[a_Association]["FDCI2"] := 
-  With[{phase = MaterialParameters[a]["Phase"] // QuantityMagnitude}, 
+  With[{phase = QuantityMagnitude[MaterialParameters[a]["Phase"], {1/"Centimeters", 1}] }, 
     With[{r = Table[
         {phase[[q, 1]], 
           LinearModelFit[Take[phase, q], {1, x}, x]["RSquared"]}, 
@@ -395,10 +415,14 @@ materialParameters[list: List[__TransmissionObject], destCl_, srcCl_, metaCl_, {
 ] ]
 
 findFDCIRanges[mFp_MaterialParameters] := 
-  With[{m = mFp["Frequencies"] // QuantityMagnitude}, 
+  With[{
+    m = QuantityMagnitude[mFp["Frequencies"], 1/"Centimeters"] ,
+    left = QuantityMagnitude[mFp["FDCI"][[1]], 1/"Centimeters"],
+    right = QuantityMagnitude[mFp["FDCI"][[2]], 1/"Centimeters"]
+  }, 
     {
-      FirstPosition[m, x_ /; (x > QuantityMagnitude[mFp["FDCI"][[1]]])] // First,
-      Length[m] - FirstPosition[m // Reverse, x_ /; (x < QuantityMagnitude[mFp["FDCI"][[2]]])] // First
+      FirstPosition[m, x_ /; (x > left)] // First,
+      Length[m] - FirstPosition[m // Reverse, x_ /; (x < right)] // First
     }
   ];
 

Kernel/Trace.wl

@@ -58,11 +58,11 @@ TDTrace /: MakeBoxes[obj: TDTrace[n_NumericArray], StandardForm] := With[{
 
 TDTrace[n_NumericArray][s_String]    := (Message[TDTrace::unknownprop, s]; $Failed)
 TDTrace[n_NumericArray]["Trace"]    := QuantityArray[n // Normal, {"Picoseconds", 1}]
-TDTrace[n_NumericArray]["Spectrum"] := QuantityArray[fourier2d[TDTrace[n]["Trace"] // QuantityMagnitude], {1/"Centimeters", 1}]
+TDTrace[n_NumericArray]["Spectrum"] := QuantityArray[fourier2d[ QuantityMagnitude[TDTrace[n]["Trace"], {"Picoseconds", 1} ] ], {1/"Centimeters", 1}]
 
 
 TDTrace[n_NumericArray]["FrequencyDomainConfidenceInterval"] := Module[{ model, x0, \[Sigma], A},
-  With[{d = TDTrace[n]["PowerSpectrum"]//QuantityMagnitude}, 
+  With[{d = QuantityMagnitude[TDTrace[n]["PowerSpectrum"], {1/"Centimeters", 1}]}, 
     model = NonlinearModelFit[Drop[d,10], A Exp[-(*FB[*)(((*SpB[*)Power[(x0 - x)(*|*),(*|*)2](*]SpB*))(*,*)/(*,*)((*SpB[*)Power[\[Sigma](*|*),(*|*)2](*]SpB*)))(*]FB*)], {\[Sigma], A, x0}, x, ConfidenceLevel->0.5, Method->"NMinimize"];
     model = Association[model["BestFitParameters"]];
     Quantity[#, 1/"Centimeters"] &/@ {Clip[model[x0] - Abs[model[\[Sigma]]], {5.0, Infinity}], 1.2 Clip[model[x0] + Abs[model[\[Sigma]]], {30.0, Infinity}]}
@@ -71,7 +71,7 @@ TDTrace[n_NumericArray]["FrequencyDomainConfidenceInterval"] := Module[{ model,
 
 TDTrace[n_NumericArray]["FDCI"] := TDTrace[n]["FrequencyDomainConfidenceInterval"]
 
-TDTrace[n_NumericArray]["PowerSpectrum"] := QuantityArray[{#[[1]], Abs[#[[2]]]^2}&/@fourier2d[TDTrace[n]["Trace"] // QuantityMagnitude], {1/"Centimeters", 1}]
+TDTrace[n_NumericArray]["PowerSpectrum"] := QuantityArray[{#[[1]], Abs[#[[2]]]^2}&/@fourier2d[ QuantityMagnitude[TDTrace[n]["Trace"], {"Picoseconds", 1} ] ], {1/"Centimeters", 1}]
 
 TDTrace[n_NumericArray]["Properties"] := {"Properties", "Spectrum", "PowerSpectrum", "Trace", "FrequencyDomainConfidenceInterval", "FDCI"}
 

Kernel/Transmission.wl

@@ -43,20 +43,21 @@ TransmissionObject[sam_TDTrace, ref_TDTrace, opts: OptionsPattern[]] := Module[{
         freqs = fsam[[All,1]]
       },
         With[{
-          pDiff = 2 Pi (1/33.356) freqs (t0Sam - t0Ref)
+          pDiff = 2 Pi (1/33.356) freqs (t0Sam - t0Ref) 
         },
 
         TransmissionObject[<|
           "\[Delta]t"->Quantity[(t0Sam-t0Ref), "Picoseconds"], 
           "T"->NumericArray[t], 
           "Thickness"->Quantity[thickness, "Centimeters"],
+          "n0" -> 1.0 + (0.029979 (t0Sam-t0Ref) / thickness),
           "Gain"->gain,
           "FDCI"->sam["FDCI"],
           "Size" -> Length[freqs],
           "PhaseShift"->0,
           "Date" ->Now,
           "Tags" -> OptionValue["Tags"],
-          "Phase"->NumericArray[Arg[fsam[[All,2]]] - Arg[fref[[All,2]]] ], 
+          "Phase"->NumericArray[   Arg[Exp[I (Arg[fsam[[All,2]]] - Arg[fref[[All,2]]]) ] Exp[ -I pDiff ] ] + pDiff ], 
           "Frequencies"->NumericArray[freqs], (Sequence @@ Options[TransmissionObject]), 
           opts
         |>]
@@ -66,6 +67,21 @@ TransmissionObject[sam_TDTrace, ref_TDTrace, opts: OptionsPattern[]] := Module[{
   ]
 ]]
 
+(* by Thies Heidecke see https://mathematica.stackexchange.com/questions/341/implementing-discrete-and-continuous-hilbert-transforms *)
+HilbertSpectrum[0] := {};
+HilbertSpectrum[n_Integer?Positive] := With[{nhalf = Quotient[n - 1, 2]},
+  Join[{0}, ConstantArray[-I, nhalf], If[EvenQ[n], {0}, {}]
+          , ConstantArray[ I, nhalf]
+  ]
+];
+
+Hilbert[data_?VectorQ, padding_Integer?NonNegative] := Module[
+  {fp = FourierParameters -> {1, -1}, n = Length[data], m, paddeddata},
+  m = n + padding;
+  paddeddata = PadRight[data, m];
+  Re @ InverseFourier[ HilbertSpectrum[m] Fourier[paddeddata, fp], fp][[;;n]]
+]
+
 TransmissionObject /: MakeBoxes[obj: TransmissionObject[a_Association], StandardForm] := With[{
   preview = ListLinePlot[ArrayResample[Select[obj["Transmission"]//QuantityMagnitude//dropHalf, #[[2]]<1.0 &], 300], PlotStyle->ColorData[97][2], PlotRange->Full,Axes -> None, ImagePadding->None],
 
@@ -81,6 +97,7 @@ TransmissionObject /: MakeBoxes[obj: TransmissionObject[a_Association], Standard
           {BoxForm`SummaryItem[{"Gain ", obj["Gain"]}]},
           {BoxForm`SummaryItem[{"PhaseShift ", 2Pi obj["PhaseShift"]}]},
           {BoxForm`SummaryItem[{"Phase ", state}]},
+          {BoxForm`SummaryItem[{"n0 ", obj["n0"]}]},
           If[Length[obj["Tags"]//Keys] > 0, {BoxForm`SummaryItem[{"Tags ", Style[#, Background->LightBlue]&/@obj["Tags"]//Keys}]}, Nothing]
         };
 
@@ -103,31 +120,98 @@ TransmissionObject[a_Association][prop_String] := If[!KeyExistsQ[a, prop],
   a[prop]
 ]
 
+TransmissionObject[a_Association]["Kramers-Kronig n"] := With[{
+  n0 = a["n0"],
+  thickness = QuantityMagnitude[a["Thickness"], "Centimeters"],
+  freqs = Normal @ a["Frequencies"]
+},
+  With[{
+    im\[Chi] =  (- (*FB[*)((1)(*,*)/(*,*)(#[[1]] thickness))(*]FB*) Log[#[[2]]]) &/@ Drop[QuantityMagnitude[TransmissionObject[a]["Transmission"], {1/"Centimeters", 1} ], 1]
+  },
+    QuantityArray[{freqs, Join[{n0}, (*SqB[*)Sqrt[n0^2 - Hilbert[im\[Chi], 0] ](*]SqB*)]} // Transpose, {1/"Centimeters", 1}]
+  ]
+]
+
 TransmissionObject[a_Association]["Transmission"] := QuantityArray[Transpose[{Normal @ a["Frequencies"], (*SpB[*)Power[(a["Gain"] Normal @ a["T"])(*|*),(*|*)2](*]SpB*)}], {1/"Centimeters", 1}]
 
-TransmissionObject /: Append[TransmissionObject[a_Association], props_Association] := TransmissionObject[Join[a, props] ]
-TransmissionObject /: Append[TransmissionObject[a_Association], prop_Rule] := TransmissionObject[Append[a, prop] ]
-TransmissionObject /: Append[TransmissionObject[a_Association], props_List] := TransmissionObject[Append[a, props] ]
+updateThicknessDependent[a_Association ] := With[{
+
+},
+  Join[a, <|
+    "n0" -> 1.0 + (0.029979 QuantityMagnitude[a["\[Delta]t"], "Picoseconds" ] / QuantityMagnitude[a["Thickness"], "Centimeters"]),
+    "Date" ->Now
+  |>] 
+]
+
+TransmissionObject /: Append[TransmissionObject[a_Association], props_Association] := TransmissionObject[Join[a, props] // updateThicknessDependent ]
+TransmissionObject /: Append[TransmissionObject[a_Association], prop_Rule] := TransmissionObject[Append[a, prop] // updateThicknessDependent ]
+TransmissionObject /: Append[TransmissionObject[a_Association], props_List] := TransmissionObject[Append[a, props] // updateThicknessDependent ]
 
 TransmissionObject[a_Association]["Frequencies"] := QuantityArray[Normal @ a["Frequencies"], 1/"Centimeters"]
 
 TransmissionObject[a_Association]["Domain"] := Quantity[#, 1/"Centimeters"] &/@ MinMax[Normal @ a["Frequencies"]]
 
 
+offsetPhase[t_] := With[{
+  offset = 2 Pi (1/33.356) QuantityMagnitude[t["\[Delta]t"], "Picoseconds"],
+  freqs = Normal[t["Frequencies"]]
+},
+  offset freqs 
+]
+
 TransmissionObject[a_Association]["Phase"] := With[{
   shift = 2 Pi (1/33.356) QuantityMagnitude[a["\[Delta]t"], "Picoseconds"] Normal[a["Frequencies"]],
   off = a["PhaseShift"]
 },
   QuantityArray[Transpose[{Normal @ a["Frequencies"], (Normal @ a["Phase"])}], {1/"Centimeters", 1}]
 ]
 
+TransmissionObject[a_Association]["Phase Features"] := With[{
+  offset = offsetPhase[ a ],
+  raw = Normal[ a["Phase"] ]
+},
+  QuantityArray[Transpose[{Normal @ a["Frequencies"], raw - offset }], {1/"Centimeters", 1}]
+]
+
+
+TransmissionObject[a_Association]["Approximated k"] := With[{
+  thickness = QuantityMagnitude[a["Thickness"], "Centimeters"],
+  gain = a["Gain"]
+},
+  QuantityArray[
+    Join[{0, 0}, {#[[1]], - 0.159152 Log[#[[2]] gain ] / (#[[1]] thickness) } &/@ Drop[Transpose[Normal /@ {a["Frequencies"], a["T"]}], 1] ]
+  , {1/"Centimeters", 1}]
+]
+
+TransmissionObject[a_Association]["Approximated \[Alpha]"] := With[{
+  thickness = QuantityMagnitude[a["Thickness"], "Centimeters"],
+  gain = a["Gain"]
+},
+  QuantityArray[
+    Join[{0, 0}, {#[[1]], ((- 0.159152 Log[#[[2]] gain ] / (#[[1]] thickness)) 4 \[Pi]  10^12 #[[1]])/(33.356 2.9979 10^10)} &/@ Drop[Transpose[Normal /@ {a["Frequencies"], a["T"]}], 1] ]
+  , {1/"Centimeters", 1/"Centimeters"}]
+]
+
+
+TransmissionObject[a_Association]["Approximated n"] := With[{
+  shift = a["PhaseShift"],
+  thickness = QuantityMagnitude[a["Thickness"], "Centimeters"],
+  n0 = a["n0"]
+},
+  QuantityArray[
+    Join[{0, n0}, {#[[1]], 1.0 + 0.159152 (#[[2]] + shift) / (#[[1]] thickness) } &/@ Drop[Transpose[Normal /@ {a["Frequencies"], a["Phase"]}],1] ]
+  , {1/"Centimeters", 1}]
+]
+
+
+
 TransmissionObject[a_Association]["FrequencyDomainConfidenceInterval"] := TransmissionObject[a]["FDCI"]
 
 TransmissionObject[a_Association]["FrequencyDomainConfidenceInterval2"] := TransmissionObject[a]["FDCI2"]
 
 
 TransmissionObject[a_Association]["FDCI2"] := 
-With[{phase = TransmissionObject[a]["Phase"] // QuantityMagnitude},
+With[{phase = QuantityMagnitude[TransmissionObject[a]["Phase"], {1/"Centimeters", 1}]},
   With[{r = Table[{phase[[q, 1]], 
       LinearModelFit[Take[phase, q], {1, x}, x]["RSquared"]}, 
      {q, Round[0.2 (phase // Length)], phase // Length, 25}]},
@@ -150,7 +234,7 @@ phaseState[phase_List] := With[{},
 
 TransmissionObject /: Keys[t_TransmissionObject] :=  t["Properties"]
 
-TransmissionObject[a_Association]["Properties"] := Join[Options[TransmissionObject][[All,1]], {"Properties", "Frequencies", "Transmission", "Phase", "\[Delta]t", "Gain", "PhaseShift", "Thickness", "Domain", "FrequencyDomainConfidenceInterval", "FDCI", "FDCI2", "FrequencyDomainConfidenceInterval2"}]
+TransmissionObject[a_Association]["Properties"] := Join[Options[TransmissionObject][[All,1]], {"Properties", "n0", "Approximated n", "Approximated k", "Approximated \[Alpha]", "Kramers-Kronig n", "Frequencies", "Transmission", "Phase", "Phase Features", "\[Delta]t", "Gain", "PhaseShift", "Thickness", "Domain", "FrequencyDomainConfidenceInterval", "FDCI", "FDCI2", "FrequencyDomainConfidenceInterval2"}]
 
 Options[TransmissionObject] = {"Thickness"->Null, "Tags"-><||>, "Gain"->1.0, "PhaseShift"->0};
 
@@ -171,9 +255,8 @@ validateOptions[OptionsPattern[] ] := With[{},
 ]
 
 TransmissionUnwrap[t: TransmissionObject[a_], "Basic", OptionsPattern[]] := With[{
-  offset = 2 Pi (1/33.356) QuantityMagnitude[t["\[Delta]t"], "Picoseconds"] ,
+  offset = offsetPhase[a],
   th = OptionValue["PhaseThreshold"]//N,
-  freqs = Normal[a["Frequencies"]],
   phaseShift = OptionValue["PhaseShift"]
 },
   If[!NumericQ[OptionValue["PhaseThreshold"] ] || !NumericQ[OptionValue["PhaseShift"] ],
@@ -182,9 +265,9 @@ TransmissionUnwrap[t: TransmissionObject[a_], "Basic", OptionsPattern[]] := With
   ];
 
   With[{unwrapped = With[{
-    phase = Normal[a["Phase"]]
+    phase = Normal[a["Phase"] ]
   },
-    clusterPhase[phase - offset freqs, 1, Length[phase]-1, th] + offset freqs
+    clusterPhase[phase - offset, 1, Length[phase]-1, th] + offset
   ]},
     Append[t, {"Phase"->NumericArray[unwrapped], "PhaseShift"->phaseShift}]
   ]

PacletInfo.wl

@@ -7,7 +7,7 @@ PacletObject[
     "Creator" -> "Kirill Vasin",
     "License" -> "MIT",
     "PublisherID" -> "JerryI",
-    "Version" -> "0.0.3",
+    "Version" -> "0.0.4",
     "WolframVersion" -> "13+",
     "PrimaryContext" -> "JerryI`TDSTools`",
     "Extensions" -> {

README.md

@@ -12,9 +12,10 @@ A small library for high-precision material parameter extraction from time-domai
 ## Features
 - Fabry-Pérot deconvolution ⭐️
 - Informed phase unwrapping
-- High precision for $n$, $\kappa$, and $\alpha$ solving
+- High precision / various approximation methods for $n$, $\kappa$, and $\alpha$ solving
 - Works for both thin and thick samples
 - **GPU Acceleration** (OpenCL) ⭐️
+- Kramers-Kronig approximation of $n$ feature (if needed)
 - Functional approach, no hidden state
 - Syntax sugar with data preview ⭐️
 - Jitter-proof
@@ -156,10 +157,18 @@ One can provide the following `opts` to the constructor
 - `"\[Delta]t"` : returns estimated time-delay between the sample and the reference signals. *Read-only*
 - `"Gain"` : returns a multiplier for `sam` signal
 - `"PhaseShift"` : returns a constanst offset of the phase in multples of `2Pi`
+- `"n0"` : returns DC refractive index, it will be updated if one change the thickness of the material using append *Read-only*
+- `"Date"` : returns date, when an object was modified last time
 ---
 - `"Frequencies"` : returns `QuantityArray` of the whole range of frequencies. *Read-only*
 - `"Transmission"` : returns the power transmission `QuantityArray`, i.e. array of $|\hat{t}(\omega)|^2$. *Read-only*
 - `"Phase"` : returns the argument of $\hat{t}$ sampled as a function of frequency in a form of `QuantityArray`. After the creation is wrapped and cannot be used before unwrapping process (see later). *Read-only*
+- `"Phase Features"` : returns quantity array of `"Phase"` with subtracted DC contribution of the refractive index
+- `"Kramers-Kronig n"` : (approximatiton!) returns quantity array of refractive index derived from the transmission using Kramers-Kronig relation (discrete)
+---
+- `"Approximated n"` : (approximatiton!) returns quantity array of refractive index derived from the phase directly
+- `"Approximated k"` : (approximatiton!) returns quantity array of imaginary refractive index derived from the transmission directly
+- `"Approximated \[Alpha]"` : (approximatiton!) returns quantity array of absorption coefficient derived from the transmission directly
 ---
 - `"Domain"` : returns the ranges of frequencies. *Read-only*
 - `"FDCI"` : a range of frequency-domain confidence interval. It is generated according to the powerspectrum and provides a coarse estimation of the working range in wavenumbers. This interval will be used later for optimizing material parameters. *Read-only*
@@ -231,8 +240,10 @@ All properties are *read-only*.
 - `"Frequencies"` : returns `QuantityArray` of frequencies
 - `"Transmission"` : returns `QuantityArray` of a power transmission. If `FabryPerotCancellation` is applied, it will return deconvoluted spectrum. 
 - `"Best Transmission"` : returns `"Transmission"` selected in the region of `"FDCI"`
-- `"Phase"` : returns `QuantityArray` of the phase of the transmissino function. 
+- `"Phase"` : returns `QuantityArray` of the phase of the transmissino function.
+- `"Phase Features"` : returns quantity array of `"Phase"` with subtracted DC contribution of the refractive index 
 - `"Best Phase"` : returns `"Phase"` selected in the region of `"FDCI"`
+- `"Best Phase Features"` : returns quantity array of `"Phase"` with subtracted DC contribution of the refractive index  in the region of `"FDCI"`
 ---
 - `"\[Alpha]"` : returns `QuantityArray` of extracted absorption coefficient.
 - `"Raw \[Alpha]"` : returns `"\[Alpha]"` without FP cancellation applied.