wafilter.java

一个用java写的地震分析软件(无源码)-used to write a seismic analysis software (without source)

package org.trinet.filters;

import java.util.*;
import javax.swing.*;
import java.awt.*;
import java.awt.event.*;

import org.trinet.jiggle.*;
import org.trinet.jasi.*;
import org.trinet.util.*;
/**
* Recursive filter to create Wood-Anderson waveform. Implimentation of the technique
* developed by Hiroo Kanamori described in "Continuous Monitoring of Ground-Motion
Parameters", BSSA, V89, pp. 311-316, Feb. 1999.<p>
*
* This is only valid for 20, 80, or 100 sample/sec data.
*/
public class WAFilter extends GenericFilter {

    /** samples per second */
    float sps;

    // WA consts set in setConstants()
    /** Wood-Anderson damping constant */
    float h;
    /** frequency constant*/
    float f0;
    /** Wood-Anderson gain */
    float gwa;
    /** Filter constant "c1" see equation (4) */
    float c1;
    /** Filter constant "c2" see equation (5) */
    float c2;

     public WAFilter() {
     }
     public WAFilter(int sampleRate) {
        this();
        setSampleRate(sampleRate);
     }

     public FilterIF getFilter(String type) {
	return new WAFilter();
     }
     public FilterIF getFilter(int sampleRate) {
	return new WAFilter(sampleRate);
     }
/*
     public FilterIF getFilter(Waveform wf) {
	return new WAFilter(wf.getSampleRate());
     }
     public WAFilter(double gainfactor) {
        super();

        setGainFactor(gainfactor);

     }
*/

/** Returns this same segment with the time-series Filtered. If the timeseries
 *  cannot be filtered (because it is at a sample rate this filter can't handle)
 *  the original WFSegment is return, unchanged. */
     public WFSegment filterSegment (WFSegment seg) {

     // set magic numbers, will return false is invalid sample rate
       if (setConstants(seg)) {

	 // note we DON'T call setSamples because we don't want to recalc bias, etc.
	   float ts[] = demean(seg.getTimeSeries());

	   // replace original with new filtered time series
	   // setTimeSeries() recalcs bias, max, min.
	   seg.setTimeSeries( filterTimeSeries(ts) );
	   seg.setIsFiltered(true);       // flag segment as filtered
       } else {

	   System.out.println("WAFilter: invalid sample rate = "+getSampleRate()+
		" for "+ seg.getChannelObj().toDelimitedString()+
		" Must be 20, 80 or 100 sps.");
       }
       return seg;
     }

     /** Set the units for the resulting waveform. */
     void setUnits(Waveform wf) {
       wf.setAmpUnits(Units.CM);         // velocity = cm/sec
     }
     private void setSampleRate(double sampleRate) {
          sps = (float) sampleRate;
     }
     /** Samples per second. */
     private float getSampleRate() {
        return sps;
     }
     /** Seconds per sample. */
     private float getSampleInterval () {
        return 1.0f/getSampleRate();
     }

     /** Set adjusted values of  h, f0, g, c1, c2 which depend on the sample rate
     * and gain. This is only valid for 20, 80, or 100 sample/sec data.
     * Returns false if sample rate is invalid.
     * These values are from Table 1: Kanamori, Maechling, and Hauksson, 1999. <p>
     *
     * h   = the damping constant (Uncorrected = 0.8) <br>
     * w0  = natural angular frequency (Uncorrected = 2(pi)/period or 2(pi) * freq. <br>
     * gwa = gain factor of a Wood-Anderson instrument (Uncorrected = 2800)<br>
     */
     boolean setConstants(WFSegment seg) {

          return setConstants(seg.getSampleInterval(), seg.getChannelObj().isHighGain());
     }
     /** Set adjusted values of  h, f0, g, c1, c2 which depend on the sample rate
     * and the gain (high or low). Sample interval is the sampling interval in
     * seconds per sample, e.g. 0.01 for 100 sample/sec data. The boolean argument is 'true'
     * for high-gains and false for low-gains. This is only valid for 20, 80, or
     * 100 sample/sec data. Returns false if sample rate is invalid.
     * These values are from Table 1: Kanamori, Maechling, and Hauksson, 1999. <p>
     *
     * h   = the damping constant (Uncorrected = 0.8) <br>
     * w0  = natural angular frequency (Uncorrected = 2(pi)/period or 2(pi) * freq. <br>
     * gwa = gain factor of a Wood-Anderson instrument (Uncorrected = 2800)<br>
     */
     boolean setConstants(double sampleInterval, boolean highGain) {

          setSampleRate(1.0f/sampleInterval);

	  // check for valid sample rate
	  if (getSampleRate() != 20.0f &&
	      getSampleRate() != 80.0f &&
	      getSampleRate() != 100.0f )  return false;


          // defaults
	   h   = 0.781f;
           f0  = 1.29f;
           gwa = 2963.0f;

          if (highGain) {
            if (getSampleRate() == 20.0f) {
               h   = 0.568f;
               f0  = 1.39f;
               gwa = 3110.0f;

            } else if (getSampleRate() == 80.0f || getSampleRate() == 100.0f) {
               h   = 0.781f;
               f0  = 1.29f;
               gwa = 2963.0f;
            }
          } else {                // low gain values
            if (getSampleRate() == 100.0f) {

               h   = 0.774f;
               f0  = 1.3f;
               gwa = 2999.0f;
            }
          }

          // natural freq. should = 2(pi)/0.8 sec^-1 OR 2(pi) * 1.25
          // the value of the f0 is the "adjusted" value to correct for
          // the recursive filter
          float w0 = (float) (2.0 * Math.PI * f0);   // natural angular frequency
          float wdt = w0 * getSampleInterval();
          c1 = 1.0f + (h*wdt);                       // equation (4)
          c2 = 1.0f + (2.0f * h*wdt) + (wdt*wdt);    // equation (5)

	  return true;
     }

/** Return recursively filtered timeseries given a timeseries.
 * If there are fewer than 2 points in the array the original array is returned.
 * The orginal is unchanged. */
       private float[] filterTimeSeries (float[] ts) {

         float wa[] = new float[ts.length];

         if (ts.length < 2) return ts;

	 float dt = getSampleInterval();

         float dtSquared = dt * dt;

         // init values
         wa[0]  = 0.0f;
         wa[1]  = 0.0f;

         float gdt = (float) ((gwa/gainFactor) * dt );

         for (int i = 2; i < ts.length -3; i++) {

             wa[i] = (float) ( gdt * (ts[i]-ts[i-1]) + (2.0 * c1 * wa[i-1]) - wa[i-2] ) / c2;

         }

         return wa;

       }

} // end of class