12
12
13
13
from math import comb , factorial
14
14
from numbers import Integral , Real
15
- from typing import Any , Callable , Literal , Optional , Tuple , Union
15
+ from typing import Any , Callable , Dict , Literal , Optional , Tuple , Union
16
16
17
17
import numpy as np
18
18
from scipy .ndimage import median_filter
@@ -380,14 +380,15 @@ def estimate_noise_stddev(
380
380
window_size : Optional [int ] = None ,
381
381
extrapolator : Callable [..., np .ndarray ] = np .pad ,
382
382
extrapolator_args : Tuple [Any , ...] = ("reflect" ,),
383
+ extrapolator_kwargs : Optional [Dict [str , Any ]] = None ,
383
384
power : Literal [- 2 , - 1 , 1 , 2 ] = 1 ,
384
385
stddev_min : Union [float , int ] = 1e-10 ,
385
386
) -> np .ndarray :
386
387
"""
387
388
EXPERIMENTAL FEATURE
388
389
389
390
Estimates the local/global noise standard deviation of a series even in the presence
390
- of trends, like baselines and peaks, as well as outliers by using forward finite
391
+ of trends, like baselines and peaks, as well as outliers by using central finite
391
392
differences.
392
393
Please see the Notes section for further details.
393
394
@@ -405,6 +406,8 @@ def estimate_noise_stddev(
405
406
diff_accuracy : int, default=2
406
407
The accuracy of the finite difference approximation, which has to be an even
407
408
integer ``>= 2``.
409
+ Higher values will enhance the effect of outliers that will corrupt the noise
410
+ estimation of their neighborhood.
408
411
window_size : int or None, default=None
409
412
The odd window size around a datapoint to estimate its local noise standard
410
413
deviation.
@@ -423,7 +426,12 @@ def estimate_noise_stddev(
423
426
It has to be a callable with the following signature:
424
427
425
428
```python
426
- series_extrap = extrapolator(series, pad_width, *extrapolator_args)
429
+ series_extrap = extrapolator(
430
+ series,
431
+ pad_width,
432
+ *extrapolator_args,
433
+ **extrapolator_kwargs,
434
+ )
427
435
```
428
436
429
437
If ``window_size`` is ``None``, only the central finite differences kernel is
@@ -432,8 +440,12 @@ def estimate_noise_stddev(
432
440
side, but of course the quality of the noise estimation can be improved by using
433
441
a more sophisticated extrapolation method.
434
442
extrapolator_args : tuple, default=("reflect",)
435
- Additional arguments that are passed to the extrapolator function as described
436
- for ``extrapolator``.
443
+ Additional positional arguments that are passed to the extrapolator function as
444
+ described for ``extrapolator``.
445
+ extrapolator_kwargs : dict or None, default=None
446
+ Additional keyword arguments that are passed to the extrapolator function as
447
+ described for ``extrapolator``.
448
+ If ``None``, no additional keyword arguments are passed.
437
449
power : {-2, -1, 1, 2}, default=1
438
450
The power to which the noise standard deviation is raised.
439
451
This can be used to compute the:
@@ -447,6 +459,8 @@ def estimate_noise_stddev(
447
459
The minimum noise standard deviation that is allowed.
448
460
Any estimated noise standard deviation below this value will be set to this
449
461
value.
462
+ Borrowing an idea from image processing, the minimum noise standard deviation
463
+ can, e.g., be estimated from one or more feature-free regions of ``series``.
450
464
It must be at least ``1e-15``.
451
465
452
466
Returns
@@ -483,7 +497,17 @@ def estimate_noise_stddev(
483
497
applying a modified version of the Median Absolute Deviation (MAD) to the
484
498
derivative/differences of the signal. By using a moving MAD filter, the local noise
485
499
level can be estimated as well.
486
- The algorithms does not work well for signals that are perfectly noise-free.
500
+
501
+ From a workflow perspective, the following steps are performed on the signal:
502
+
503
+ - The signal is extrapolated to avoid edge effects.
504
+ - The central finite differences are computed.
505
+ - Their absolute values are taken.
506
+ - The median (global) or median filter (local) is applied to these absolute
507
+ differences. With proper scaling, this will give an estimate of the noise level.
508
+
509
+ There is one limitation, namely that the algorithm does not work well for signals
510
+ that are perfectly noise-free, but this is a rare case in practice.
487
511
488
512
The kernel size for the central finite difference kernel is given by
489
513
``2 * floor((differences + 1) / 2) - 1 + diff_accuracy``.
@@ -506,7 +530,7 @@ def estimate_noise_stddev(
506
530
)
507
531
if window_size % 2 == 0 :
508
532
raise ValueError (
509
- "Got window_size = {window_size}, expected an odd integer."
533
+ f "Got window_size = { window_size } , expected an odd integer."
510
534
)
511
535
512
536
# power
@@ -543,6 +567,13 @@ def estimate_noise_stddev(
543
567
"size)."
544
568
)
545
569
570
+ ### Preparation ###
571
+
572
+ # the keyword arguments for the extrapolator are set up
573
+ extrapolator_kwargs = (
574
+ extrapolator_kwargs if extrapolator_kwargs is not None else dict ()
575
+ )
576
+
546
577
### Noise Standard Deviation Estimation ###
547
578
548
579
# the signal is extrapolated to avoid edge effects
@@ -552,9 +583,10 @@ def estimate_noise_stddev(
552
583
series ,
553
584
pad_width ,
554
585
* extrapolator_args ,
586
+ ** extrapolator_kwargs ,
555
587
)
556
588
557
- # the absolute forward finite differences are computed ...
589
+ # the absolute central finite differences are computed ...
558
590
abs_diff_series = np .abs (
559
591
np .convolve (series_extrap , np .flip (diff_kernel ), mode = "valid" )
560
592
)
0 commit comments