codageneratortn.java

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

                    System.out.println("  CodaGeneratorTN starting sample index < 0" );
                    windowCount++;
                    badWindowCount++;
                    break CheckStatus; // before start of input time series data do next window
                }

                if (windowStartingSampleIdx >= lastSampleIdxInTimeSeries) {               // end of data
                    if (debug) System.out.println("  CodaGeneratorTN starting sample index at end of series" );
                    currentWindowSamples = lastSampleIdxInTimeSeries + samplesPerWindow - windowStartingSampleIdx;
                    endOfCoda(windowEndingSampleIdx, currentWindowSamples, averagedAbsAmp);
                    return doCodaStatistics();
                }

                if (windowEndingSampleIdx > lastSampleIdxInTimeSeries) {
                    windowEndingSampleIdx = lastSampleIdxInTimeSeries; // a truncated window
                }

                currentWindowSamples = windowEndingSampleIdx - windowStartingSampleIdx + 1;   // length of current window
                if (2 * currentWindowSamples < samplesPerWindow) {  // don't calc ave amp just add this length to coda duration
                    //System.out.println("  CodaGeneratorTN window samples < window size" );
                    currentWindowSamples = lastSampleIdxInTimeSeries - windowStartingSampleIdx;
                    endOfCoda(windowEndingSampleIdx, currentWindowSamples, averagedAbsAmp);
                    return doCodaStatistics();
                }

                // Usual path for coda calcuation, determines averaged absolute amplitude of window samples;
                boolean hasNoClippedSamples = true;
                for (int idx = windowStartingSampleIdx; idx <= windowEndingSampleIdx; idx++) {
                    //if( hasResponseClippingAmp && filteredAmps[idx] >= responseClippingAmp) hasNoClippedSamples = false;
                    averagedAbsAmp += Math.abs((double) filteredAmps[idx] - wfBias);
                }
                averagedAbsAmp = averagedAbsAmp/currentWindowSamples;
                if( hasResponseClippingAmp && (averagedAbsAmp >= responseClippingAmp)) hasNoClippedSamples = false;

                if (debug) {
                    System.out.print("  CodaGeneratorTN.calcCoda windowStartIdx: " +
                        windowStartingSampleIdx + " endIdx: " + windowEndingSampleIdx);
                    System.out.println(" clipAmp: " + responseClippingAmp +  " " +
                        " unClipped: " +  hasNoClippedSamples + " aveAmp: " + averagedAbsAmp);
                }
                if (averagedAbsAmp < 1.0) averagedAbsAmp = 1.0;

                validStartingAmpWindowCount++;
                boolean hasLowAmplitude = false;               // not yet finished
                int startOfCodaWindowCount = lowAmpWindowCount;// low-amplitude for 1st two windows
                if (validStartingAmpWindowCount <= 2) {        // execute only on first 2 windows
                    if (validStartingAmpWindowCount == 1) firstWindowAveAbsAmp = averagedAbsAmp;
                    if (firstWindowAveAbsAmp < responseNormalizedCodaCutoffAmp) lowAmpWindowCount++;
                    if (averagedAbsAmp < minCodaStartingAmp) lowAmpWindowCount++;
                    if (lowAmpWindowCount == 1) {
                        aveAbsAmpAtEndOfCoda = Math.rint(minCodaStartingAmp);
                        quality = 0.1;
                    }
                    else aveAbsAmpAtEndOfCoda = averagedAbsAmp;
                    if (lowAmpWindowCount > startOfCodaWindowCount) hasLowAmplitude = true; // terminate and do coda statistics
                }

                if (lowAmpWindowCount > 1 && validStartingAmpWindowCount > 2) {
                    quality = 0.0;
                    exitStatus = -1;
                    exitStatusString = "lowStartingAmp";
                    codaDesc = LOW_AMP_CODA_DESC;
                    return false;
                }

                if (! hasLowAmplitude) lowAmpWindowCount = 0;

                if (lowAmpWindowCount > 0) {       // bad coda window
                    windowCount++;
                    badWindowCount++;
                    break CheckStatus;            // do next window
                }

                // normal "60mv" coda (1 window less than total windows calculated);
                if (averagedAbsAmp < responseNormalizedCodaCutoffAmp) {
                    tauFinal = getWindowTauSeconds(windowCount);
                    aveAbsAmpAtEndOfCoda = aveAbsAmpOfLastGoodWindow;
                    codaDesc = BELOW_CUTOFF_AMP_CODA_DESC;
                    exitStatusString = "responseCutoffAmp";
                    exitStatus = 3;
                    return doCodaStatistics();
                }

                // high noise coda, background noise greater than cutoff counts;
                else if (averagedAbsAmp < maxCodaNoiseCutoffAmp) {
                    tauFinal = getWindowTauSeconds(windowCount);
                    aveAbsAmpAtEndOfCoda = aveAbsAmpOfLastGoodWindow;
                    codaDesc = NOISY_CODA_DESC;
                    exitStatusString = "noiseCutoffAmp";
                    exitStatus = 4;
                    return doCodaStatistics();
                }

                // Check for clipped amp;
                if (goodWindowCount > 0) {        // check window average amp for clipping
                    if (hasNoClippedSamples && averagedAbsAmp < responseClippingAmp &&
                        averagedAbsAmp < aveAbsAmpOfLastGoodWindow*passThruNoiseDecayAmpRatio ) {
                        if (endOfCoda(windowEndingSampleIdx, currentWindowSamples, averagedAbsAmp)) return doCodaStatistics();
                        else break CheckStatus;
                    }
                    else {                       // reset processing of coda windows
                        codaDesc = HIGH_AMP_CODA_DESC;
                        if (goodWindowCount >= minGoodWindowsToEndCalcOnClipping) {   // end of processing
                            aveAbsAmpAtEndOfCoda = aveAbsAmpOfLastGoodWindow;
                            tauFinal = getWindowTauSeconds(windowCount);
                            exitStatusString = "minGoodClipped";
                            exitStatus = 1;
                            return doCodaStatistics();
                        }
                        else if (resetOnClipping) {                  // reset counts and start again - WHY should this be done?
                            badWindowCount += goodWindowCount + 1;
                            windowCount = badWindowCount;
                            goodWindowCount = 0;
                            isFirstWindow = true;
                            break CheckStatus; // do next window
                        }
                        else return doCodaStatistics();       // added this logic aww 9/00
                    }
                }
                else {   // no valid windows yet
                    if (hasNoClippedSamples && averagedAbsAmp < responseClippingAmp) {   // got a good one
                        if (endOfCoda(windowEndingSampleIdx, currentWindowSamples, averagedAbsAmp)) {
                            return doCodaStatistics();
                        }
                        else break CheckStatus;
                    }
                    else { // discard current coda windows;
                        windowCount++;
                        badWindowCount++;
                        break CheckStatus; // do next window
                    }
                }
            } // end of CheckStatus block

            if (goodWindowCount >= maxGoodWindowCount) {  // terminate prematurely may still have unanalyzed time series data
                    tauFinal = getWindowTauSeconds(windowCount);
                    aveAbsAmpAtEndOfCoda = aveAbsAmpOfLastGoodWindow;
                    codaDesc = COUNT_CODA_DESC;
                    exitStatusString = "maxGoodCount";
                    exitStatus = 5;
                    return doCodaStatistics();
            }

            if (windowEndingSampleIdx == lastSampleIdxInTimeSeries &&  windowCount > 0) { // series ended at last window done
                currentWindowSamples = lastSampleIdxInTimeSeries - windowStartingSampleIdx + 1;
                endOfCoda(windowEndingSampleIdx, currentWindowSamples, averagedAbsAmp);
                return doCodaStatistics();
            }


        } // end of iteration while loop
        return false;  // abnormal shouldn't get here, ran off end of data
    }

    public boolean isFilterEnabled() { return filterEnabled;}
    public boolean hasFilter() { return (timeSeriesFilter != null);}
    public void enableFilter() { filterEnabled = true; }
    public void disableFilter() { filterEnabled = false; }
    public void setFilter(String seedchan) {
        if (timeSeriesFilter != null) setFilter(timeSeriesFilter.getFilter(seedchan));
        else System.out.println("Warning : CodaGeneratorTN.setFilter(String) timeSeriesFilter is null");
    }

/** Sets the filter to apply to input waveform time series. */
    public void setFilter(FilterIF filter) {
        timeSeriesFilter = filter;
    }

    protected float [] filterTimeSeries(float [] timeSeriesAmps) { // de-spiking then bandpass would be nice
        if (timeSeriesFilter == null ) return timeSeriesAmps;
        return (filterEnabled) ? (float []) timeSeriesFilter.filter(timeSeriesAmps) : timeSeriesAmps;
    }

    protected double getWindowTauSeconds(int windowCount) {
        return secsPastPTimeAtStartingIdx + ((double) windowCount + .5)*windowSize;
    }

// Method invoked when time series scan is ended because of amplitude tests or end of data
    private boolean endOfCoda(int windowEndingSampleIdx, int currentWindowSamples, double averagedAbsAmp) {
        // Test for end of data coda time series scan
        if (windowEndingSampleIdx >= lastSampleIdxInTimeSeries || goodWindowCount >= maxGoodWindowCount) {
              tauFinal = secsPastPTimeAtStartingIdx
                         + getWindowStartingSampleSecsOffset(windowCount)
                         + currentWindowSamples*secsPerSample;
            aveAbsAmpAtEndOfCoda = aveAbsAmpOfLastGoodWindow;       // note scalar fudge factor below -aww
            if (windowEndingSampleIdx >= lastSampleIdxInTimeSeries) codaDesc = TRUNCATED_CODA_DESC;
            else codaDesc = COUNT_CODA_DESC;
            exitStatusString = "endOfSeries";
            exitStatus = 2;
            return true; // do coda statistics
        }

        // not at coda end thus set values for current iteration
        if (averagedAbsAmp <= 0.) {
            System.err.println("Coda error bad averagedAbsAmp:" + averagedAbsAmp + " ; reset to 1. count.");
            averagedAbsAmp = 1.;
        }
        aveAbsAmpOfLastGoodWindow = averagedAbsAmp;

        double tauCurrentWindow = getWindowTauSeconds(windowCount);
        if (tauCurrentWindow <= 0.) {
            System.err.println("Coda error bad tauCurrentWindow:"+ tauCurrentWindow+ " ; reset to 1.0 seconds.");
            tauCurrentWindow = 1.;
        }

        windowCount++;
        goodWindowCount++;

        saveBufferValues(tauCurrentWindow, averagedAbsAmp);

        if (isFirstWindow) {
            firstGoodWindow = windowCount;
            isFirstWindow = false;
        }
        return false;
    }

    protected void saveBufferValues(double tauCurrentWindow, double averagedAbsAmp) {
        int idx = goodWindowCount - 1;
        windowWeight[idx]    = getWindowWeight(tauCurrentWindow, averagedAbsAmp);
        log10WindowAmp[idx]  =  windowWeight[idx]*MathTN.log10(averagedAbsAmp);
        log10WindowTime[idx] = -windowWeight[idx]*MathTN.log10(tauCurrentWindow);
        goodWindowTau[idx]   = tauCurrentWindow;
        goodWindowAmp[idx]   = averagedAbsAmp;
        if (debug) System.out.println("  CodaGeneratorTN saveBufferValues fit goodWindowCount: " +
            goodWindowCount + " " + tauCurrentWindow + " " +  averagedAbsAmp + " Wt: " + windowWeight[idx]);
    }

    protected double getWindowWeight(double tauCurrentWindow, double averagedAbsAmp) {
        return 1.0;
    }

// Method is invoked at end coda data collection for each time series.
    private boolean doCodaStatistics() {

        tauFinal             = Math.abs(tauFinal);  // when is it less than zero if starting idx is negative?
        aveAbsAmpAtEndOfCoda = Math.rint(aveAbsAmpAtEndOfCoda);

        if (goodWindowCount < 1) {    // no valid coda windows
            exitStatusString = "noGoodWindows";
            exitStatus = -2;
            return false;
        }
        else if (qFix < .01) {       // why, if ever true?
            exitStatusString = "badQFix";
            exitStatus = -3;
            return false;
        }
        calcL1Norm();
        return true;
    }

// l1 norm fit of data to determine Q and A parameters
    protected void calcL1Norm() {
        double smallestResidualSum = 999999.;
        for (int idx = 0; idx < goodWindowCount; idx++) {  // solve for Afix using Qfix
            double Q0 = qFix;
            double A0 = (log10WindowAmp[idx] - log10WindowTime[idx] * Q0)/windowWeight[idx];
            double sum = 0.;

            for (int idx3 = 0; idx3 < goodWindowCount; idx3++) {
                double residual = windowWeight[idx3]*A0 + log10WindowTime[idx3]*Q0 - log10WindowAmp[idx3];
                sum += Math.abs(residual);
            }

            if (sum <= smallestResidualSum) {
                aFix = A0;
                smallestResidualSum = sum;
            }
        }