diff --git a/src/changes/changes.xml b/src/changes/changes.xml index 2ce9a18a3..adf9d6cb0 100644 --- a/src/changes/changes.xml +++ b/src/changes/changes.xml @@ -51,6 +51,10 @@ If the output is not quite correct, check for invisible trailing spaces! + + Added method "MathArrays#convolve(double[], double[])" to compute the + discrete, linear convolution of two sequences. + Added low-discrepancy random generator "HaltonSequenceGenerator". diff --git a/src/main/java/org/apache/commons/math3/util/MathArrays.java b/src/main/java/org/apache/commons/math3/util/MathArrays.java index 47ba15f98..ff87a7f2c 100644 --- a/src/main/java/org/apache/commons/math3/util/MathArrays.java +++ b/src/main/java/org/apache/commons/math3/util/MathArrays.java @@ -29,6 +29,7 @@ import org.apache.commons.math3.exception.DimensionMismatchException; import org.apache.commons.math3.exception.MathArithmeticException; import org.apache.commons.math3.exception.MathIllegalArgumentException; import org.apache.commons.math3.exception.MathInternalError; +import org.apache.commons.math3.exception.NoDataException; import org.apache.commons.math3.exception.NonMonotonicSequenceException; import org.apache.commons.math3.exception.NotPositiveException; import org.apache.commons.math3.exception.NotStrictlyPositiveException; @@ -1352,4 +1353,51 @@ public class MathArrays { return array; } + /** + * Calculates the convolution between two sequences. + *

+ * The solution is obtained via straightforward computation of the convolution sum (and not via FFT; for longer sequences, + * the performance of this method might be inferior to an FFT-based implementation). + * + * @param x the first sequence (double array of length {@code N}); the sequence is assumed to be zero elsewhere + * (i.e. {x[i]}=0 for i<0 and i>={@code N}). Typically, this sequence will represent an input signal to a system. + * @param h the second sequence (double array of length {@code M}); the sequence is assumed to be zero elsewhere + * (i.e. {h[i]}=0 for i<0 and i>={@code M}). Typically, this sequence will represent the impulse response of the system. + * @return the convolution of {@code x} and {@code h} (double array of length {@code N} + {@code M} -1) + * @throws NullArgumentException if either {@code x} or {@code h} is null + * @throws NoDataException if either {@code x} or {@code h} is empty + * + * @see Convolution (Wikipedia) + * @since 4.0 + */ + public static double[] convolve(double[] x, double[] h) throws NullArgumentException, NoDataException { + MathUtils.checkNotNull(x); + MathUtils.checkNotNull(h); + + final int N = x.length; + final int M = h.length; + + if (N == 0 || M == 0) { + throw new NoDataException(); + } + + // initialize the output array + final int totalLength = N + M - 1; + final double[] y = new double[totalLength]; + + // straightforward implementation of the convolution sum + for (int n = 0; n < totalLength; n++) { + double yn = 0; + for (int k = 0; k < M; k++) { + final int j = n - k; + if ((j > -1) && (j < N) ) { + yn = yn + x[j] * h[k]; + } + } + y[n] = yn; + } + + return y; + } + } diff --git a/src/test/java/org/apache/commons/math3/util/MathArraysTest.java b/src/test/java/org/apache/commons/math3/util/MathArraysTest.java index 30d247537..4c47f1d59 100644 --- a/src/test/java/org/apache/commons/math3/util/MathArraysTest.java +++ b/src/test/java/org/apache/commons/math3/util/MathArraysTest.java @@ -13,12 +13,15 @@ */ package org.apache.commons.math3.util; +import static org.junit.Assert.fail; + import java.util.Arrays; import org.apache.commons.math3.TestUtils; import org.apache.commons.math3.exception.DimensionMismatchException; import org.apache.commons.math3.exception.MathArithmeticException; import org.apache.commons.math3.exception.MathIllegalArgumentException; +import org.apache.commons.math3.exception.NoDataException; import org.apache.commons.math3.exception.NonMonotonicSequenceException; import org.apache.commons.math3.exception.NotPositiveException; import org.apache.commons.math3.exception.NotStrictlyPositiveException; @@ -833,4 +836,62 @@ public class MathArraysTest { Assert.fail("expecting MathIllegalArgumentException"); } catch (MathIllegalArgumentException ex) {} } + + @Test + public void testConvolve() { + /* Test Case (obtained via SciPy) + * x=[1.2,-1.8,1.4] + * h=[1,0.8,0.5,0.3] + * convolve(x,h) -> array([ 1.2 , -0.84, 0.56, 0.58, 0.16, 0.42]) + */ + double[] x1 = { 1.2, -1.8, 1.4 }; + double[] h1 = { 1, 0.8, 0.5, 0.3 }; + double[] y1 = { 1.2, -0.84, 0.56, 0.58, 0.16, 0.42 }; + double tolerance = 1e-13; + + double[] yActual = MathArrays.convolve(x1, h1); + Assert.assertArrayEquals(y1, yActual, tolerance); + + double[] x2 = { 1, 2, 3 }; + double[] h2 = { 0, 1, 0.5 }; + double[] y2 = { 0, 1, 2.5, 4, 1.5 }; + + yActual = MathArrays.convolve(x2, h2); + Assert.assertArrayEquals(y2, yActual, tolerance); + + try { + MathArrays.convolve(new double[]{1, 2}, null); + fail("an exception should have been thrown"); + } catch (NullArgumentException e) { + // expected behavior + } + + try { + MathArrays.convolve(null, new double[]{1, 2}); + fail("an exception should have been thrown"); + } catch (NullArgumentException e) { + // expected behavior + } + + try { + MathArrays.convolve(new double[]{1, 2}, new double[]{}); + fail("an exception should have been thrown"); + } catch (NoDataException e) { + // expected behavior + } + + try { + MathArrays.convolve(new double[]{}, new double[]{1, 2}); + fail("an exception should have been thrown"); + } catch (NoDataException e) { + // expected behavior + } + + try { + MathArrays.convolve(new double[]{}, new double[]{}); + fail("an exception should have been thrown"); + } catch (NoDataException e) { + // expected behavior + } + } }