codageneratortn.java

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

package org.trinet.jasi.coda.TN;

import java.util.*;

import org.trinet.jasi.coda.*;
import org.trinet.jasi.*;
import org.trinet.util.*;
import org.trinet.jasi.TN.*;
/** Base class for calculating time series coda data and get results in a CodaTN object.
* Concrete subclasses set parameters and override default methods.
*/
public abstract class CodaGeneratorTN implements CodaGeneratorIF {

    public static final String  LOW_AMP_CODA_DESC          = "lowStartAmp";
    public static final String  BELOW_CUTOFF_AMP_CODA_DESC = "cutoffAmp";
    public static final String  COUNT_CODA_DESC            = "count";
    public static final String  TRUNCATED_CODA_DESC        = "truncated";
    public static final String  NOISY_CODA_DESC            = "noisy";
    public static final String  HIGH_AMP_CODA_DESC         = "highAmp";

    public static final String  UNITS                      = "c";

    protected static boolean       debug = false;

    //protected Concat        fmt;                   // string buffer formater for output
    protected FilterIF      timeSeriesFilter;
    protected boolean       filterEnabled = false;

    protected String        algorithm;            // output to coda object
    protected int           goodWindowCount;      // output to coda object or total windows?
    protected double        tauFinal;             // output to coda object
    protected double        aveAbsAmpAtEndOfCoda; // output to coda object good or not windows
    protected double        qFree;                // output to coda object
    protected double        aFree;                // output to coda object
    protected double        aFix;                 // output to coda object
    protected double        rms;                  // output to coda object
    protected double        quality;              // output to coda object
    protected String        codaDesc;             // output to coda object
    protected String        codaPhase = "S";      // output to coda object

    protected int     firstGoodWindow;       // optional output
    protected double  firstWindowAveAbsAmp;  // optional output
    protected int     badWindowCount;        // optional output
    protected int     exitStatus;            // output via method
    protected String  exitStatusString;

    protected double[] goodWindowAmp;        // output to coda object
    protected double[] goodWindowTau;        // output to coda object

    protected double[] windowWeight;
    protected double[] log10WindowTime; // shared buffer used to cache window log10 time
    protected double[] log10WindowAmp;  // shared buffer used to cache window log10 amps

    protected int      noiseBiasSamples;
    protected int      maxGoodWindowCount;
    protected int      minGoodWindowsToEndCalcOnClipping;

    protected double   codaStartSignalToNoiseRatio;
    protected double   codaCutoffSignalToNoiseRatio;
    protected double   passThruNoiseDecayAmpRatio;

    protected double   qFix;

    protected boolean  resetOnClipping;

    protected double  windowSize;                // shared set from method input
    protected double  secsPerSample;             // shared set from method input
    protected double  secsPastPTimeAtStartingIdx;// shared set from method input
    protected int     lastSampleIdxInTimeSeries; // shared set from method input
    protected double  minCodaStartingAmp;        // shared set from method input
    protected double  maxCodaNoiseCutoffAmp;     // shared set from method input

    protected int     windowCount;               // shared calc counter could be output
    protected boolean isFirstWindow;             // shared calc flag
    protected double  aveAbsAmpOfLastGoodWindow; // shared calc value

    protected CodaGeneratorTN() {}

    public static void setDebug(boolean value) {
        debug = value;
    }

    protected void initBuffers(int maxCount) {
        this.maxGoodWindowCount = maxCount;
        windowWeight     = new double[maxCount];
        goodWindowTau    = new double[maxCount]; // output
        log10WindowTime  = new double[maxCount]; // shared buffer used to cache window log10 time
        goodWindowAmp    = new double[maxCount]; // output to coda object
        log10WindowAmp   = new double[maxCount]; // shared buffer used to cache window log10 amps
    }

    private void initOutputData() {

        codaDesc                  = "";
        qFree                     = 0.;
        aFree                     = 0.;
        aFix                      = 0.;
        rms                       = 0.;
        tauFinal                  = 0.;
        aveAbsAmpAtEndOfCoda      = 1.;
        firstWindowAveAbsAmp      = 1.;

        firstGoodWindow           = 0;
        badWindowCount            = 0;
        windowCount               = 0;
        goodWindowCount           = 0;
        quality                   = 0.;
        exitStatus                = 0;
    }

    protected TimeAmp [] getGoodWindowTimeAmps() {
        TimeAmp [] timeAmps = new TimeAmp [goodWindowCount];
        for (int idx = 0; idx < goodWindowCount; idx++) {
            timeAmps[idx] = new TimeAmp(goodWindowTau[idx], goodWindowAmp[idx]);
        }
        return timeAmps;
    }

    protected double getWindowStartingSampleSecsOffset(int windowCount) {
        return windowCount*getWindowSize();
    }

    protected int getWindowStartingSampleIndexOffset(int windowCount, double secsPerSample) {
        return (int) Math.round(getWindowStartingSampleSecsOffset(windowCount)/secsPerSample);
    }

// Coda object also needs to have its channel desc, datetime (Ptime) set by the invoking class
    private CodaTN resultsToCoda(CodaTN coda) {
        CodaTN resultsCoda = coda;
        if (resultsCoda == null) resultsCoda = new CodaTN();

        resultsCoda.changeDescriptor(new CodaPhaseDescriptorTN(codaPhase));

        resultsCoda.algorithm.setValue(algorithm);
        resultsCoda.ampUnits.setValue(UNITS);
        resultsCoda.windowSize.setValue(windowSize);
        resultsCoda.timeS.setValue(secsPastPTimeAtStartingIdx);

        resultsCoda.weightIn.setValue(quality);
        resultsCoda.tau.setValue(tauFinal);

        timeAmpsToResults(resultsCoda);

        if (goodWindowCount > 0) {
            resultsCoda.windowCount.setValue(goodWindowCount);
            resultsCoda.aFix.setValue(aFix);
            resultsCoda.aFree.setValue(aFree);
            resultsCoda.qFix.setValue(qFix);
            resultsCoda.qFree.setValue(qFree);
            resultsCoda.residual.setValue(rms);
         }

        return resultsCoda;
    }

    protected CodaTN timeAmpsToResults(CodaTN resultsCoda) {
        Collection timeAmpPairs = resultsCoda.getWindowTimeAmpPairs();
        if (debug) System.out.print("  CodaGeneratorTN.timeAmpsToResults goodWindowCount: " + goodWindowCount);
        if (goodWindowCount > 0) {
            TimeAmp [] timeAmps = getGoodWindowTimeAmps();
            int size = timeAmps.length;
            if (size > 6) size = 6;
            for (int idx = 0; idx < size; idx++) {
                timeAmpPairs.add(timeAmps[idx]);
                if (debug) System.out.print(" " + timeAmps[idx].toString() + " WT: " + windowWeight[idx]);
            }
        }
        else {
            if (debug) System.out.print(" zero good windows final timeamp:" + tauFinal + " " + aveAbsAmpAtEndOfCoda);
            timeAmpPairs.add(new TimeAmp(tauFinal, aveAbsAmpAtEndOfCoda));
        }
        if (debug) System.out.println();
        return resultsCoda;
    }

// Method to collect coda window data along time series;
// adapted from  CUSP CODAT_TM_R4_SUB.FOR
// Add method to use waveform object and Coda object???
    public boolean calcCoda(Coda result, double secsPerSample, double secsAfterPTime,
                            float [] timeSeriesAmps, int startingIdxOffset, double responseNormalizedCodaCutoffAmp,
                            double responseClippingAmp, double inLTA) {

        boolean status = calcCoda(secsPerSample, secsAfterPTime,
                            timeSeriesAmps, startingIdxOffset, responseNormalizedCodaCutoffAmp,
                            responseClippingAmp, inLTA);

        if (status) setCodaResults(result);
        return status;
    }

    public boolean calcCoda(double secsPerSample, double secsAfterPTime,
                            float [] timeSeriesAmps, int startingIdxOffset, double responseNormalizedCodaCutoffAmp,
                            double responseClippingAmp, double inLTA) {

        int numberOfSamples = timeSeriesAmps.length - startingIdxOffset;
        return calcCoda(secsPerSample, secsAfterPTime, timeSeriesAmps, startingIdxOffset, numberOfSamples,
                               responseNormalizedCodaCutoffAmp, responseClippingAmp, inLTA);
    }

    public boolean calcCoda( double secsPerSample, double secsAfterPTime,
                             float [] timeSeriesAmps, int startingIdxOffset, int timeSeriesSamples,
                             double responseNormalizedCodaCutoffAmp, double responseClippingAmp,
                             double inLTA
                           ) {

        initOutputData();

        if (timeSeriesSamples < 0 ) throw new IllegalArgumentException("input number of samples < 0");

        this.secsPerSample = secsPerSample;
        this.secsPastPTimeAtStartingIdx   = secsAfterPTime;

        int firstSampleIdxInTimeSeries = startingIdxOffset;

        if (timeSeriesSamples > timeSeriesAmps.length) timeSeriesSamples = timeSeriesAmps.length;
        if (firstSampleIdxInTimeSeries > 0 ) timeSeriesSamples -= firstSampleIdxInTimeSeries;
        if (timeSeriesSamples < 1) {
            System.out.println("  CodaGeneratorTN calcCoda input starting index beyond times series end.");
            return false;
        }
        this.lastSampleIdxInTimeSeries = firstSampleIdxInTimeSeries + timeSeriesSamples - 1;
        if (debug) System.out.println("  CodaGeneratorTN.calcCoda.secsPastPTimeAtStartingIdx: " +
            secsPastPTimeAtStartingIdx + " startingIdxOffset: " + startingIdxOffset +
            " lastSampleIdx: " + lastSampleIdxInTimeSeries);

        boolean hasResponseClippingAmp = false;
        if (responseClippingAmp > 0.) hasResponseClippingAmp = true;

        aveAbsAmpOfLastGoodWindow = 1.;
        isFirstWindow = true;

        int samplesPerWindow = (int) Math.rint(getWindowSize()/secsPerSample);

        int lowAmpWindowCount = 0;
        int validStartingAmpWindowCount = 0;       // count for amplitude check on first 2 windows

        int windowStartingSampleIdx = 0;
        int windowEndingSampleIdx = 0;
        int currentWindowSamples = 0;

        double averagedAbsAmp = 0.;

        float [] filteredAmps = filterTimeSeries(timeSeriesAmps);

        int nsamples = Math.min(timeSeriesSamples, noiseBiasSamples);
        double wfBias = calcBias(filteredAmps, 0, nsamples); // option change to pre-calculated LTA and bias for waveform input
        double longTermAverage =
             (inLTA == 0.) ? calcLongTermAbsAverage(filteredAmps, 0, nsamples, wfBias, startingIdxOffset): Math.abs(inLTA);


        minCodaStartingAmp = codaStartSignalToNoiseRatio * longTermAverage ;
        maxCodaNoiseCutoffAmp = codaCutoffSignalToNoiseRatio * longTermAverage;

        int iterations = (timeSeriesSamples + samplesPerWindow - 1)/samplesPerWindow;
        int iterCount = 0;

        if (debug) {
            System.out.println("  CodaGeneratorTN longTermAve: " + longTermAverage +
            " minCodaStartingAmp: " +  minCodaStartingAmp  + " maxCodaNoiseCutoffAmp: " + maxCodaNoiseCutoffAmp);
            System.out.println("  CodaGeneratorTN calcCoda bias:" + wfBias + " iterations:" + iterations);
            System.out.println("  CodaGeneratorTN calcCoda:" + inputToString());
        }
        MainLoop: while (iterCount < iterations) {
            iterCount++;
            CheckStatus: {   // break out of this block to check coda termination return conditions

                windowStartingSampleIdx = firstSampleIdxInTimeSeries +
                                              getWindowStartingSampleIndexOffset(windowCount, secsPerSample);
                windowEndingSampleIdx = windowStartingSampleIdx + samplesPerWindow - 1;       // end of current window

                if (windowStartingSampleIdx < 0) {