Skip to content

Despike utils

goring_nikora 

goring_nikora(u, remaining_spikes=5, max_iter=10, robust_statistics=False)

Implements the Goring & Nikora (2002) phase-space de-spiking algorithm, returning modified velocity array.

Parameters:

Name Type Description Default
u ndarray

N-D array with ndim >= 1. The last axis is the time axis; all other axes are treated as independent rows.

 required
remaining_spikes int

Iterations will stop once there are remaining_spikes or fewer bad samples

 5
max_iter int

Maximum number of iterations

 10
robust_statistics bool

If True, ellipsoid centers will be based on the median and axis lengths will be based on median absolute deviation as suggested by Wahl (2003). If False, mean and standard deviation are used, consistent with the original Goring & Nikora implementation.

 False

Returns:

Type Description
ndarray

float64 array of the same shape as u with spikes removed and interpolated over.
References

Goring, D. G., & Nikora, V. I. (2002). Despiking acoustic Doppler velocimeter data. Journal of hydraulic engineering, 128(1), 117-126.

Wahl, T. L. (2003). Discussion of "Despiking acoustic doppler velocimeter data" by Derek G. Goring and Vladimir I. Nikora. Journal of Hydraulic Engineering, 129(6), 484-487.

Source code in src/pytoast/utils/despike_utils.py 
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
def goring_nikora(
    u: np.ndarray,
    remaining_spikes: int = 5,
    max_iter: int = 10,
    robust_statistics: bool = False,
) -> np.ndarray:
    """Implements the Goring & Nikora (2002) phase-space de-spiking algorithm, returning modified velocity array.

    Parameters
    ----------
    u : np.ndarray
        N-D array with `ndim >= 1`. The last axis is the time axis; all other
        axes are treated as independent rows.

    remaining_spikes : int
        Iterations will stop once there are `remaining_spikes` or fewer bad samples

    max_iter : int
        Maximum number of iterations

    robust_statistics : bool
        If True, ellipsoid centers will be based on the median and axis lengths will be based on median absolute
        deviation as suggested by Wahl (2003). If False, mean and standard deviation are used, consistent with the
        original Goring & Nikora implementation.

    Returns
    -------
    np.ndarray
        `float64` array of the same shape as `u` with spikes removed and interpolated over.

    References
    ----------
    Goring, D. G., & Nikora, V. I. (2002). Despiking acoustic Doppler velocimeter data. Journal of hydraulic
        engineering, 128(1), 117-126.

    Wahl, T. L. (2003). Discussion of "Despiking acoustic doppler velocimeter data" by
        Derek G. Goring and Vladimir I. Nikora. Journal of Hydraulic Engineering, 129(6), 484-487.
    """

    def flag_bad_indices(u: np.ndarray) -> np.ndarray:
        """Flag spikes in a 2D array (n_rows, n_samples) using phase-space
        method."""
        # Gradients along time axis
        du = np.gradient(u, axis=1) / 2
        du2 = np.gradient(du, axis=1) / 2

        # Per-row statistics
        if robust_statistics:
            sigma_u = median_abs_deviation(u, axis=1, nan_policy="omit")
            sigma_du = median_abs_deviation(du, axis=1, nan_policy="omit")
            sigma_du2 = median_abs_deviation(du2, axis=1, nan_policy="omit")
            u_bar = np.nanmedian(u, axis=1)
            du_bar = np.nanmedian(du, axis=1)
            du2_bar = np.nanmedian(du2, axis=1)
        else:
            sigma_u = np.nanstd(u, axis=1)
            sigma_du = np.nanstd(du, axis=1)
            sigma_du2 = np.nanstd(du2, axis=1)
            u_bar = np.nanmean(u, axis=1)
            du_bar = np.nanmean(du, axis=1)
            du2_bar = np.nanmean(du2, axis=1)

        # Expected absolute maximum
        n = u.shape[1]
        lam = np.sqrt(2 * np.log(n))

        # Rotation angle per row
        theta = np.arctan(np.nansum(u * du2, axis=1) / np.nansum(u**2, axis=1))

        # Ellipse axes (unrotated)
        a1 = lam * sigma_u
        b1 = lam * sigma_du
        a3 = lam * sigma_du
        b3 = lam * sigma_du2

        # Rotated ellipse axes via batched 2x2 solve
        cos_t = np.cos(theta)
        sin_t = np.sin(theta)
        A = np.empty((u.shape[0], 2, 2))
        A[:, 0, 0] = cos_t**2
        A[:, 0, 1] = sin_t**2
        A[:, 1, 0] = sin_t**2
        A[:, 1, 1] = cos_t**2
        b_vec = np.stack([(lam * sigma_u) ** 2, (lam * sigma_du2) ** 2], axis=1)  # (n_rows, 2)
        x = np.linalg.solve(A, b_vec[:, :, None]).squeeze(-1)  # (n_rows, 2)
        a2 = np.sqrt(x[:, 0])
        b2 = np.sqrt(x[:, 1])

        # Broadcast all (n_rows,) stats to (n_rows, 1)
        u_bar = u_bar[:, np.newaxis]
        du_bar = du_bar[:, np.newaxis]
        du2_bar = du2_bar[:, np.newaxis]
        a1 = a1[:, np.newaxis]
        b1 = b1[:, np.newaxis]
        a2 = a2[:, np.newaxis]
        b2 = b2[:, np.newaxis]
        a3 = a3[:, np.newaxis]
        b3 = b3[:, np.newaxis]
        cos_t = cos_t[:, np.newaxis]
        sin_t = sin_t[:, np.newaxis]

        # u vs du
        bad_u_du = (u - u_bar) ** 2 / a1**2 + (du - du_bar) ** 2 / b1**2 > 1

        # u vs du2 (rotated ellipse)
        bad_u_du2 = (
            (cos_t * (u - u_bar) + sin_t * (du2 - du2_bar)) ** 2 / a2**2
            + (sin_t * (u - u_bar) - cos_t * (du2 - du2_bar)) ** 2 / b2**2
        ) > 1

        # du vs du2
        bad_du_du2 = (du - du_bar) ** 2 / a3**2 + (du2 - du2_bar) ** 2 / b3**2 > 1

        return np.asarray(bad_u_du | bad_u_du2 | bad_du_du2)

    flat, orig_shape = _flatten_to_2d(u)
    bad_index = flag_bad_indices(flat)
    total_bad = np.sum(bad_index, axis=1)
    iterations = 0

    while np.any(total_bad > remaining_spikes) and iterations < max_iter:
        flat[bad_index] = np.nan
        flat = interp_rows(flat)
        bad_index = flag_bad_indices(flat)
        total_bad = np.sum(bad_index, axis=1)
        iterations += 1

    flat = interp_rows(flat)
    return flat.reshape(orig_shape)

recursive_gaussian 

recursive_gaussian(u, alpha=3.0, max_iter=10)

Despike each row of u by iteratively flagging samples beyond alpha * sigma of a Gaussian fit, replacing them with NaN, and interpolating along the time axis.

Parameters:

Name Type Description Default
u ndarray

N-D array with ndim >= 1. The last axis is the time axis.

 required
alpha float

Standard-deviation multiplier defining the spike threshold.

 3.0
max_iter int

Maximum number of refit/replace iterations per row.

 10

Returns:

Type Description
ndarray

float64 array of the same shape as u.
Source code in src/pytoast/utils/despike_utils.py 
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
def recursive_gaussian(
    u: np.ndarray,
    alpha: float = 3.0,
    max_iter: int = 10,
) -> np.ndarray:
    """Despike each row of `u` by iteratively flagging samples beyond
    `alpha * sigma` of a Gaussian fit, replacing them with NaN, and
    interpolating along the time axis.

    Parameters
    ----------
    u : np.ndarray
        N-D array with `ndim >= 1`. The last axis is the time axis.
    alpha : float, optional
        Standard-deviation multiplier defining the spike threshold.
    max_iter : int, optional
        Maximum number of refit/replace iterations per row.

    Returns
    -------
    np.ndarray
        `float64` array of the same shape as `u`.
    """
    flat, orig_shape = _flatten_to_2d(u)
    m, _ = flat.shape
    for ii in range(m):
        u_row = flat[ii, :]
        k = 0
        num_bad = np.inf
        while (k < max_iter) and (num_bad > 0):
            mean, std = norm.fit(u_row)
            bad_cols = (u_row < mean - alpha * std) | (u_row > mean + alpha * std)
            u_row[bad_cols] = np.nan
            num_bad = np.sum(bad_cols)
            u_row = naninterp(u_row)

        flat[ii, :] = u_row

    return np.asarray(flat.reshape(orig_shape))

threshold 

threshold(u, threshold_min=-3.0, threshold_max=3.0)

Threshold-based despiking.

Samples outside [threshold_min, threshold_max] are set to NaN and linearly interpolated over along the time axis.

Parameters:

Name Type Description Default
u ndarray

N-D array with ndim >= 1. The last axis is the time axis.

 required
threshold_min float

Value below which samples are flagged.

 -3.0
threshold_max float

Value above which samples are flagged.

 3.0

Returns:

Type Description
ndarray

float64 array of the same shape as u with spikes removed.
Source code in src/pytoast/utils/despike_utils.py 
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
def threshold(
    u: np.ndarray,
    threshold_min: float = -3.0,
    threshold_max: float = 3.0,
) -> np.ndarray:
    """Threshold-based despiking.

    Samples outside `[threshold_min, threshold_max]` are set to NaN and
    linearly interpolated over along the time axis.

    Parameters
    ----------
    u : np.ndarray
        N-D array with `ndim >= 1`. The last axis is the time axis.
    threshold_min : float, optional
        Value below which samples are flagged.
    threshold_max : float, optional
        Value above which samples are flagged.

    Returns
    -------
    np.ndarray
        `float64` array of the same shape as `u` with spikes removed.
    """
    flat, orig_shape = _flatten_to_2d(u)
    bad = (flat < threshold_min) | (flat > threshold_max)
    flat[bad] = np.nan
    flat = interp_rows(flat)
    return flat.reshape(orig_shape)