DP44/Common/DTS.Common.Utilities/NHTSASubSample.cs

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace DTS.Common.Utilities
{
    public class NHTSASubSample<T> : Exceptional
    {
        public UInt64 preTriggerSamples;
        public T[] data;
        public UInt16 subSampleInterval;
        public double sampleRate;
        public int iNpre;

        public void SubSample()
        {
            /*----------------------------------------------------------------------------

             * This is translated from Tim's C code translation.  Variable names are retained for tracability.
             * 
                This function is translated from the original NHTSA FORTRAN code by
                Tim Kippen of Diversified Technical Systems, Inc. 

                The original comments appear below:

                      SUBROUTINE SUBSMP (AOLD, DTOLD,DTNEW, NFP,NLP) 
            C**************************************************************************
            C.
            C.          This routine takes a signal contained in array AOLD,
            C.          which was sampled with a time step of DTOLD seconds, and 
            C.          converts it to a signal sampled at DTNEW seconds.  The revised
            C.          signal is returned in AOLD.  The routine uses an interpolation
            C.          algorithm which will set the sampling rate to exactly that 
            C.          specified by DTNEW, regardless of whether or not  DTNEW
            C.          is an integer multiple of DTOLD.
            C.          
            C.          This routine does not worry about aliasing, so the input
            C.          signal should have been appropiately filtered prior to 
            C.          invoking SUBSMP.  
            C.
            C.          Implemented on the NHTSA's Office of Crashworthiness
            C.          VAX 11/780.  The DO WHILE are not standard FORTRAN 77
            C.          an may have to be simulated with IF statements and GO TO's
            C.          on other machines
            C.
            C.          Start by setting the values at time 0 equal to each other
            C.          Then create the data with the new sampling rate before time 0
            C.          Will not create more than 50 data points before time 0.  
            C.
            C      Revision Record
            C      ---------------
            C.          Original   - Jeff Marcus (NHTSA)
            C.          July, 1987   Jeff Marcus (NHTSA)
            C.                           Revised to guarantee no shift of time 0
            C      Sept. 1987 - Clay Gabler
            C                   Revised to handle case of NFP=0
            C.
            C**************************************************************************
            */
            var pdOld = data;
            var dDeltaOld = 1.0 / sampleRate;
            var dtNew = 1.0 / sampleRate * subSampleInterval;
            iNpre = (int)preTriggerSamples;
            var iNpts = data.Length;


            double dTimeOld, dTimeNew, dTimeNext;
            int iNew, iOld;

            // start at t=0 (index is iPre)
            var iTZero = iNpre;
            var iNewPts = iNpts;

            // allocate array for new data points
            var pdNew = new T[1 + iNpts];

            var dtOld = dDeltaOld;

            // Anchor Signal at Time Zero
            pdNew[iTZero] = pdOld[iTZero];

            // Subsample Portion of Signal Prior to Time Zero
            if (iNpre > 0)
            {

                // start with first pre-trigger point index and calculate the time
                iNew = iTZero - 1;
                dTimeNew = (iNew - iTZero) * dtNew;

                // for every old point
                for (iOld = iTZero; iOld > 0; iOld--)
                {
                    // calculate the time for the old point and the next old point
                    dTimeOld = (iOld - iTZero) * dtOld;

                    if (iNew <= 0)
                        break;

                    dTimeNext = (iOld - iTZero - 1) * dtOld;

                    if (dTimeNew <= dTimeNext || iNew <= 0)
                    {
                        continue;
                    }

                    // Templates can't be used with math operators.  .NET 4.0 dynamic is a work around.
                    dynamic dOld = pdOld[iOld];
                    dynamic dOldMinusOne = pdOld[iOld - 1];

                    // calculate the new value
                    pdNew[iNew] = (T)(pdOld[iOld] -                      // current old value
                            ((dTimeNew - dTimeOld) / dtOld *     // delta Time
                            (dOldMinusOne - dOld))); // delta value for old data

                    // setup next new point and calculate the new deltaT
                    iNew--;
                    dTimeNew = (iNew - iTZero) * dtNew;
                }

                // calculate the new first point
                iNpre = iTZero - (iNew + 1);
            }

            if (iNpts > iNpre)
            {
                // start with first pre-trigger point index and calculate the time
                iNew = iTZero + 1;
                dTimeNew = (iNew - iTZero) * dtNew;

                // for every old point
                for (iOld = iTZero; iOld < iNpts; iOld++)
                {
                    // calculate the time for the old point and the next old point
                    dTimeOld = (iOld - iTZero) * dtOld;

                    if (iNew >= iNpts)
                        break;

                    dTimeNext = (iOld - iTZero + 1) * dtOld;

                    if (dTimeNew >= dTimeNext || iNew >= iNpts)
                    {
                        continue;
                    }

                    // Templates can't be used with math operators.  .NET 4.0 dynamic is a work around.
                    dynamic dOld = pdOld[iOld];
                    dynamic dOldMinusOne = pdOld[iOld - 1];

                    // calculate the new value
                    pdNew[iNew] = (T)(pdOld[iOld] -						 // current old value
                            ((dTimeNew - dTimeOld) / dtOld *	 // delta Time
                            (dOldMinusOne - dOld))); // delta value for old data

                    // setup next new point and calculate the new deltaT
                    iNew++;
                    dTimeNew = (iNew - iTZero) * dtNew;
                }

                // calculate the new number of points
                iNewPts = iNpre + (iNew - iTZero - 1);
            }

            // zero old array
            //memset( pdOld , 0, *iNpts*sizeof(double));

            // copy values for return to calling function
            iNpts = iNewPts;
            //*dDeltaOld = dtNew;

            data = new T[iNpts];
            for (var i = 0; i < iNpts; i++)
            {
                data[i] = pdNew[System.Math.Max(0, iTZero - (iNpre + 1) + i)];
            }

            sampleRate = 1.0 / dtNew;
            //free (pdNew);
            return;
        }
    }
}