noise_filter.java

java的小波分析程序

      x = pointz[i].val - mean;
      stdDevSum = stdDevSum + (x * x);
    }
    double sigmaSquared = stdDevSum / (size-1);
    double sigma = Math.sqrt( sigmaSquared );

    info.mean = mean;
    info.sigma = sigma;

    return pointz;
  } // alloc_points



  /**
    <p>
    normal_interval
    </p>

    <p>
    Numerically integreate the normal curve with mean
    <i>info.mean</i> and standard deviation <i>info.sigma</i>
    over the range <i>low</i> to <i>high</i>.
    </P>

    <p>
    There normal curve equation that is integrated is:
     </p>
     <pre>
       f(y) = (1/(s * sqrt(2 * pi)) e<sup>-(1/(2 * s<sup>2</sup>)(y-u)<sup>2</sup></sup>
     </pre>

     <p>
     Where <i>u</i> is the mean and <i>s</i> is the standard deviation.
     </p>

     <p>
     The area under the section of this curve from <i>low</i> to
     <i>high</i> is returned as the function result.
     </p>

     <p>
     The normal curve equation results in a curve expressed as
     a probability distribution, where probabilities are expressed
     as values greater than zero and less than one.  The total area
     under a normal curve with a mean of zero and a standard
     deviation of one is is one.
     </p>

     <p>
     The integral is calculated in a dumb fashion (e.g., we're not
     using anything fancy like simpson's rule).  The area in
     the interval <b>x</b><sub>i</sub> to <b>x</b><sub>i+1</sub>
     is 
     </P>

     <pre>
     area = (<b>x</b><sub>i+1</sub> - <b>x</b><sub>i</sub>) * g(<b>x</b><sub>i</sub>)
     </pre>

     <p>
     where the function g(<b>x</b><sub>i</sub>) is the point on the
     normal curve probability distribution at <b>x</b><sub>i</sub>.
     </p>

     @param info       This object encapsulates the mean and standard deviation
     @param low        Start of the integral
     @param high       End of the integral
     @param num_points Number of points to calculate (should be even)

   */
  private double normal_interval(bell_info info,
				 double low, 
				 double high, 
				 int num_points )
  {
    double integral = 0;

    if (info != null) {
      double s = info.sigma;
      // calculate 1/(s * sqrt(2 * pi)), where <i>s</i> is the stddev
      double sigmaSqrt = 1.0 / (s * (Math.sqrt(2 * Math.PI)));
      double oneOverTwoSigmaSqrd = 1.0 / (2 * s * s);

      double range = high - low;
      double step = range / num_points;
      double x = low;
      double f_of_x;
      double area;
      double t;
      for (int i = 0; i < num_points-1; i++) {
	t = x - info.mean;
	f_of_x = sigmaSqrt * Math.exp( -(oneOverTwoSigmaSqrd * t * t) );
	area = step * f_of_x; // area of one rectangle in the interval
	integral = integral + area;  // sum of the rectangles
	x = x + step;
      } // for
    }

    return integral;
  } // normal_interval


  /**
    <p>
    Set num_points values in the histogram bin <i>b</i>
    to zero.  Or, if the number of values is less
    than num_zero, set all values in the bin to zero.
    </p>

    <p>
    The num_zero argument is derived from the area under
    the normal curve in the histogram bin interval.  This
    area is a fraction of the total curve area.  When multiplied
    by the total number of coefficient points we get
    num_zero.
    </p>

    <p>
    The noise coefficients are preserved (returned) in the noise
    array argument.
    </p>

   */
  private void zero_points( bin b, int num_zero, double noise[] )
  {
    int num = b.vals.size();
    int end = num_zero;
    if (end > num)
      end = num;

    point p;
    for (int i = 0; i < end; i++) {
      p = (point)b.vals.elementAt( i );
      noise[ p.index ] = p.val;
      p.val = 0;
    }
  } // zero_points


  /**
    <p>
    Subtract the gaussian (or normal) curve from the histogram
    of the coefficients.  This is done by integrating the 
    gaussian curve over the range of a bin.  If the number of
    items in the bin is less than or equal to the area under the
    curve in that interval, all items in the bin are set to
    zero.  If the number of items in the bin is greater than
    the area under the curve, then a number of bin items equal
    to the curve area is set to zero.
    </p>
    <p>
    The area under a normal curve is always less than or equal
    to one.  So the area returned by normal_interval is the 
    fraction of the total area.  This is multiplied by
    the total number of coefficients.
    </p>
    <p>
    The function returns the number of coefficients that
    are set to zero (e.g., the number of coefficients that
    fell within the gaussian curve).  These coefficients are
    the noise coefficients.  The noise coefficients are returned
    in the noise argument.
    </p>
   */
  private int subtract_gauss_curve( bin binz[], 
			            bell_info info,
				    int total_points,
				    double noise[] )
  {
    int points_in_interval = total_points / binz.length;
    double start = binz[0].start;
    double end = binz[1].start;
    double step = end - start;
    double percent;
    int num_points;
    int total_zeroed = 0;

    for (int i = 0; i < binz.length; i++) {
      percent = normal_interval( info, start, end, points_in_interval );
      num_points = (int)(percent * (double)total_points);
      total_zeroed = total_zeroed + num_points;
      if (num_points > 0) {
	zero_points( binz[i], num_points, noise );
      }
      start = end;
      end = end + step;
    } // for

    return total_zeroed;
  } // subtract_gauss_curve

  
  /**
   <p>
   This function is passed the section of the Haar
   coefficients that correspond to a single spectrum.
   It compares this spectrum to a gaussian
   curve and zeros out the coefficients within the
   gaussian curve.
   </p>
   <p>
   The function returns the number of points filtered out as
   the function result.  The noise spectrum is also returned
   in the <i>noise</i> argument.
   </p>
   */
  private int filter_spectrum( double coef[], int start, int end,
				double noise[] )
  {
    final int num_bins = 32;
    int num_filtered;

    bell_info info = new bell_info();
    point pointz[] = alloc_points( coef, start, end, info );
    bin binz[] = calc_histo( pointz, num_bins );
    num_filtered = subtract_gauss_curve( binz, info, pointz.length, noise );
    
    int zero_count = 0;
    // copy filtered coefficients back into the coefficient array
    int ix = start;
    for (int i = 0; i < pointz.length; i++) {
      coef[ix] = pointz[i].val;
      ix++;
    }

    return num_filtered;
  } // filter_spectrum


  /**
    Normalize the noise array to zero by subtracting
    the smallest value from all points.
   */
  private void normalize_to_zero( double noise[] )
  {
    double min = noise[0];
    for (int i = 1; i < noise.length; i++) {
      if (min > noise[i])
	min = noise[i];
    } // for

    // normalize
    for (int i = 0; i < noise.length; i++) {
      noise[i] = noise[i] - min;
    }
  } // normalize_to_zero


  /**
   <p>
   This function is passed a set of Haar wavelet
   coefficients that result from the Haar wavelet
   transform.  It applies a gaussian noise filter
   to each frequency spectrum.  This filter zeros
   out coefficients that fall within a gaussian
   curve.  This alters the input data (the coef array).
   </p>
   <p>
   The <i>coef</i> argument is the input argument and
   contains the coefficients.  The <i>noise</i> argument
   is an output argument and contains the coefficients
   that have been filtered out.  This allows a noise
   spectrum to be rebuilt.
   </p>
   */
  public void gaussian_filter( double coef[], double noise[] )
  {
    final int min_size = 64;  // minimum spectrum size
    
    int total_filtered = 0;
    int num_filtered;
    int end = coef.length;
    int start = end >> 1;
    while (start >= min_size) {
      num_filtered = filter_spectrum( coef, start, end, noise );
      total_filtered = total_filtered + num_filtered;
      end = start;
      start = end >> 1;
    }

    // Note that coef[0] is the average across the
    // time series.  This value is needed to regenerate
    // the noise spectrum time series.
    noise[0] = coef[0];

    System.out.println("gaussian_filter: total points filtered out = " +
		       total_filtered );
  } // gaussian_filter



  /**
    Calculate the Haar tranform on the time series (whose
    length must be a factor of two) and filter it.  Then
    calculate the inverse transform and write the result
    to a file whose name is <i>file_name</i>.  A noise
    spectrum is written to <i>file_name</i>_noise.
   */
  public void filter_time_series( String file_name, double ts[] )
  {
    double noise[] = new double[ ts.length ];
    wavelets.inplace_haar haar = new wavelets.inplace_haar();
    haar.wavelet_calc( ts );
    haar.order();
    gaussian_filter( ts, noise );
    haar.inverse();

    PrintWriter prStr = OpenFile( file_name );
    if (prStr != null) {
      for (int i = 0; i < ts.length; i++) {
	prStr.println( i + " " + ts[i] );
      }
      prStr.close();
    }

    if (noise != null) {
      // calculate the inverse Haar function for the noise
      haar.setWavefx( noise );
      haar.setIsOrdered();
      haar.inverse();

      normalize_to_zero( noise );

      // write the noise spectrum out to a file
      prStr = OpenFile( file_name + "_noise" );
      if (prStr != null) {
	for (int i = 0; i < noise.length; i++) {
	  prStr.println( i + " " + noise[i] );
	}
	prStr.close();
      }
    }
  } // filter_time_series

} // noise_filter