001 package util.wavelet; 002 import java.util.List; 003 import static java.lang.Math.*; 004 005 /** An enumeration of noise computation methods for thresholding in wavelet space. 006 * Each entry in the enumeration implements a variation of the compute() method 007 * which takes as input the signal and its non-uniform noise map and returns 008 * its output in the noise component of the scales list. 009 * 010 * @author John Talbot 011 */ 012 public enum NoiseComputation { 013 UpperLimit() { 014 /** 015 * Compute the upper limit of the noise for each wavelet coefficient using 016 * a box car type range over the noise values affecting a particular wavelet 017 * coefficient and finding the maximim noise value in each range. 018 * The computed noise values are returned in the scales list. 019 * <b>Reference</b> : 020 * <a href="http://adsabs.harvard.edu/cgi-bin/nph-bib_query?bibcode=1997ApJ...482.1011S">Starck <i>et al.</i> : 1997</a>, 021 * <i>Astrophysics Journal</i> <b>482</b>, 1011. 022 * 023 * @param signal <b>input:</b> the signal with non-uniform noise 024 * @param scales <b>output:</b> noise values returned in noise component of each scale in list 025 */ 026 public void compute(Signal signal, java.util.List<Scale> scales) { 027 float[] sigma = getStandardDeviation(signal.size()); 028 for (int scaleNumber = 0; scaleNumber < scales.size(); scaleNumber++) { 029 Scale scale = scales.get(scaleNumber); 030 int step2 = (int)(pow(2.0d, (double) scaleNumber) + 0.5d) * 2; 031 for (int i = 0; i < signal.size(); i++) { 032 float maxNoise = 0.0f; 033 for (int j = i - step2; j <= i + step2; j++) maxNoise = max(maxNoise, signal.getNoise(j)); 034 scale.setNoise(i, sigma[scaleNumber] * maxNoise); 035 } 036 } 037 } 038 039 }, 040 Dirac() { 041 /** 042 * Compute the exact noise threshold values by convolving the square of the wavelet transform of the 043 * dirac delta function with the square of the non-uniform noise of a signal. 044 * (see Stark et al. 2002, Astronomical Image an data Analysis p.43 eq. 2.22) 045 * 046 * @param signal <b>input:</b> the signal with noise 047 * @param scales <b>output:</b> the pre-computed scales of the wavelet transform of the signal 048 */ 049 public void compute(Signal signal, java.util.List<Scale> scales) { 050 SimpleWaveletSpace diracWaveletSpace = new SimpleWaveletSpace(SimpleSignal.getDirac(signal.size())); 051 SimpleSignal noise2 = signal.getNoise().square(); 052 noise2.setBoundaryCondition(BoundaryCondition.Periodic); 053 for (Scale scale : scales) { 054 SimpleSignal result = new SimpleSignal(noise2); 055 SimpleSignal dirac2 = diracWaveletSpace.getScale(scale.getScaleNumber()).square(); 056 scale.setNoise(result.convolve(dirac2)); 057 } 058 /* TODO: check to see if the dirac function we are using works by using it to 059 * compute a simple wavelet transform 060 */ 061 } 062 }; 063 064 public static void main(String[] args) { 065 System.out.println("Testing the dirac method to compute noise thresholds ..."); 066 // TODO: 067 } 068 069 /** 070 * Compute the noise assuming unit variance for each point in the signal. 071 * 072 * @param signal <b>input:</b> the signal with noise 073 * @param scales <b>output:</b> the pre-computed scales of the wavelet transform of the signal 074 */ 075 public void compute(Signal signal, java.util.List<Scale> scales) { 076 float[] sigma = getStandardDeviation(signal.size()); 077 for (int scaleNumber = 0; scaleNumber < scales.size(); scaleNumber++) { 078 Scale scale = scales.get(scaleNumber); 079 int step2 = (int)(pow(2.0d, (double) scaleNumber) + 0.5d) * 2; 080 for (int i = 0; i < signal.size(); i++) { 081 float maxNoise = 0.0f; 082 for (int j = i - step2; j <= i + step2; j++) maxNoise = max(maxNoise, 1f); 083 scale.setNoise(i, sigma[scaleNumber] * maxNoise); 084 } 085 } 086 } 087 088 089 090 /** 091 * Get the standard deviation (sigma) at each scale in the wavelet transform 092 * of a random noise signal comprised of a normally distributed gaussian deviates. 093 * These sigmaNoise values can be used to estimate the sigma of a signal using : 094 * <ul> 095 * <li><b>sigma[scale] = sigma[0] * sigmaNoise[scale] / sigmaNoise[0]</b></li> 096 * </ul> 097 * where sigma[0] is estimated by equating it with the standard deviation of scale 0 098 * of the wavelet transform of a signal (i.e signalWavelet[0][i]). 099 * Then when the absolute value of signalWavelet[scale][i] is greater than 100 * 3*sigma[scale] we have a significant coefficient. 101 * If the signal to noise is non-uniform then another technique is used. 102 * (test data: when n=1200 sigma[] = {0.750, 0.304, 0.177, 0.112, 0.077, 0.062, 0.045, 0.030, 0.039)) 103 104 * @param size the size of the random signal 105 * @return a list of standard deviations at each scale of the wavelet transform of the noise signal 106 */ 107 public static float[] getStandardDeviation(int size) { 108 SimpleSignal noise = SimpleSignal.getGaussianNoise(size); 109 SimpleWaveletSpace simpleWaveletSpace = new SimpleWaveletSpace(noise); 110 int maximumScale = simpleWaveletSpace.getNumberOfScales(); 111 float[] sigma = new float[maximumScale]; 112 for (int scaleNumber = 0; scaleNumber < maximumScale; scaleNumber++) { 113 float sum = 0f; 114 float sum2 = 0f; 115 SimpleSignal scale = simpleWaveletSpace.getScale(scaleNumber); 116 for (int i = 0; i < size; i++) { 117 float coefficient = scale.getValue(i); 118 sum += coefficient; 119 sum2 += coefficient*coefficient; 120 } 121 sigma[scaleNumber] = (float) sqrt((sum2 - sum * sum / (float)size )/ ((float)size - 1f) ); 122 } 123 return sigma; 124 } 125 } 126 127 /* Wavelet transform of the dirac delta function : 128 129 16*h0 a-trous B3-spline smoothing kernel (h0) 1 4 6 4 1 130 c0 dirac delta funtion 0 0 1 0 0 131 16*c1 smoothed array 1 1 4 6 4 1 132 16*h1 smoothing kernel h1 1 0 4 0 6 0 4 0 1 133 256*c2 smoothed array 2 1 4 10 20 31 40 44 40 31 20 10 4 1 134 16*h2 1 0 0 0 4 0 0 0 6 0 0 0 4 0 0 0 1 135 4096*c3 1 4 10 20 35 56 84 120 161 204 246 284 315 336 344 336 315 284 246 204 161 120 84 56 35 20 10 4 1 136 -------------------------------------------------------------------------------------------------------------------- 137 16*w0 (c0-c1) wavelet scale 0 -1 -4 10 -4 -1 138 256*w1 (c1-c2) wavelet scale 1 -1 -4 -10 -20 -15 24 52 24 -15 -20 -10 -4 -1 139 w2 (c2-c3) wavelet scale 2 140 141 256 *w0^2 square of w0 1 16 100 16 1 142 65536*w1^2 square of w1 1 16 100 400 225 576 2704 576 225 400 100 16 1 143 w2^2 square of w2 144 */