correct voigt amplitude

deconvolution.py

@@ -1,13 +1,12 @@
 from dataclasses import dataclass, field
-from pathlib import Path
 from typing import Callable
 
 import click
 import labmodule as lab
 import numpy as np
 from customfunc.app import PlotApp
-from numpy.fft import fft, fftshift, ifft
 from scipy.optimize import curve_fit
+from scipy.interpolate import interp1d
 
 LN2 = np.log(2)
 
@@ -32,6 +31,7 @@ class VoigtFitResult:
     name: str
     curve: np.ndarray = field(repr=False)
     covariance: np.ndarray
+    amp: float
     gaussian_width: float
     lorentzian_width: float
 
@@ -114,21 +114,23 @@ def fft_gaussian(f: np.ndarray, dx: float, amp: float, width: float) -> np.ndarr
     return gaussian(f, freq_amp, freq_width, 0)
 
 
-def voigt(x: np.ndarray, gaussian_width: float, lorentzian_width: float) -> np.ndarray:
+def voigt(x: np.ndarray, amp: float, gaussian_width: float, lorentzian_width: float) -> np.ndarray:
     """Voigt profile of max 1 with specified Gaussian and Lorentzian FWHM"""
     conv = convolve(unit_gaussian(x, gaussian_width, 0), unit_lorentzian(x, lorentzian_width, 0))
     if conv.max() > 0:
-        return conv / conv.max()
+        return conv * amp / conv.max()
     return conv
 
 
-def log_voigt(x: np.ndarray, gaussian_width: float, lorentzian_width: float) -> np.ndarray:
+def log_voigt(
+    x: np.ndarray, amp: float, gaussian_width: float, lorentzian_width: float
+) -> np.ndarray:
     """log of `voigt`"""
-    return np.log(voigt(x, gaussian_width, lorentzian_width))
+    return np.log(voigt(x, amp, gaussian_width, lorentzian_width))
 
 
 def fft_voigt(
-    f: np.ndarray, dx: float, gaussian_width: float, lorentzian_width: float
+    f: np.ndarray, dx: float, amp: float, gaussian_width: float, lorentzian_width: float
 ) -> np.ndarray:
     """
     returns the product of a `peak_func` and a `gaussian`
@@ -136,7 +138,8 @@ def fft_voigt(
     given by `gaussian_width` and `lorentzian_width`
     """
     profile = fft_gaussian(f, dx, 1, gaussian_width) * fft_lorentzian(f, dx, 1, lorentzian_width)
-    return profile / profile.max()
+    fac = amp / np.abs(np.fft.ifft(np.fft.fftshift(profile))).max()
+    return profile * fac
 
 
 def single_func_fit(
@@ -192,13 +195,13 @@ def log_direct_voigt_fit(x: np.ndarray, y: np.ndarray) -> VoigtFitResult:
     y_fit = np.log(y[ind])
     x_fit = x[ind]
     params, cov = curve_fit(log_voigt, x_fit, y_fit, bounds=[1e-12, np.inf])
-    curve = np.exp(log_voigt(x, *params))
+    curve = voigt(x, *params)
     return VoigtFitResult("Log direct fit", curve, cov, *params)
 
 
 def linear_fft_voigt_fit(f: np.ndarray, A: np.ndarray, x: np.ndarray, dx: float) -> VoigtFitResult:
     params, cov = curve_fit(
-        lambda _f, _g, _l: fft_voigt(_f, dx, _g, _l), f, A, bounds=[1e-12, np.inf]
+        lambda _f, _a, _g, _l: fft_voigt(_f, dx, _a, _g, _l), f, A, bounds=[1e-12, np.inf]
     )
     curve = voigt(x, *params)
     return VoigtFitResult("Linear FFT fit", curve, cov, *params)
@@ -335,25 +338,33 @@ def demo3():
 
 
 @click.argument("file", type=click.Path(True, dir_okay=False))
+@click.option("--all", "show_all", default=False, is_flag=True)
 @main.command(name="fit")
-def fit_measurement(file):
+def fit_measurement(file, show_all):
-    wl, intens = lab.osa.load_osa_spectrum(file)
+    raw_wl, intens = lab.osa.load_osa_spectrum(file)
-    wl *= 1e9
+    raw_wl *= 1e9
-    gauss_fit = single_func_fit(wl, intens, elevated_gaussian, extra_params=1)
+    gauss_fit = single_func_fit(raw_wl, intens, elevated_gaussian, extra_params=1)
     wl0 = gauss_fit.x0
     baseline = gauss_fit.params[-1]
     intens -= baseline
     gauss_fit.curve -= baseline
     intens[intens < 0] = 0
-    wl -= wl0
+
-    with PlotApp(n=np.arange(len(wl) // 2)[1:]) as app:
+    lim = np.max(np.abs(raw_wl - wl0))
+    new_wl = np.linspace(-lim, lim, 2048)
+    intens = interp1d(raw_wl - wl0, intens, bounds_error=False, fill_value=0)(new_wl)
+    wl = new_wl
+
+    with PlotApp(n=np.arange(len(wl) // 2)) as app:
         app.set_antialiasing(True)
 
         ax1 = app["Measurement"]
         ax2 = app["Fourier domain"]
         ax1.set_line_data("Measured spectrum", wl, intens, width=2, color="red")
-        ax1.set_line_data("Initial Guassian fit", wl, gauss_fit.curve, width=2, color="blue")
+        ax1.set_line_data(
-        ax1.plot_widget.plotItem.setLogMode(y=True)
+            "Initial Guassian fit", raw_wl - wl0, gauss_fit.curve, width=2, color="blue"
+        )
+        # ax1.plot_widget.plotItem.setLogMode(y=True)
 
         @app.update
         def draw(n):
@@ -366,35 +377,63 @@ def fit_measurement(file):
             ax1.set_line_data("Fitted spectrum", wl_fit, intens_fit)
 
             trans = np.abs(np.fft.fft(intens_fit))[ind]
-            m = trans.max()
             ax2.set_line_data("Transformed", f, trans)
+            ax2.lines["Transformed"].setZValue(10)
 
-            fft_fit = linear_fft_voigt_fit(f, trans / m, wl_fit, dwl)
+            fft_fit = linear_fft_voigt_fit(f, trans, wl_fit, dwl)
             lin_fit = linear_direct_voigt_fit(wl_fit, intens_fit)
             log_fit = log_direct_voigt_fit(wl_fit, intens_fit)
             for fit in fft_fit, lin_fit, log_fit:
-                display_fit(wl_fit, dwl, f, fit, m)
+                display_fit(wl_fit, dwl, f, fit)
 
-        def display_fit(x: np.ndarray, dx: float, f: np.ndarray, fit: VoigtFitResult, amp: float):
+        def display_fit(x: np.ndarray, dx: float, f: np.ndarray, fit: VoigtFitResult):
             ax2.set_line_data(
-                fit.name, f, amp * fft_voigt(f, dx, fit.gaussian_width, fit.lorentzian_width)
+                fit.name,
+                f,
+                fft_voigt(f, dx, fit.amp, fit.gaussian_width, fit.lorentzian_width),
             )
             ax1.set_line_data(fit.name, x, fit.curve)
+            if show_all:
                 ax1.set_line_data(
                     fit.name + "gauss",
                     x,
                     gaussian(x, 1, fit.gaussian_width, 0),
-                label=f"{fit.name} - Gaussian {fit.gaussian_width:.2f}",
+                    label=f"{fit.name} - Gaussian {fit.gaussian_width:.3f}",
                     width=1.5,
                 )
                 ax1.set_line_data(
                     fit.name + "lorentz",
                     x,
                     lorentzian(x, 1, fit.lorentzian_width, 0),
-                label=f"{fit.name} - Lorentzian {fit.lorentzian_width:.2f}",
+                    label=f"{fit.name} - Lorentzian {fit.lorentzian_width:.3f}",
                     width=1.5,
                 )
 
 
+@main.command(name="test", help="test scaling of voigt profile")
+def test_amp():
+    x, dx, f, df = setup_axes(1001, (-20, 20))
+    with PlotApp(
+        lorentz_fwhm=np.geomspace(0.1, 1), gauss_fwhm=np.linspace(1, 10), amp=np.linspace(1, 7)
+    ) as app:
+        app.set_antialiasing(True)
+
+        ax1 = app["Main"]
+        ax2 = app["Fourrier"]
+        ind = np.argsort(f)
+
+        @app.update
+        def draw(lorentz_fwhm, gauss_fwhm, amp):
+            vo = voigt(x, amp, gauss_fwhm, lorentz_fwhm)
+            trans = np.fft.fft(np.fft.fftshift(vo))
+            expected = fft_voigt(f, dx, amp, gauss_fwhm, lorentz_fwhm)
+            back_transformed = np.fft.fftshift(np.fft.ifft(expected).real)
+
+            ax1.set_line_data("original", x, vo)
+            ax1.set_line_data("back tranformed", x, back_transformed)
+            ax2.set_line_data("expected", f[ind], expected[ind])
+            ax2.set_line_data("original transformed", f[ind], trans[ind].real)
+
+
 if __name__ == "__main__":
     main()