diff --git a/deconvolution.py b/deconvolution.py index 72f87c1..def8257 100644 --- a/deconvolution.py +++ b/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): - wl, intens = lab.osa.load_osa_spectrum(file) - wl *= 1e9 - gauss_fit = single_func_fit(wl, intens, elevated_gaussian, extra_params=1) +def fit_measurement(file, show_all): + raw_wl, intens = lab.osa.load_osa_spectrum(file) + raw_wl *= 1e9 + 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.plot_widget.plotItem.setLogMode(y=True) + ax1.set_line_data( + "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,34 +377,62 @@ 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) - ax1.set_line_data( - fit.name + "gauss", - x, - gaussian(x, 1, fit.gaussian_width, 0), - label=f"{fit.name} - Gaussian {fit.gaussian_width:.2f}", - 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}", - width=1.5, - ) + 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:.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:.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__":