Coverage for larch/xafs/estimate_noise.py: 25%
28 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"""
3 Estimate Noise in an EXAFS spectrum
4"""
5from numpy import pi, sqrt, where
7from larch import parse_group_args, Group, isgroup, Make_CallArgs
9from larch.math import index_of, realimag
10from .xafsutils import set_xafsGroup
11from .xafsft import xftf, xftr
13@Make_CallArgs(["k","chi"])
14def estimate_noise(k, chi=None, group=None, rmin=15.0, rmax=30.0,
15 kweight=1, kmin=0, kmax=20, dk=4, dk2=None, kstep=0.05,
16 kwindow='kaiser', nfft=2048, _larch=None, **kws):
17 """
18 estimate noise levels in EXAFS spectrum and estimate highest k
19 where data is above the noise level
20 Parameters:
21 -----------
22 k: 1-d array of photo-electron wavenumber in Ang^-1 (or group)
23 chi: 1-d array of chi
24 group: output Group [see Note below]
25 rmin: minimum R value for high-R region of chi(R)
26 rmax: maximum R value for high-R region of chi(R)
27 kweight: exponent for weighting spectra by k**kweight [1]
28 kmin: starting k for FT Window [0]
29 kmax: ending k for FT Window [20]
30 dk: tapering parameter for FT Window [4]
31 dk2: second tapering parameter for FT Window [None]
32 kstep: value to use for delta_k ( Ang^-1) [0.05]
33 window: name of window type ['kaiser']
34 nfft: value to use for N_fft [2048].
36 Returns:
37 ---------
38 None -- outputs are written to supplied group. Values (scalars) written
39 to output group:
40 epsilon_k estimated noise in chi(k)
41 epsilon_r estimated noise in chi(R)
42 kmax_suggest highest estimated k value where |chi(k)| > epsilon_k
44 Notes:
45 -------
47 1. This method uses the high-R portion of chi(R) as a measure of the noise
48 level in the chi(R) data and uses Parseval's theorem to convert this noise
49 level to that in chi(k). This method implicitly assumes that there is no
50 signal in the high-R portion of the spectrum, and that the noise in the
51 spectrum s "white" (independent of R) . Each of these assumptions can be
52 questioned.
53 2. The estimate for 'kmax_suggest' has a tendency to be fair but pessimistic
54 in how far out the chi(k) data goes before being dominated by noise.
55 3. Follows the 'First Argument Group' convention, so that you can either
56 specifiy all of (an array for 'k', an array for 'chi', option output Group)
57 OR pass a group with 'k' and 'chi' as the first argument
58 """
59 k, chi, group = parse_group_args(k, members=('k', 'chi'),
60 defaults=(chi,), group=group,
61 fcn_name='esitmate_noise')
64 tmpgroup = Group()
65 rmax_out = min(10*pi, rmax+2)
67 xftf(k, chi, kmin=kmin, kmax=kmax, rmax_out=rmax_out,
68 kweight=kweight, dk=dk, dk2=dk2, kwindow=kwindow,
69 nfft=nfft, kstep=kstep, group=tmpgroup)
71 chir = tmpgroup.chir
72 rstep = tmpgroup.r[1] - tmpgroup.r[0]
74 irmin = int(0.01 + rmin/rstep)
75 irmax = min(nfft/2, int(1.01 + rmax/rstep))
76 highr = realimag(chir[irmin:irmax])
78 # get average of window function value, scale eps_r scale by this
79 # this is imperfect, but improves the result.
80 kwin_ave = tmpgroup.kwin.sum()*kstep/(kmax-kmin)
81 eps_r = sqrt((highr*highr).sum() / len(highr)) / kwin_ave
83 # use Parseval's theorem to convert epsilon_r to epsilon_k,
84 # compensating for kweight
85 w = 2 * kweight + 1
86 scale = sqrt((2*pi*w)/(kstep*(kmax**w - kmin**w)))
87 eps_k = scale*eps_r
89 # do reverse FT to get chiq array
90 xftr(tmpgroup.r, tmpgroup.chir, group=tmpgroup, rmin=0.5, rmax=9.5,
91 dr=1.0, window='parzen', nfft=nfft, kstep=kstep)
93 # sets kmax_suggest to the largest k value for which
94 # | chi(q) / k**kweight| > epsilon_k
95 iq0 = index_of(tmpgroup.q, (kmax+kmin)/2.0)
96 tst = tmpgroup.chiq_mag[iq0:] / ( tmpgroup.q[iq0:])**kweight
97 kmax_suggest = tmpgroup.q[iq0 + where(tst < eps_k)[0][0]]
99 # restore original _sys.xafsGroup, set output variables
101 group.epsilon_k = eps_k
102 group.epsilon_r = eps_r
103 group.kmax_suggest = kmax_suggest