A {@link MedianFilter} rewrites every sample, which throws away detail even where the signal + * was perfectly clean. The Hampel filter is the conservative version. For each incoming sample it + * looks at the window that ends at that sample and computes two robust statistics: + * + *
The MAD is rescaled by {@code 1.4826}, the factor that makes it match the standard deviation + * for normally distributed data, and the sample is declared an outlier when + * + *
+ * |x - median| > threshold * 1.4826 * MAD + *+ * + *
Outliers are reported and replaced by the window median; everything else passes through + * unchanged. Both statistics are robust, so the spike being tested cannot inflate the very + * yardstick it is measured against - which is exactly what happens if one uses a mean and a standard + * deviation instead. + * + *
A threshold of 3 is the usual starting point, but how eager the filter turns out to be depends + * just as much on the window, because the MAD of a handful of samples is itself a noisy estimate of + * the spread. Measured on clean Gaussian noise with a threshold of 3, the filter flags about 7% of + * the samples with a window of 5, 3% with a window of 11 and 0.5% with a window of 101, slowly + * approaching the 0.3% one would get from a perfect estimate of sigma. Widen the window, or raise + * the threshold, when too many good samples are being rewritten. + * + *
{@code
+ * HampelFilter filter = new HampelFilter(7, 3.0);
+ * for (double sample : sensorReadings) {
+ * double clean = filter.accept(sample);
+ * if (filter.lastWasOutlier()) {
+ * log.warn("spike of {} replaced by {}", sample, clean);
+ * }
+ * }
+ * }
+ *
+ * The window is causal and includes the sample being judged, so the filter introduces no delay, + * and the raw sample - not its replacement - is what enters the window, which keeps the statistics + * honest. Before the window has filled up the statistics are computed over the samples seen so far. + * If more than half of the window holds one and the same value the MAD is zero, and any sample that + * differs from the median at all is then flagged; that is the textbook behaviour of the identifier. + * + *
Each sample costs O(w) time for a window of {@code w} samples and nothing is allocated after + * construction. This class is not thread-safe. + * + * @see MedianFilter + * @see Median absolute deviation + */ +public final class HampelFilter { + + /** Makes the MAD of normally distributed data an unbiased estimator of its standard deviation. */ + public static final double GAUSSIAN_MAD_SCALE = 1.4826; + + /** Threshold, in robust standard deviations, used when none is given. */ + public static final double DEFAULT_THRESHOLD = 3.0; + + private final MedianFilter window; + private final double threshold; + private final double[] sortedWindow; + private final double[] deviations; + + private boolean lastWasOutlier; + private long outlierCount; + + /** + * Creates a filter with the customary threshold of three robust standard deviations. + * + * @param windowSize how many recent samples the statistics are computed over; odd sizes are the usual choice + * @throws IllegalArgumentException if {@code windowSize} is not positive + */ + public HampelFilter(int windowSize) { + this(windowSize, DEFAULT_THRESHOLD); + } + + /** + * Creates a filter. + * + * @param windowSize how many recent samples the statistics are computed over + * @param threshold how many robust standard deviations a sample may deviate before it counts as an + * outlier; must be non-negative + * @throws IllegalArgumentException if {@code windowSize} is not positive or {@code threshold} is negative or not finite + */ + public HampelFilter(int windowSize, double threshold) { + if (!(threshold >= 0.0) || !Double.isFinite(threshold)) { + throw new IllegalArgumentException("The threshold must be finite and non-negative, but was " + threshold); + } + this.window = new MedianFilter(windowSize); + this.threshold = threshold; + this.sortedWindow = new double[windowSize]; + this.deviations = new double[windowSize]; + } + + /** + * Feeds one sample into the filter. + * + * @param value the incoming sample + * @return the sample itself, or the window median if the sample was judged to be an outlier + * @throws IllegalArgumentException if {@code value} is NaN or infinite + */ + public double accept(double value) { + double median = window.accept(value); + double scaledDeviation = threshold * GAUSSIAN_MAD_SCALE * medianAbsoluteDeviation(); + lastWasOutlier = Math.abs(value - median) > scaledDeviation; + if (lastWasOutlier) { + outlierCount++; + return median; + } + return value; + } + + /** + * Filters a whole signal, one sample at a time, starting from the current state. + * + * @param signal the samples to filter + * @return a new array of the same length in which flagged samples are replaced by the local median + * @throws IllegalArgumentException if any sample is NaN or infinite + * @throws NullPointerException if {@code signal} is {@code null} + */ + public double[] filter(double[] signal) { + double[] filtered = new double[signal.length]; + for (int i = 0; i < signal.length; i++) { + filtered[i] = accept(signal[i]); + } + return filtered; + } + + /** + * Runs the filter over a signal and reports which samples were flagged, without altering them. + * + * @param signal the samples to inspect + * @return a new array of the same length, {@code true} where the sample was judged an outlier + * @throws IllegalArgumentException if any sample is NaN or infinite + * @throws NullPointerException if {@code signal} is {@code null} + */ + public boolean[] detectOutliers(double[] signal) { + boolean[] flags = new boolean[signal.length]; + for (int i = 0; i < signal.length; i++) { + accept(signal[i]); + flags[i] = lastWasOutlier; + } + return flags; + } + + /** + * Returns the median absolute deviation of the current window, the robust counterpart of the + * standard deviation. + * + * @return the MAD, {@code 0} when the window holds a single sample + * @throws IllegalStateException if no sample has been accepted yet + */ + public double medianAbsoluteDeviation() { + int size = window.copySortedWindow(sortedWindow); + if (size == 0) { + throw new IllegalStateException("The window is empty"); + } + double median = window.median(); + + // The window is sorted, so the absolute deviations form two already sorted runs that meet at + // the median: descending to its left, ascending to its right. Merging them is linear. + int split = size; + for (int i = 0; i < size; i++) { + if (sortedWindow[i] > median) { + split = i; + break; + } + } + int left = split - 1; + int right = split; + for (int i = 0; i < size; i++) { + boolean takeLeft = right >= size || (left >= 0 && median - sortedWindow[left] <= sortedWindow[right] - median); + if (takeLeft) { + deviations[i] = median - sortedWindow[left]; + left--; + } else { + deviations[i] = sortedWindow[right] - median; + right++; + } + } + + int middle = size / 2; + return size % 2 != 0 ? deviations[middle] : 0.5 * (deviations[middle - 1] + deviations[middle]); + } + + /** + * Returns the median of the current window. + * + * @return the window median + * @throws IllegalStateException if no sample has been accepted yet + */ + public double median() { + return window.median(); + } + + /** + * Tells whether the most recently accepted sample was flagged. + * + * @return {@code true} if the last sample was replaced by the median + */ + public boolean lastWasOutlier() { + return lastWasOutlier; + } + + /** + * Returns how many samples have been flagged since the last reset. + * + * @return the number of detected outliers + */ + public long outlierCount() { + return outlierCount; + } + + /** + * Returns the configured window length. + * + * @return the window size given at construction time + */ + public int windowSize() { + return window.windowSize(); + } + + /** + * Returns the configured threshold. + * + * @return the threshold in robust standard deviations + */ + public double threshold() { + return threshold; + } + + /** + * Returns how many samples the window currently holds. + * + * @return the current fill level, at most {@link #windowSize()} + */ + public int size() { + return window.size(); + } + + /** + * Empties the window and forgets the outlier count. + */ + public void reset() { + window.reset(); + lastWasOutlier = false; + outlierCount = 0; + } + + @Override + public String toString() { + return "HampelFilter{windowSize=" + windowSize() + ", threshold=" + threshold + ", outliers=" + outlierCount + '}'; + } +} diff --git a/src/main/java/com/thealgorithms/streaming/MedianFilter.java b/src/main/java/com/thealgorithms/streaming/MedianFilter.java new file mode 100644 index 000000000000..fa876a4feece --- /dev/null +++ b/src/main/java/com/thealgorithms/streaming/MedianFilter.java @@ -0,0 +1,229 @@ +package com.thealgorithms.streaming; + +import java.util.Arrays; + +/** + * A streaming median filter: it reports the median of the last {@code windowSize} samples of + * a signal. + * + *
The median filter is the standard first line of defence against impulsive noise. A single + * corrupted sample - a dropped bit, a spike on an analogue line, a sensor glitch - drags a moving + * average with it, but it cannot drag a median: as long as fewer than half of the samples in the + * window are corrupted, the output stays on the true signal. Unlike a linear low-pass filter, a + * median filter also preserves sharp edges instead of smearing them. + * + *
Two arrays of {@code windowSize} doubles are kept: the samples in arrival order, so the oldest + * one can be evicted, and the same samples in ascending order, so the median is the middle element. + * Insertion and removal in the sorted copy are binary search plus a memory move, which makes an + * update O(w) in the worst case but with a very small constant, and nothing is allocated after + * construction. + * + *
{@code
+ * MedianFilter filter = new MedianFilter(5);
+ * for (double sample : signal) {
+ * double clean = filter.accept(sample);
+ * }
+ * }
+ *
+ * The filter is causal: the value it reports for a sample is the median of that sample and the + * {@code windowSize - 1} preceding ones, so the output lags the input by roughly half a window. + * Before the window is full the median is taken over however many samples have arrived. With an even + * window size the median is the average of the two middle samples. + * + *
This class is not thread-safe. + * + * @see HampelFilter for a variant that only replaces samples identified as outliers + * @see Median filter + */ +public final class MedianFilter { + + private final double[] arrivalOrder; + private final double[] ascending; + private final int windowSize; + + private int size; + private int oldest; + + /** + * Creates an empty filter. + * + * @param windowSize how many recent samples the median is taken over; odd sizes are the usual choice + * @throws IllegalArgumentException if {@code windowSize} is not positive + */ + public MedianFilter(int windowSize) { + if (windowSize <= 0) { + throw new IllegalArgumentException("The window size must be positive, but was " + windowSize); + } + this.windowSize = windowSize; + this.arrivalOrder = new double[windowSize]; + this.ascending = new double[windowSize]; + } + + /** + * Feeds one sample into the filter. + * + * @param value the incoming sample + * @return the median of the window that now ends at this sample + * @throws IllegalArgumentException if {@code value} is NaN or infinite + */ + public double accept(double value) { + if (!Double.isFinite(value)) { + throw new IllegalArgumentException("Samples must be finite, but was " + value); + } + if (size == windowSize) { + removeFromSorted(arrivalOrder[oldest]); + size--; + } + insertIntoSorted(value); + arrivalOrder[oldest] = value; + oldest = (oldest + 1) % windowSize; + size++; + return median(); + } + + /** + * Filters a whole signal, one sample at a time, starting from the current state. + * + * @param signal the samples to filter + * @return a new array of the same length holding the running medians + * @throws IllegalArgumentException if any sample is NaN or infinite + * @throws NullPointerException if {@code signal} is {@code null} + */ + public double[] filter(double[] signal) { + double[] filtered = new double[signal.length]; + for (int i = 0; i < signal.length; i++) { + filtered[i] = accept(signal[i]); + } + return filtered; + } + + /** + * Returns the median of the samples currently in the window. + * + * @return the median, averaging the two middle samples when the window holds an even number of them + * @throws IllegalStateException if no sample has been accepted yet + */ + public double median() { + if (size == 0) { + throw new IllegalStateException("The window is empty"); + } + int middle = size / 2; + if (size % 2 != 0) { + return ascending[middle]; + } + return 0.5 * (ascending[middle - 1] + ascending[middle]); + } + + /** + * Returns the smallest sample currently in the window. + * + * @return the minimum of the window + * @throws IllegalStateException if no sample has been accepted yet + */ + public double min() { + requireNonEmpty(); + return ascending[0]; + } + + /** + * Returns the largest sample currently in the window. + * + * @return the maximum of the window + * @throws IllegalStateException if no sample has been accepted yet + */ + public double max() { + requireNonEmpty(); + return ascending[size - 1]; + } + + /** + * Returns the configured window length. + * + * @return the window size given at construction time + */ + public int windowSize() { + return windowSize; + } + + /** + * Returns how many samples the window currently holds. + * + * @return the current fill level, at most {@link #windowSize()} + */ + public int size() { + return size; + } + + /** + * Tells whether the window has warmed up. + * + * @return {@code true} once the window holds {@link #windowSize()} samples + */ + public boolean isFull() { + return size == windowSize; + } + + /** + * Tells whether the window holds no samples. + * + * @return {@code true} if nothing has been accepted since the last reset + */ + public boolean isEmpty() { + return size == 0; + } + + /** + * Returns the window contents in ascending order. + * + * @return a new array holding the sorted window + */ + public double[] sortedWindow() { + return Arrays.copyOf(ascending, size); + } + + /** + * Empties the window. + */ + public void reset() { + size = 0; + oldest = 0; + } + + @Override + public String toString() { + return "MedianFilter{windowSize=" + windowSize + ", size=" + size + ", median=" + (size == 0 ? Double.NaN : median()) + '}'; + } + + /** + * Copies the sorted window into a caller supplied buffer, which lets tight loops avoid allocating. + * + * @param destination buffer of at least {@link #windowSize()} elements + * @return the number of values written + */ + int copySortedWindow(double[] destination) { + System.arraycopy(ascending, 0, destination, 0, size); + return size; + } + + private void insertIntoSorted(double value) { + int position = Arrays.binarySearch(ascending, 0, size, value); + if (position < 0) { + position = -(position + 1); + } + System.arraycopy(ascending, position, ascending, position + 1, size - position); + ascending[position] = value; + } + + private void removeFromSorted(double value) { + int position = Arrays.binarySearch(ascending, 0, size, value); + System.arraycopy(ascending, position + 1, ascending, position, size - position - 1); + } + + private void requireNonEmpty() { + if (size == 0) { + throw new IllegalStateException("The window is empty"); + } + } +} diff --git a/src/test/java/com/thealgorithms/streaming/HampelFilterTest.java b/src/test/java/com/thealgorithms/streaming/HampelFilterTest.java new file mode 100644 index 000000000000..fe33656eb5da --- /dev/null +++ b/src/test/java/com/thealgorithms/streaming/HampelFilterTest.java @@ -0,0 +1,207 @@ +package com.thealgorithms.streaming; + +import static org.junit.jupiter.api.Assertions.assertEquals; +import static org.junit.jupiter.api.Assertions.assertFalse; +import static org.junit.jupiter.api.Assertions.assertThrows; +import static org.junit.jupiter.api.Assertions.assertTrue; + +import java.util.Arrays; +import java.util.Random; +import org.junit.jupiter.api.Assertions; +import org.junit.jupiter.api.DisplayName; +import org.junit.jupiter.api.Test; +import org.junit.jupiter.params.ParameterizedTest; +import org.junit.jupiter.params.provider.ValueSource; + +class HampelFilterTest { + + private static double median(double[] values) { + double[] sorted = values.clone(); + Arrays.sort(sorted); + int middle = sorted.length / 2; + return sorted.length % 2 != 0 ? sorted[middle] : 0.5 * (sorted[middle - 1] + sorted[middle]); + } + + /** + * Median absolute deviation of the window ending at {@code index}, computed the obvious way. + */ + private static double bruteForceMad(double[] signal, int index, int windowSize) { + int from = Math.max(0, index - windowSize + 1); + double[] window = Arrays.copyOfRange(signal, from, index + 1); + double windowMedian = median(window); + double[] deviations = new double[window.length]; + for (int i = 0; i < window.length; i++) { + deviations[i] = Math.abs(window[i] - windowMedian); + } + return median(deviations); + } + + @ParameterizedTest + @ValueSource(ints = {0, -1}) + void rejectsNonPositiveWindowSizes(int windowSize) { + assertThrows(IllegalArgumentException.class, () -> new HampelFilter(windowSize)); + } + + @ParameterizedTest + @ValueSource(doubles = {-1.0, Double.NaN, Double.POSITIVE_INFINITY}) + void rejectsInvalidThresholds(double threshold) { + assertThrows(IllegalArgumentException.class, () -> new HampelFilter(5, threshold)); + } + + @Test + void queriesBeforeTheFirstSampleFail() { + HampelFilter filter = new HampelFilter(5); + assertEquals(5, filter.windowSize()); + assertEquals(HampelFilter.DEFAULT_THRESHOLD, filter.threshold()); + assertEquals(0, filter.size()); + assertEquals(0L, filter.outlierCount()); + assertFalse(filter.lastWasOutlier()); + assertThrows(IllegalStateException.class, filter::median); + assertThrows(IllegalStateException.class, filter::medianAbsoluteDeviation); + } + + @ParameterizedTest + @ValueSource(doubles = {Double.NaN, Double.POSITIVE_INFINITY, Double.NEGATIVE_INFINITY}) + void rejectsNonFiniteSamples(double value) { + HampelFilter filter = new HampelFilter(5); + assertThrows(IllegalArgumentException.class, () -> filter.accept(value)); + } + + @ParameterizedTest + @ValueSource(ints = {1, 2, 3, 5, 8, 15}) + @DisplayName("the median absolute deviation agrees with a brute force computation") + void madMatchesBruteForce(int windowSize) { + Random random = new Random(20240517L + windowSize); + double[] signal = new double[2_000]; + for (int i = 0; i < signal.length; i++) { + signal[i] = random.nextInt(100) - 50; + } + + HampelFilter filter = new HampelFilter(windowSize, 3.0); + for (int i = 0; i < signal.length; i++) { + filter.accept(signal[i]); + assertEquals(bruteForceMad(signal, i, windowSize), filter.medianAbsoluteDeviation(), 1e-12, "at index " + i); + } + } + + @Test + @DisplayName("a single spike is flagged and replaced by the local median") + void replacesASpike() { + HampelFilter filter = new HampelFilter(5, 3.0); + double[] signal = {5.0, 5.1, 4.9, 5.0, 5.05, 50.0, 5.0, 5.1}; + double[] filtered = filter.filter(signal); + + assertEquals(1L, filter.outlierCount()); + assertTrue(filtered[5] < 6.0, "the spike should have been replaced, but stayed at " + filtered[5]); + for (int i = 0; i < signal.length; i++) { + if (i != 5) { + assertEquals(signal[i], filtered[i], 0.0, "clean sample at index " + i + " must pass through"); + } + } + } + + @Test + void detectOutliersReportsWithoutChangingTheSignal() { + HampelFilter filter = new HampelFilter(5, 3.0); + double[] signal = {5.0, 5.1, 4.9, 5.0, 5.05, 50.0, 5.0, 5.1}; + boolean[] flags = filter.detectOutliers(signal); + + boolean[] expected = new boolean[signal.length]; + expected[5] = true; + Assertions.assertArrayEquals(expected, flags); + assertEquals(1L, filter.outlierCount()); + } + + @Test + @DisplayName("a clean ramp contains no outliers, since the MAD grows with the slope") + void leavesASmoothRampAlone() { + HampelFilter filter = new HampelFilter(5, 3.0); + double[] ramp = new double[200]; + for (int i = 0; i < ramp.length; i++) { + ramp[i] = i; + } + Assertions.assertArrayEquals(ramp, filter.filter(ramp), 0.0); + assertEquals(0L, filter.outlierCount()); + } + + @Test + void leavesAConstantSignalAlone() { + HampelFilter filter = new HampelFilter(7); + for (int i = 0; i < 100; i++) { + assertEquals(3.0, filter.accept(3.0)); + } + assertEquals(0L, filter.outlierCount()); + assertEquals(0.0, filter.medianAbsoluteDeviation()); + } + + @Test + @DisplayName("only a few samples of clean Gaussian noise are flagged") + void staysQuietOnCleanNoise() { + Random random = new Random(4242L); + HampelFilter filter = new HampelFilter(11, 3.0); + for (int i = 0; i < 5_000; i++) { + filter.accept(random.nextGaussian()); + } + assertTrue(filter.outlierCount() < 500, "flagged " + filter.outlierCount() + " of 5000 clean samples"); + } + + @Test + @DisplayName("spikes buried in noise are caught") + void catchesSpikesInNoisyData() { + Random random = new Random(7L); + double[] signal = new double[1_000]; + for (int i = 0; i < signal.length; i++) { + signal[i] = 20.0 + random.nextGaussian(); + } + int[] spikes = {100, 300, 700}; + for (int spike : spikes) { + signal[spike] = 200.0; + } + + HampelFilter filter = new HampelFilter(9, 3.0); + double[] filtered = filter.filter(signal); + for (int spike : spikes) { + assertTrue(filtered[spike] < 30.0, "spike at " + spike + " survived as " + filtered[spike]); + } + } + + @Test + @DisplayName("a zero MAD makes the identifier maximally strict") + void flagsAnyDeviationWhenTheMadIsZero() { + HampelFilter filter = new HampelFilter(5, 3.0); + for (int i = 0; i < 5; i++) { + filter.accept(5.0); + } + assertEquals(0.0, filter.medianAbsoluteDeviation()); + assertEquals(5.0, filter.accept(5.5)); + assertTrue(filter.lastWasOutlier()); + } + + @Test + void aThresholdOfZeroFlagsEverythingOffTheMedian() { + HampelFilter filter = new HampelFilter(3, 0.0); + filter.accept(1.0); + filter.accept(2.0); + assertTrue(filter.lastWasOutlier()); + assertEquals(1.5, filter.median(), 1e-12); + } + + @Test + void resetForgetsTheWindowAndTheCounters() { + HampelFilter filter = new HampelFilter(5, 3.0); + filter.filter(new double[] {5.0, 5.1, 4.9, 5.0, 5.05, 50.0}); + assertEquals(1L, filter.outlierCount()); + + filter.reset(); + assertEquals(0, filter.size()); + assertEquals(0L, filter.outlierCount()); + assertFalse(filter.lastWasOutlier()); + assertThrows(IllegalStateException.class, filter::median); + } + + @Test + void toStringMentionsTheState() { + HampelFilter filter = new HampelFilter(5, 2.5); + assertTrue(filter.toString().contains("threshold=2.5"), filter.toString()); + } +} diff --git a/src/test/java/com/thealgorithms/streaming/MedianFilterTest.java b/src/test/java/com/thealgorithms/streaming/MedianFilterTest.java new file mode 100644 index 000000000000..a176000022cc --- /dev/null +++ b/src/test/java/com/thealgorithms/streaming/MedianFilterTest.java @@ -0,0 +1,172 @@ +package com.thealgorithms.streaming; + +import static org.junit.jupiter.api.Assertions.assertEquals; +import static org.junit.jupiter.api.Assertions.assertFalse; +import static org.junit.jupiter.api.Assertions.assertThrows; +import static org.junit.jupiter.api.Assertions.assertTrue; + +import java.util.Arrays; +import java.util.Random; +import org.junit.jupiter.api.Assertions; +import org.junit.jupiter.api.DisplayName; +import org.junit.jupiter.api.Test; +import org.junit.jupiter.params.ParameterizedTest; +import org.junit.jupiter.params.provider.ValueSource; + +class MedianFilterTest { + + /** + * Median of the last {@code windowSize} samples of {@code signal} ending at {@code index}, + * computed the obvious way. + */ + private static double bruteForceMedian(double[] signal, int index, int windowSize) { + int from = Math.max(0, index - windowSize + 1); + double[] window = Arrays.copyOfRange(signal, from, index + 1); + Arrays.sort(window); + int middle = window.length / 2; + return window.length % 2 != 0 ? window[middle] : 0.5 * (window[middle - 1] + window[middle]); + } + + @ParameterizedTest + @ValueSource(ints = {0, -1, Integer.MIN_VALUE}) + void rejectsNonPositiveWindowSizes(int windowSize) { + assertThrows(IllegalArgumentException.class, () -> new MedianFilter(windowSize)); + } + + @Test + void queriesBeforeTheFirstSampleFail() { + MedianFilter filter = new MedianFilter(3); + assertTrue(filter.isEmpty()); + assertFalse(filter.isFull()); + assertEquals(0, filter.size()); + assertEquals(3, filter.windowSize()); + assertEquals(0, filter.sortedWindow().length); + assertThrows(IllegalStateException.class, filter::median); + assertThrows(IllegalStateException.class, filter::min); + assertThrows(IllegalStateException.class, filter::max); + } + + @ParameterizedTest + @ValueSource(doubles = {Double.NaN, Double.POSITIVE_INFINITY, Double.NEGATIVE_INFINITY}) + void rejectsNonFiniteSamples(double value) { + MedianFilter filter = new MedianFilter(3); + assertThrows(IllegalArgumentException.class, () -> filter.accept(value)); + } + + @Test + @DisplayName("before the window fills up the median is taken over what has arrived") + void warmsUpGracefully() { + MedianFilter filter = new MedianFilter(5); + assertEquals(4.0, filter.accept(4.0)); + assertEquals(3.0, filter.accept(2.0), 1e-12); + assertEquals(4.0, filter.accept(9.0), 1e-12); + assertEquals(3, filter.size()); + assertFalse(filter.isFull()); + } + + @Test + void aWindowOfOneIsTheIdentity() { + MedianFilter filter = new MedianFilter(1); + assertEquals(5.0, filter.accept(5.0)); + assertEquals(-3.0, filter.accept(-3.0)); + assertTrue(filter.isFull()); + assertEquals(-3.0, filter.min()); + assertEquals(-3.0, filter.max()); + } + + @Test + @DisplayName("an even window averages the two middle samples") + void handlesEvenWindows() { + MedianFilter filter = new MedianFilter(4); + filter.accept(1.0); + filter.accept(2.0); + filter.accept(3.0); + assertEquals(2.5, filter.accept(4.0), 1e-12); + assertEquals(3.5, filter.accept(5.0), 1e-12); + } + + @Test + @DisplayName("a lone spike is removed, an edge is preserved") + void removesSpikesButKeepsEdges() { + MedianFilter filter = new MedianFilter(3); + double[] signal = {1.0, 1.0, 100.0, 1.0, 1.0, 9.0, 9.0, 9.0}; + double[] filtered = filter.filter(signal); + Assertions.assertArrayEquals(new double[] {1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 9.0, 9.0}, filtered, 1e-12); + } + + @Test + void reportsTheExtremesOfTheWindow() { + MedianFilter filter = new MedianFilter(3); + filter.filter(new double[] {5.0, 1.0, 9.0}); + assertEquals(1.0, filter.min()); + assertEquals(9.0, filter.max()); + Assertions.assertArrayEquals(new double[] {1.0, 5.0, 9.0}, filter.sortedWindow(), 0.0); + + filter.accept(7.0); + assertEquals(1.0, filter.min()); + assertEquals(9.0, filter.max()); + Assertions.assertArrayEquals(new double[] {1.0, 7.0, 9.0}, filter.sortedWindow(), 0.0); + } + + @ParameterizedTest + @ValueSource(ints = {1, 2, 3, 4, 7, 16, 31}) + @DisplayName("agrees with a brute force median for every window size") + void matchesBruteForce(int windowSize) { + Random random = new Random(20240517L + windowSize); + double[] signal = new double[3_000]; + for (int i = 0; i < signal.length; i++) { + signal[i] = random.nextInt(50); + } + + MedianFilter filter = new MedianFilter(windowSize); + for (int i = 0; i < signal.length; i++) { + assertEquals(bruteForceMedian(signal, i, windowSize), filter.accept(signal[i]), 1e-12, "at index " + i); + assertEquals(Math.min(i + 1, windowSize), filter.size()); + } + } + + @Test + @DisplayName("handles duplicates, which is where a sorted mirror of the window can go wrong") + void handlesHeavyDuplication() { + Random random = new Random(11L); + double[] signal = new double[2_000]; + for (int i = 0; i < signal.length; i++) { + signal[i] = random.nextInt(3); + } + + MedianFilter filter = new MedianFilter(5); + for (int i = 0; i < signal.length; i++) { + assertEquals(bruteForceMedian(signal, i, 5), filter.accept(signal[i]), 1e-12, "at index " + i); + } + } + + @Test + void resetEmptiesTheWindow() { + MedianFilter filter = new MedianFilter(3); + filter.filter(new double[] {1.0, 2.0, 3.0}); + filter.reset(); + + assertTrue(filter.isEmpty()); + assertThrows(IllegalStateException.class, filter::median); + assertEquals(8.0, filter.accept(8.0)); + } + + @Test + void toStringMentionsTheState() { + MedianFilter filter = new MedianFilter(3); + assertTrue(filter.toString().contains("size=0"), filter.toString()); + filter.accept(2.0); + assertTrue(filter.toString().contains("median=2.0"), filter.toString()); + } + + @Test + void reportsWhetherItHoldsSamples() { + MedianFilter filter = new MedianFilter(2); + assertTrue(filter.isEmpty()); + filter.accept(1.0); + assertFalse(filter.isEmpty()); + assertFalse(filter.isFull()); + filter.accept(2.0); + assertTrue(filter.isFull()); + } +}