diff --git a/geomagio/algorithm/SqDistAlgorithm.py b/geomagio/algorithm/SqDistAlgorithm.py index 5ebca90dbe9dea8f332135b98bcea44cac380de0..2d0a36c91bbb88e9dbb117c3cc60ba106c353bcc 100644 --- a/geomagio/algorithm/SqDistAlgorithm.py +++ b/geomagio/algorithm/SqDistAlgorithm.py @@ -362,22 +362,21 @@ class SqDistAlgorithm(Algorithm): raise AlgorithmException( "sigma0 must have length %d" % (hstep + 1)) - # generate a vector of weights that will "smooth" seasonal # variations locally...the quotes are because what we really # do is distribute the error correction across a range of # seasonal corrections; we do NOT convovle this filter with # a signal to dampen noise, as is typical. -EJR 10/2016 - + # smooth parameter should specify the required cut-off period in terms # of discrete samples; for now, we generate a Gaussian filter according # to White et al. (USGS SIR 2014-5045). - fom = 10**(-3/20.) # halve power at corner frequency - omg = np.pi / np.float64(smooth) # corner angular frequency - sig = np.sqrt(-2*np.log(fom)/omg**2) + np.finfo(float).eps # prevent sig==0 - ts = np.linspace(np.max((-m, -3*np.round(sig))), - np.min((m, 3*np.round(sig))), - np.min((2*m, 6*np.round(sig))) + 1) # this is alwyays symmetric + fom = 10**(-3 / 20.) # halve power at corner frequency + omg = np.pi / np.float64(smooth) # corner angular frequency + sig = np.sqrt(-2 * np.log(fom) / omg**2) + np.finfo(float).eps # sig>0 + ts = np.linspace(np.max((-m, -3 * np.round(sig))), + np.min((m, 3 * np.round(sig))), + np.min((2 *m, 6 * np.round(sig))) + 1) nts = ts.size weights = np.exp(-0.5 * (ts / sig)**2) weights = weights / np.sum(weights) @@ -461,14 +460,18 @@ class SqDistAlgorithm(Algorithm): # performance. r[i + 1] = gamma * (1 - alpha) * et / m - ## update and append to s using equation-error formulation + # update and append to s using equation-error formulation #s[i + m] = s[i] + gamma * (1 - alpha) * et # distribute error correction across range of seasonal # corrections according to weights calculated above - s[i + m] = s[i] + gamma * (1-alpha) * et * weights[nts/2] - s[i+m-nts/2:i+m] = s[i+m-nts/2:i+m] + gamma * (1-alpha) * et * weights[:nts/2] - s[i+1:i+nts/2+1] = s[i+1:i+nts/2+1] + gamma * (1-alpha) * et * weights[nts/2+1:] + s[i + m] = s[i] + gamma * (1 - alpha) * et * weights[nts / 2] + s[i + m - nts / 2:i + m] = (s[i + m - nts / 2:i + m] + + gamma * (1 - alpha) * et * + weights[:nts / 2]) + s[i + 1:i + nts / 2 + 1] = (s[i + 1:i + nts / 2 + 1] + + gamma * (1 - alpha) * et * + weights[nts / 2 + 1:]) # update l and b using equation-error formulation l = l + phi * b + alpha * et