Continuous Wavelet Transforms, a Java Implementation

/**
 * Author Mark Bishop; 2014
 * License GNU v3; 
 * This class is free software: you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with this program.  If not, see <http://www.gnu.org/licenses/>.
 */

import java.util.Random;

/**
 * Class responsibility: Provide examples of certain signals, some of which are
 * not analytic. All public members return xy as a double[][] such that xy[0] =
 * x and xy[1] = y. X is computed on [min:max] in increments dx = (max - min +
 * 1) / (count);
 */

public class Examples {

	/**
	 * An impulse at x = nearly 0, elsewhere 0. Min, max, and count must be chosen so
	 * that x = 0 +/- dt occurs in the sequence, eg., min = -8, max = 7, count = 16
	 * produces a sequence that includes x = 0 at x[8]
	 */
	public static double[][] delta(double min, double max, int count) {
		double dt = Signal.dxFromCountRange(max, min, count);
		double[][] xy = new double[2][count];
		double x = min;
		for (int i = 0; i < count; i++) {
			xy[0][i] = x;
			if (x<dt/2&&x>-dt/2) {
				xy[1][i] = 1.0;
			}
			x += dt;
		}
		return xy;
	}

	public static double[][] triangle(double t0, double A, double min,
			double max, int count) {
		double dx = Signal.dxFromCountRange(max, min, count);
		double[][] xy = new double[2][count];
		double slope = -A * A / (-2 * t0);
		double x = min;
		for (int i = 0; i < count; i++) {
			if (x > -2 * t0 && x < 2 * t0) {
				if (x < 0) {
					xy[1][i] = slope * x+A*A;
				} else {
					xy[1][i] = -slope * x+A*A;
				}
			} else {
				xy[1][i] = 0;
			}
			xy[0][i] = x;
			x += dx;
		}
		return xy;
	}
	
	public static double[][] box(double t0, double A, double min,
			double max, int count) {
		double dx = Signal.dxFromCountRange(max, min, count);
		double[][] xy = new double[2][count];		
		double x = min;
		for (int i = 0; i < count; i++) {
			if (x > -t0 && x < t0) {				
					xy[1][i] = A;
				
			} else {
				xy[1][i] = 0;
			}
			xy[0][i] = x;
			x += dx;
		}
		return xy;
	}


	public static double[][] chirp(double min, double max, int count) {
		double dx = Signal.dxFromCountRange(max, min, count);
		double[][] xy = new double[2][count];
		double x = min;
		for (int i = 0; i < count; i++) {
			xy[0][i] = x;
			xy[1][i] = Math.sin(10 * x * x);
			x += dx;
		}
		return xy;
	}

	/**
	 * To show the impulses; min, max, and count must be chosen such that there
	 * will exist x values where x[i] mod T = 0.
	 * 
	 * @param min
	 * @param max
	 * @param count
	 * @param T
	 *            the interval between impulses
	 * @return An impulse whenever x mod T = 0, elsewhere 0.
	 */
	public static double[][] shah(double min, double max, int count, int T) {
		double[][] xy = new double[2][count];
		double x = min;
		double dx = Signal.dxFromCountRange(max, min, count);
		for (int i = 0; i < count; i++) {
			xy[0][i] = x;
			if (Math.abs(xy[0][i] % T) < dx) {
				xy[1][i] = 1.0;
			}
			x += dx;
		}
		return xy;
	}

	/**
	 * Gaussian modified mixed sinusoids. Gaussian centered at 0 with stdev =
	 * 128, sinusoids = sin(x) +sin(x/64)
	 */
	public static double[][] gaussModSin(double min, double max, int count) {
		double[][] xy = new double[2][count];
		double mean = 0.0;
		double stdev = 128;
		String function = "(1/(pow(2*pi,0.5)*" + stdev + "))*exp(-pow((" + mean
				+ "-x),2)/(2*pow(" + stdev + ",2)))";
		xy = new Signal(min, max, count, function, Signal.RangeOption.MinMax).xy;
		double dx = Signal.dxFromCountRange(max, min, count);
		double x = min;
		for (int i = 0; i < count; i++) {
			double y = 0.2 * Math.sin(x) + Math.sin(x / 64);
			xy[1][i] *= y;
			x += dx;
		}
		return xy;
	}

	public static double[][] gaussian(double mean, double stdev, double min,
			double max, int count) {
		double[][] xy = new double[2][count];
		String function = "(1/(pow(2*pi,0.5)*" + stdev + "))*exp(-pow((" + mean
				+ "-x),2)/(2*pow(" + stdev + ",2)))";
		xy = new Signal(min, max, count, function, Signal.RangeOption.MinMax).xy;
		return xy;
	}

	public static double[][] gaussianDerivative(double mean, double stdev,
			double min, double max, int count, int m) {
		double[][] xy = new double[2][count];
		double[][] gauss = gaussian(mean, stdev, min, max, count);
		double[] x = gauss[0];
		int n = x.length;
		double f1 = Math.pow(-1, m + 1);
		double gamma = 1 / Gamma.gamma(m + .5);
		double hermite = 0;
		for (int i = 0; i < n; i++) {
			hermite = Hermite.probabilistHermitePoly(m, x[i] / (stdev));
			xy[1][i] = f1 * gamma * hermite * gauss[1][i] / Math.exp(m * m);

		}
		xy[0] = x;
		return xy;
	}

	/**
	 * n = 1024, x[0] = 0, x[n-1] = 40 Sinusoidal frequencies at 1.4, 1/2, and
	 * 1/9. However, but y is momentarily forced to zero at x[459], producing a
	 * singularity.
	 */
	public static double[][] signalDropout() {
		double min = 0;
		double max = 40;
		int count = 1024;
		double[] freq = new double[] { 1.4, 0.5, 1.0 / 3.0 };
		double[][] xy = new double[2][count];
		String function = "sin(2*pi*" + freq[0] + "*x)+sin(2*pi*" + freq[1]
				+ "*x)+0.3*sin(2*pi*" + freq[2] + "*x)";
		xy = new Signal(min, max, count, function, Signal.RangeOption.MinMax).xy;
		xy[1][459] = 0;
		return xy;
	}

	/**
	 * Devil's staircase
	 * 
	 * @param min
	 *            prefer 0
	 * @param max
	 *            prefer 1
	 * @param count
	 * @return Cantor function
	 */
	public static double[][] cantor(double min, double max, int count) {
		double[][] xy = new double[2][count];
		double dx = Signal.dxFromCountRange(max, min, count);
		double x = min;
		for (int i = 0; i < count; i++) {
			xy[0][i] = x;
			x += dx;
		}
		xy[1] = cantorFn(xy[0]);
		return xy;
	}

	/**
	 * Singular everywhere.
	 * 
	 * @param min
	 * @param max
	 * @param count
	 * @return a random sequence.
	 */
	public static double[][] random(double min, double max, int count) {
		Random r = new Random();
		double[][] xy = new double[2][count];
		double dx = Signal.dxFromCountRange(max, min, count);
		double x = min;
		for (int i = 0; i < count; i++) {
			xy[0][i] = x;
			x += dx;
			xy[1][i] = r.nextDouble();
		}
		return xy;
	}

	/**
	 * Notes for the following methods: In string representations of fractional
	 * values in bases other than 10, a colon (:) is used instead of a decimal
	 * point. All methods target positive bases only.
	 */

	/**
	 * Compute the Cantor function using a method of typographic substitution.
	 * 
	 * @param x
	 *            sequence for x axis
	 * @return
	 */
	private static double[] cantorFn(double[] x) {
		int n = x.length;
		double[] y = new double[n];
		String one = "1";
		String zero = "0";
		String point = ":";
		boolean isNegative = false;
		for (int i = 0; i < n; i++) {
			double base10 = x[i];
			// Express in base 3.
			String strB3 = base10ToBase(base10, 8, 3);
			if (strB3.contains("-")) {
				isNegative = true;
				strB3 = strB3.replaceFirst("-", "");
			} else {
				isNegative = false;
			}
			char[] charB3 = strB3.toCharArray();
			int l = charB3.length;
			boolean isFirst1 = true;
			for (int j = 0; j < l; j++) {
				// If it contains a 1, replace every digit after the first 1 by
				// 0.
				if (charB3[j] == one.charAt(0) && isFirst1) {
					isFirst1 = false;
					continue;
				} else if (!isFirst1 && charB3[j] != point.charAt(0)) {
					charB3[j] = zero.charAt(0);
				}
			}
			String strTransit = new String(charB3);
			// Replace all 2s with 1s
			strTransit = strTransit.replaceAll("2", "1");
			// Interpret the result as a binary number.
			y[i] = baseToBase10(strTransit, 2);
			if (isNegative) {
				y[i] *= -1.0;
			}
		}
		return y;
	}

	/**
	 * @return x Using the first "precision" places to the right of the decimal,
	 *         express x in base = newBase
	 */
	private static String base10ToBase(double x, int precision, int newBase) {
		String result = "";
		if (x < 0) {
			result += "-";
		}
		int characteristic = (int) x;
		char[] cCharacteristic = toBase((int) Math.abs(x), newBase);
		result += (new String(cCharacteristic) + ":");
		// discard digits below precision
		double mantissa = Math
				.abs((int) (Math.pow(10, precision) * (x - characteristic))
						/ Math.pow(10, precision));
		double interimProduct = mantissa;
		for (int i = 0; i < precision; i++) {
			double nextProduct = interimProduct * newBase;
			if ((int) nextProduct < 10) {
				result += (int) nextProduct;
			} else {
				result += Integer.toString((int) nextProduct, newBase);
			}
			interimProduct = nextProduct - (int) nextProduct;
		}
		return result.trim();
	}

	private static double baseToBase10(String strB2, int initialBase) {
		boolean isNegative = false;
		if (strB2.contains("-")) {
			isNegative = true;
			strB2 = strB2.replaceFirst("-", "");
		}
		double b10 = 0;
		String[] split = strB2.split(":");
		int characteristic = Integer.valueOf(split[0], initialBase);
		b10 += characteristic;
		if (split.length > 1) {
			char[] mantissa = split[1].toCharArray();
			int n = mantissa.length;
			for (int i = 0; i < n; i++) {
				double d = (double) (Character.getNumericValue(mantissa[i]));
				if (d != 0) {
					b10 += d * Math.pow(initialBase, -(i + 1));
				}
			}
		}
		if (isNegative) {
			b10 *= -1.0;
		}
		return b10;
	}

	/**
	 * Integers only.
	 * 
	 * @param x
	 *            an integer in base 10
	 * @param base
	 *            the base in which to re-express x
	 * @return x in base as a character array
	 */
	private static char[] toBase(int x, int base) {
		String strB = Integer.toString(x, base);
		char[] charB = strB.toCharArray();
		return charB;
	}

}