floatbutterworthhighpassfilter.java

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

package org.trinet.util;

/** Butterworth highpass filter. Only works (correctly) with 100, 80 and 40 sps data.
* Coefficients were calculated for 4 poles and a 1 Hz corner.
* A cosine taper is applied. */
public class FloatButterworthHighPassFilter implements FloatFilterIF {

    protected static boolean reverse = false;
    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];

/*
    private double gain;
    private double coef1;
    private double coef2;
    private double coef3;
    private double coef4;
*/
// default to values for 100 sps
    private double  gain  =  1.085574782;
    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 FloatButterworthHighPassFilter () { }

    public FloatButterworthHighPassFilter (boolean enableReverseFiltering) {
         this();
         reverse = enableReverseFiltering;
    }

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

    private FloatButterworthHighPassFilter (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);
        }
    }

    private FloatButterworthHighPassFilter (int sampleRate) {
       setCoefs(sampleRate);
    }

    /** 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 highpass filter with one corner at 1.0 Hz.
    */
    public void setCoefs(int sampleRate) {

        if (sampleRate == 80) {
// 4 pole butterworth hp 80 sps corner at 1 Hz
            gain  =  1.108101438;
            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) {
            gain  =  1.228117167;
            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
            gain  =  1.085574782;
            coef1 = -0.8485559993;
            coef2 =  3.5335352195;
            coef3 = -5.5208191366;
            coef4 =  3.8358255406;
        }
    }
    /** Return the correct filter for the channel. */
    public FilterIF getFilter(int sampleRate) {
        return new FloatButterworthHighPassFilter(sampleRate);
    }
    /** 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 FloatButterworthHighPassFilter("H");
        }
        else if (substr.equals("E")) {
             return new FloatButterworthHighPassFilter("E");
        } else {
             return null;
        }
    }

    /** Return a filtered time series. */
    public Object filter(Object in) {
       return filter((float []) in);
    }

    /** Filter the time-series. 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;
    }
    /** Remove bias (mean) from the time series. */
    public float [] demean(float [] values) {

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

        double sum = 0.;
        for (int idx = 0; idx < size; idx++) {
            sum += values[idx];
        }
        double mean = Math.round(sum/size);
        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];
    }
    /** If set true time-series will be reverse filtered to remove phase shift. */
    public static void setReverseFiltering(boolean value) {
        reverse = value;
    }

}