generalbutterworthfilter.java

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

package org.trinet.util;

import org.trinet.jasi.*;

/** Butterworth filter. Only works (correctly) with 100, 80 and 40 sps data.
* You may choose either Highpass or Lowpass at 1 Hz corner.
* You may choose to reverse filter or not.
* Coefficients were calculated for 4 poles and a 1 Hz corner.
* A cosine taper is applied. */
public class GeneralButterworthFilter implements FloatFilterIF, WaveformFilterIF {

    public static final int NZEROS = 4;
    public static final int NPOLES = 4;

    public static final double cosineTaperAlpha = 0.05;

    private double [] xv = new double [NZEROS+1];
    private double [] yv = new double [NPOLES+1];

    public static final int HIGHPASS = 0;
    public static final int LOWPASS  = 1;
//    public static final int BANDPASS = 2;

    protected static boolean reverse = false;
    protected static int type = HIGHPASS;
    protected static int sps = 100;

    private static final double  gainH100  =  1.085574782;
    private static final double  gainH80   =  1.108101438;
    private static final double  gainH40   =  1.228117167;

    private static final double  gainL100  =  1.112983215e+06;
    private static final double  gainL80   =  4.649932749e+05;
    private static final double  gainL40   =  3.201129162e+04;

// default to values for Highpass 100 sps
    private double  gain  =  gainH100;
    private double  coef1 = -0.8485559993;
    private double  coef2 =  3.5335352195;
    private double  coef3 = -5.5208191366;
    private double  coef4 =  3.8358255406;

//    protected static final FloatButterworthHighPassFilter H_SEEDCHAN_FILTER =
//                     new FloatButterworthHighPassFilter("H");
//    protected static final FloatButterworthHighPassFilter E_SEEDCHAN_FILTER =
//                     new FloatButterworthHighPassFilter("E");

    public GeneralButterworthFilter () { }

    public GeneralButterworthFilter (boolean enableReverseFiltering) {
         this(type, sps, reverse);
    }

    public GeneralButterworthFilter (int type, int sampleRate) {
       this(type, sampleRate, reverse);
    }
    public GeneralButterworthFilter (int type, boolean enableReverseFiltering) {
       this(type, sps, enableReverseFiltering);
    }
    public GeneralButterworthFilter (int type, int sampleRate, boolean enableReverseFiltering) {
       setType(type);        // must set type BEFORE coefs
       setCoefs(sampleRate);
       setReverseFiltering(enableReverseFiltering);
    }

    /** Filter (4 pole) coefficents are determined by the stream type.
    xx@deprecated use FloatButterworthHighPassFilter (int sampleRate)
    */

    public GeneralButterworthFilter (String seedchan) {
        String substr = seedchan.substring(0,1).toUpperCase();

        // assume all "H" and "E" are 100sps
        if (substr.equals("H") || substr.equals("E")) {
          setCoefs(100);
        }
        else if (substr.equals("B")) {
          setCoefs(40);
        }
    }

    /** Set filter type: HIGHPASS or LOWPASS, corner at 1Hz. */
    public void setType (int filterType) {
       if (filterType > 0 && filterType < 3) type = filterType;
    }

    public void setHighpass() {
      setType(HIGHPASS);
    }

    public void setLowpass() {
      setType(LOWPASS);
    }

    /** Set gain and filter coefficients for the given sample rate.
    Values are derived from filter design code of Fisher.
    See: http://www-users.cs.york.ac.uk/~fisher/mkfilter/trad.html<br>
    Fourth order filters with one corner at 1.0 Hz.  <br>
    If this is not called class defaults to 100 sps filter.
    */
    public void setCoefs(int sampleRate) {

        if (sampleRate == 80) {
// 4 pole butterworth hp 80 sps corner at 1 Hz
            coef1 = -0.8144059977;
            coef2 =  3.4247473473;
            coef3 = -5.4051668617;
            coef4 =  3.7947911031;
// 4 pole butterworth hp 40 sps corner at 1 Hz
        }
        else if (sampleRate == 40) {
            coef1 = -0.6630104844;
            coef2 =  2.9240526562;
            coef3 = -4.8512758825;
            coef4 =  3.5897338871;
        }
        else if (sampleRate == 100) {
// 4 pole butterworth hp 100 sps corner at 1 Hz
            coef1 = -0.8485559993;
            coef2 =  3.5335352195;
            coef3 = -5.5208191366;
            coef4 =  3.8358255406;
        }

        setGain(sampleRate);
    }
    /** Set gain for the given sample rate and filter type.
    Values are derived from filter design code of Fisher.
    See: http://www-users.cs.york.ac.uk/~fisher/mkfilter/trad.html<br>
    Fourth order filters with one corner at 1.0 Hz.
    */
    private void setGain (int sampleRate) {
       if (type == HIGHPASS) {
         if (sampleRate == 100) {
            gain = gainH100;
         } else if (sampleRate == 80) {
            gain = gainH80;
         } else if (sampleRate == 40) {
            gain = gainH40;
         }
       } else {
         if (sampleRate == 100) {
            gain = gainL100;
         } else if (sampleRate == 80) {
            gain = gainL80;
         } else if (sampleRate == 40) {
            gain = gainL40;
         }
       }
    }

    /** If set true time-series will be reverse filtered to remove phase shift. */
    public static void setReverseFiltering(boolean value) {
        reverse = value;
    }

    /** Return the correct highpass filter for the sample rate. */
    public FilterIF getFilter(int sampleRate) {
        return new GeneralButterworthFilter(HIGHPASS, sampleRate, reverse);
    }

    /** Return the correct filter for the sample rate. */
    public FilterIF getFilter(int type, int sampleRate) {
        return new GeneralButterworthFilter(type, sampleRate, reverse);
    }

    /** Return the correct filter for the channel.
    @deprecated Use getFilter (int sampleRate)*/
    public FilterIF getFilter(String seedchan) {
        String substr = seedchan.substring(0,1).toUpperCase();
//        if (substr.equals("H")) return H_SEEDCHAN_FILTER;
//        else if (substr.equals("E"))  return E_SEEDCHAN_FILTER;

        if (substr.equals("H")) {
             return new GeneralButterworthFilter("H");
        }
        else if (substr.equals("E")) {
             return new GeneralButterworthFilter("E");
        }
        else if (substr.equals("B")) {
             return new GeneralButterworthFilter("B");
        } else {
             return null;
        }
    }

    /** Returns passed object with timeseries filtered in place.
     *  If the object passed is not float[] or Waveform the original object is
     *  returned unchanged. */
    public Object filter(Object in) {
       if (in instanceof float[]) return filter((float []) in);
       if (in instanceof Waveform) return filter((Waveform)in);

       return in;
    }

    /** Filter the time-series in place . Applies a cosine taper. Will reverse filter
    if setReverse(true). */
    public float[] filter(float [] in) {
       float [] values = demean(in);
       values = org.trinet.util.AmplitudeCosineTaper.taper(values, cosineTaperAlpha);
       values = doFilter(values);
       return (reverse) ? doReverseFilter(values) : values;
    }

    /** Returns Waveform with time series filtered. Note that the original timeseries
     *  in the Waveform is altered so make a copy if you want to retain original.
     * Returns 'null' if filtering fails.  */
     public Waveform filter (Waveform wf) {

       WFSegment seg[] = wf.getArray();
       if (seg != null) {
	    for (int i=0; i < seg.length; i++) {
	      filter(seg[i].getTimeSeries());  // filter the float[] in place
	      seg[i].scanTimeSeries();
	      seg[i].setIsFiltered(true);
	      wf.setIsFiltered(true);
	    }
	}

       // Recalculate bias, max, min, etc.
       wf.scanForBias();
       wf.scanForAmps();

       return wf;
     }


    /** Remove bias (mean) from the time series in place. */
    public float [] demean(float [] values) {

        if (values == null) return values;
        int size = values.length;
        if (size == 0) return values;

	// calculate mean
        double sum = 0.;
        for (int idx = 0; idx < size; idx++) {
            sum += values[idx];
        }
        double mean = Math.round(sum/size);

	// remove mean
        for (int idx = 0; idx < size; idx++) {
            values[idx] = (float) (values[idx] - mean);
        }
        return values;

    }

    private void zeroBuffers() {
        java.util.Arrays.fill(xv, 0.);
        java.util.Arrays.fill(yv, 0.);
    }

    /** Forward filter. */
    private float [] doFilter(float [] in) {
        if (in == null || gain == 0.) return in;
        zeroBuffers();
        int size = in.length;
        for (int idx=0; idx < size; idx++) {
            in[idx] = filterSample(in[idx]);
        }
        return in;
    }

    /** Reverse filter to undo phase shift. */
    private float [] doReverseFilter(float [] in) {
        if (in == null || gain == 0.) return in;
        zeroBuffers();
        int size = in.length;
        for (int idx=size-1; idx > -1; idx--) {
            in[idx] = filterSample(in[idx]);
        }
        return in;
    }

    /** Filter a single sample in the time series. */
    private float filterSample(float value) {
        double inputValue = value;
        xv[0] = xv[1];
        xv[1] = xv[2];
        xv[2] = xv[3];
        xv[3] = xv[4];
        xv[4] = inputValue / gain;
        yv[0] = yv[1];
        yv[1] = yv[2];
        yv[2] = yv[3];
        yv[3] = yv[4];
        yv[4] = (xv[0] + xv[4]) - 4 * (xv[1] + xv[3]) + 6 * xv[2]
                 + ( coef1 * yv[0]) + (  coef2 * yv[1])
                 + ( coef3 * yv[2]) + (  coef4 * yv[3]);
        return  (float) yv[4];
    }


}