1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
|
"""Baseline estimation algorithms."""
import numpy as np
import scipy.linalg as LA
import math
def baseline(y, deg=None, max_it=None, tol=None):
"""
Computes the baseline of a given data.
Iteratively performs a polynomial fitting in the data to detect its
baseline. At every iteration, the fitting weights on the regions with
peaks are reduced to identify the baseline only.
Parameters
----------
y : ndarray
Data to detect the baseline.
deg : int (default: 3)
Degree of the polynomial that will estimate the data baseline. A low
degree may fail to detect all the baseline present, while a high
degree may make the data too oscillatory, especially at the edges.
max_it : int (default: 100)
Maximum number of iterations to perform.
tol : float (default: 1e-3)
Tolerance to use when comparing the difference between the current
fit coefficients and the ones from the last iteration. The iteration
procedure will stop when the difference between them is lower than
*tol*.
Returns
-------
ndarray
Array with the baseline amplitude for every original point in *y*
"""
# for not repeating ourselves in `envelope`
if deg is None: deg = 3
if max_it is None: max_it = 100
if tol is None: tol = 1e-3
order = deg + 1
coeffs = np.ones(order)
# try to avoid numerical issues
cond = math.pow(abs(y).max(), 1. / order)
x = np.linspace(0., cond, y.size)
base = y.copy()
vander = np.vander(x, order)
vander_pinv = LA.pinv(vander)
for _ in range(max_it):
coeffs_new = np.dot(vander_pinv, y)
if LA.norm(coeffs_new - coeffs) / LA.norm(coeffs) < tol:
break
coeffs = coeffs_new
base = np.dot(vander, coeffs)
y = np.minimum(y, base)
return base
def envelope(y, deg=None, max_it=None, tol=None):
"""
Computes the upper envelope of a given data.
It is implemented in terms of the `baseline` function.
Parameters
----------
y : ndarray
Data to detect the baseline.
deg : int
Degree of the polynomial that will estimate the envelope.
max_it : int
Maximum number of iterations to perform.
tol : float
Tolerance to use when comparing the difference between the current
fit coefficients and the ones from the last iteration.
Returns
-------
ndarray
Array with the envelope amplitude for every original point in *y*
"""
return y.max() - baseline(y.max() - y, deg, max_it, tol)
|
