Coverage for larch/xafs/deconvolve.py: 16%
57 statements
« prev ^ index » next coverage.py v7.6.0, created at 2024-10-16 21:04 +0000
« prev ^ index » next coverage.py v7.6.0, created at 2024-10-16 21:04 +0000
1#!/usr/bin/env python
2# XAS spectral decovolution
3#
5import numpy as np
6from scipy.signal import deconvolve
7from larch import parse_group_args, Make_CallArgs
9from larch.math import (gaussian, lorentzian, interp,
10 index_of, index_nearest, remove_dups,
11 savitzky_golay)
13from .xafsutils import set_xafsGroup, TINY_ENERGY
15@Make_CallArgs(["energy","norm"])
16def xas_deconvolve(energy, norm=None, group=None, form='lorentzian',
17 esigma=1.0, eshift=0.0, smooth=True,
18 sgwindow=None, sgorder=3, _larch=None):
19 """XAS spectral deconvolution
21 de-convolve a normalized mu(E) spectra with a peak shape, enhancing the
22 intensity and separation of peaks of a XANES spectrum.
24 The results can be unstable, and noisy, and should be used
25 with caution!
27 Arguments
28 ----------
29 energy: array of x-ray energies (in eV) or XAFS data group
30 norm: array of normalized mu(E)
31 group: output group
32 form: functional form of deconvolution function. One of
33 'gaussian' or 'lorentzian' [default]
34 esigma energy sigma to pass to gaussian() or lorentzian()
35 [in eV, default=1.0]
36 eshift energy shift to apply to result. [in eV, default=0]
37 smooth whether to smooth result with savitzky_golay method [True]
38 sgwindow window size for savitzky_golay [found from data step and esigma]
39 sgorder order for savitzky_golay [3]
41 Returns
42 -------
43 None
44 The array 'deconv' will be written to the output group.
46 Notes
47 -----
48 Support See First Argument Group convention, requiring group
49 members 'energy' and 'norm'
51 Smoothing with savitzky_golay() requires a window and order. By
52 default, window = int(esigma / estep) where estep is step size for
53 the gridded data, approximately the finest energy step in the data.
54 """
55 energy, mu, group = parse_group_args(energy, members=('energy', 'norm'),
56 defaults=(norm,), group=group,
57 fcn_name='xas_deconvolve')
58 eshift = eshift + 0.5 * esigma
60 en = remove_dups(energy, tiny=TINY_ENERGY)
61 estep1 = int(0.1*en[0]) * 2.e-5
62 en = en - en[0]
63 estep = max(estep1, 0.01*int(min(en[1:]-en[:-1])*100.0))
64 npts = 1 + int(max(en) / estep)
65 if npts > 25000:
66 npts = 25001
67 estep = max(en)/25000.0
69 x = np.arange(npts)*estep
70 y = interp(en, mu, x, kind='cubic')
72 kernel = lorentzian
73 if form.lower().startswith('g'):
74 kernel = gaussian
76 yext = np.concatenate((y, np.arange(len(y))*y[-1]))
78 ret, err = deconvolve(yext, kernel(x, center=0, sigma=esigma))
79 nret = min(len(x), len(ret))
81 ret = ret[:nret]*yext[nret-1]/ret[nret-1]
82 if smooth:
83 if sgwindow is None:
84 sgwindow = int(1.0*esigma/estep)
86 sqwindow = int(sgwindow)
87 if sgwindow < (sgorder+1):
88 sgwindow = sgorder + 2
89 if sgwindow % 2 == 0:
90 sgwindow += 1
91 ret = savitzky_golay(ret, sgwindow, sgorder)
93 out = interp(x+eshift, ret, en, kind='cubic')
94 group = set_xafsGroup(group, _larch=_larch)
95 group.deconv = out
98@Make_CallArgs(["energy","norm"])
99def xas_convolve(energy, norm=None, group=None, form='lorentzian',
100 esigma=1.0, eshift=0.0, _larch=None):
101 """
102 convolve a normalized mu(E) spectra with a Lorentzian or Gaussian peak
103 shape, degrading separation of XANES features.
105 This is provided as a complement to xas_deconvolve, and to deliberately
106 broaden spectra to compare with spectra measured at lower resolution.
108 Arguments
109 ----------
110 energy: array of x-ray energies (in eV) or XAFS data group
111 norm: array of normalized mu(E)
112 group: output group
113 form: form of deconvolution function. One of
114 'lorentzian' or 'gaussian' ['lorentzian']
115 esigma energy sigma (in eV) to pass to gaussian() or lorentzian() [1.0]
116 eshift energy shift (in eV) to apply to result [0]
118 Returns
119 -------
120 None
121 The array 'conv' will be written to the output group.
123 Notes
124 -----
125 Follows the First Argument Group convention, using group members named
126 'energy' and 'norm'
127 """
129 energy, mu, group = parse_group_args(energy, members=('energy', 'norm'),
130 defaults=(norm,), group=group,
131 fcn_name='xas_convolve')
132 eshift = eshift + 0.5 * esigma
134 en = remove_dups(energy, tiny=TINY_ENERGY)
135 en = en - en[0]
136 estep = max(0.001, 0.001*int(min(en[1:]-en[:-1])*1000.0))
138 npad = 1 + int(max(estep*2.01, 50*esigma)/estep)
140 npts = npad + int(max(en) / estep)
142 x = np.arange(npts)*estep
143 y = interp(en, mu, x, kind='cubic')
145 kernel = lorentzian
146 if form.lower().startswith('g'):
147 kernel = gaussian
149 k = kernel(x, center=0, sigma=esigma)
150 ret = np.convolve(y, k, mode='full')
152 out = interp(x-eshift, ret[:len(x)], en, kind='cubic')
154 group = set_xafsGroup(group, _larch=_larch)
155 group.conv = out / k.sum()