# Patent

## Title
Resource adaptive spectrum estimation of streaming data

## Abstract

Introduced herein is a “resource adaptive” algorithm for spectrum and periodicity estimation on a continuous stream of data. The formulation is based on the derivation of a closed-form incremental computation of the spectrum, augmented by a load-shedding scheme that can adapt to available CPU resources to provide a resource efficient solution for real time spectrum estimation.

## Background
This invention was made with Government support under Contract No.: H98230-05-3-0001 awarded by the U.S. Department of Defense. The Government has certain rights in this invention.


FIELD OF THE INVENTION
The present invention relates to spectrum or periodogram estimation in streaming data under conditions of limited resources.
BACKGROUND OF THE INVENTION
Spectrum estimation, that is, analysis of the frequency content of a signal, is a core operation in numerous applications, such as data compression, medical data analysis (ECG data), pitch detection of musical content, and other applications. Widely used estimators of the frequency content are the periodogram and the autocorrelation of a sequence. For statically stored sequences, both methods have an O(nlogn) complexity using the Fast Fourier Transform (FFT). For dynamically updated sequences (streaming case), the same estimators can be computed incrementally, by continuous update of the summation in the FFT computation, through the use of Momentary Fourier Transform.
However, in a high-rate, data streaming environment with multiple processes ‘competing’ over computational resources, there is no guarantee that each running process will be allotted sufficient processing time to fully complete its operation. Instead of blocking or abandoning the execution of processing threads that cannot fully complete, a desirable compromise would be for the system to make provisions for adaptive process computation. Under this processing model every analytic unit (e.g., in this case the ‘periodogram estimation unit’) can provide partial (coarser') results under tight processing constraints.
Under the aforementioned processing model and given limited processing time, one may not be seeking for results that are accurate or perfect, but only ‘good enough’. Even so, since a typical streaming application will require fast, ‘on-the-fly’ decisions, an intelligent sampling procedure of exemplary efficiency would appear to represent a significant improvement over conventional efforts. A need has thus been recognized in connection with effecting such an improvement, among others.
SUMMARY OF THE INVENTION
There is broadly contemplated herein a method and apparatus for periodogram estimation based on resource (such as CPU, memory etc.) availability, in accordance with at least one presently preferred embodiment of the present invention. Also broadly contemplated herein is an intelligent sampling procedure that can decide whether to retain or discard an examined sample, based on a “lightweight” linear predictor whereby a sample is recorded only if its value cannot be predicted by previously seen sequence values.
Also, considering that in view of the sampling process, the retained data samples (a subset of the examined data window) are not guaranteed to be equi-spaced, there is also contemplated herein an elaboration of a closed-form periodogram estimation in the context of unevenly spaced samples.
In summary, one aspect of the invention provides a method of providing a spectrum estimation for data in a data stream, the method comprising the step of providing a spectrum estimation based on resource availability.
Another aspect of the invention provides an apparatus for providing a spectrum estimation for data in a data stream, the apparatus comprising an arrangement for providing a spectrum estimation based on resource availability.
Furthermore, an additional aspect of the present invention provides a program storage device readable by machine, tangibly embodying a program of instructions executable by the machine to perform method steps for providing a spectrum estimation for data in a data stream, the method comprising the step of providing a spectrum estimation based on resource availability.



BRIEF DESCRIPTION OF THE DRAWINGS
The above and other objects, features, and advantages of the present invention will become more apparent from the following detailed description when taken in conjunction with the accompanying drawings, in which:
FIG. 1 provides a general visual depiction of a presently inventive methodology.
FIG. 2 illustrates an unevenly sampled signal.
FIGS. 3A-3B illustrate a linear interpolation scheme for the adaptive pruning of samples.
FIGS. 4A-4B illustrates a comparison of spectrum estimation errors for intelligent sampling vs. equi-sampling techniques.
FIGS. 5A-5H provide a comparison of estimations for additional datasets.
FIGS. 6A-6F illustrate spectrum approximation for different threshold values, for weblog data.
FIGS. 7A-7F illustrate spectrum approximation for different threshold values, for chaotic data.
FIG. 8 illustrates a training phase for a threshold estimator.
FIGS. 9A-9H provide a comparison of spectrum estimations for different compression rates.
FIGS. 10A-10F provide a comparison of sampling techniques for different compression rates.
FIGS. 11A-11C illustrate experimental results for a presently inventive algorithm as applied to streaming automotive measurements.
FIGS. 12A-12B provide comparative prediction error histograms.



DETAILED DESCRIPTION OF THE INVENTION
Generally speaking, in considering a data streaming scenario, one objective addressed herein is the provision of efficient mechanisms for estimating and updating the spectrum within a current data window. As such, a periodogram may be used as an estimate of the spectrum. A schematic illustration of resource-adaptive methodology in accordance with a preferred embodiment of the present invention is provided in FIG. 1. (It should be noted that throughout the present disclosure, the terms “periodicity estimation”, “spectrum estimation” and “periodogram estimation” are used interchangeably and are intended to relate to essentially the same phenomenon.)
Briefly, in the context of an examined window 102 of a data stream 100, a load shedding arrangement will provide an intelligent sampling scheme (104). At a decision point 106, if there is insufficient CPU time (108), then more points will be removed (110) to yield a spectrum estimation (112) in that context. But given sufficient CPU time (114), then no intermediate step is essentially needed to arrive at spectrum estimation (116). This process will be described in greater detail below.
Essentially, at any given time, there might not be enough processing capacity to provide a periodogram update using all the samples within the data window. The first step toward tackling this problem is the reduction of points using an ‘on-the-fly’ load-shedding scheme. Sub-sampling can lead to data aliasing and deteriorate the quality of the estimated periodogram. Therefore, the sampling should not only be fast but also intelligent, mitigating the impact of the sub-sampling on the squared error of the estimated periodogram. Sampling is based on a linear predictor, which retains a sample only if its value cannot be predicted by its neighbors. An estimator unit is also employed, which changes over time the ‘elasticity’ of the linear predictor, for proper adaptation to the current CPU load.
If there is enough CPU time to process the final number of retained samples, the spectrum is computed. Otherwise, more samples are dropped randomly and the new estimate is computed on the remaining samples.
The computation of the approximate periodogram is based on a formulation of the DFT and the periodogram using unevenly spaced samples, a highly desirable step due to the sampling process. Under a sliding window model, some of the previously used samples are discarded, while new samples are added in the window. A periodicity estimation algorithm proposed herein possesses a very simple update structure, requiring only subtraction of contributions from discarded samples and addition of contributions due to the newly included samples.
Turning now to specific implementations, the Discrete Fourier Transform is used to analyze the frequency content in a discrete and evenly sampled signal. In particular for a discrete time signal x[n] the DFT X [m] is defined for all samples 0≦m, n≦N−1 as:




  
    
      
        
          X
          ⁡
          
            [
            m
            ]
          
        
        =
        
          
            1
            
              N
            
          
          ⁢
          
            
              ∑
              
                n
                =
                0
              
              
                N
                -
                1
              
            
            ⁢
            
              
            
            ⁢
            
              
                x
                ⁡
                
                  [
                  n
                  ]
                
              
              ⁢
              
                e
                
                  
                    -
                    j
                  
                  ⁢
                  
                    
                      2
                      ⁢
                      π
                      ⁢
                      
                        
                      
                      ⁢
                      nm
                    
                    N
                  
                
              
            
          
        
      
    
    
      
        (
        1
        )
      
    
  




The periodogram P of a signal corresponds to the energy of its DFT:

P[m]=∥x[m]∥2  (2)

Consider now, a continuous signal x(t) sampled unevenly at discrete time instants {t0, t1, . . . , tN−1}. An example of this is shown in FIG. 2.
One may write this unevenly sampled signal using the discrete notation as x[kn] where ti=kiT (kiε+) and T corresponds to the sampling interval with all sampling instants as multiples. This is also shown in FIG. 2. (Herethroughout there will be described an adaptive load-shedding algorithm that retains unevenly spaced samples and while an incremental DFT estimation is provided for such discrete signals.)
One can measure the complexity of algorithms in terms of the number of additions (subtractions), multiplications and divisions involved in the computations. Thus, one may label the complexity of a single multiplication as ξMul, of a division as ξDiv and of a sum/subtraction as ξSub.
In connection with a load shedding scheme, one may consider the typical problem of running spectral analysis where we slide a window across the temporal signal and incrementally update the signal's DFT (and the respective periodogram). Preferably, one starts with an evenly sampled signal, with sampling interval T. Consider that the window slides by a fixed amount Width×T. As a result of this sliding we discard n1 points from the beginning of the signal and add n2 points to the end. However, if the available CPU cycles do not allow us to update the DFT using all the points, one can adaptively prune the set of added points using uneven sub-sampling to meet the CPU constraint while minimizing the impact on the accuracy of the updated DFT.
The disclosure now turns to an algorithm (with linear complexity) for the adaptive pruning of the newly added samples. In order to decide whether one can retain a particular sample, one may preferably determine whether it can be linearly predicted from its neighbors. (Higher order predictors are also possible, but would clearly result in higher complexity.) In particular, to make a decision for sample ki one preferably compares the interpolated value xint[ki] with the actual value x[ki], where the interpolated value is computed as:




  
    
      
        
          
            x
            int
          
          ⁡
          
            [
            
              k
              i
            
            ]
          
        
        =
        
          
            
              
                x
                ⁡
                
                  [
                  
                    k
                    
                      i
                      -
                      1
                    
                  
                  ]
                
              
              ⁢
              
                (
                
                  
                    k
                    
                      i
                      +
                      1
                    
                  
                  -
                  
                    k
                    i
                  
                
                )
              
            
            +
            
              
                x
                ⁡
                
                  [
                  
                    k
                    
                      i
                      +
                      1
                    
                  
                  ]
                
              
              ⁢
              
                (
                
                  
                    k
                    i
                  
                  -
                  
                    k
                    
                      i
                      -
                      1
                    
                  
                
                )
              
            
          
          
            
              k
              
                i
                +
                1
              
            
            -
            
              k
              
                i
                -
                1
              
            
          
        
      
    
    
      
        (
        3
        )
      
    
  




where sample ki−1 is the last retained sample before sample ki and sample ki+1 is the immediately following sample. If





  
    
    
      
        
          x
          int
        
        ⁡
        
          [
          
            k
            i
          
          ]
        
      
      -
      
        x
        ⁡
        
          [
          
            k
            i
          
          ]
        
      
    
    
  
  ≤
  
    
      Thresh
      ⨯
      
        
        
          x
          ⁡
          
            [
            
              k
              i
            
            ]
          
        
        
      
    
    100
  




one can discard the sample ki, otherwise it is retained. The parameter Thresh is an adaptive threshold that determines the quality of the approximation. If the threshold is large, more samples are discarded, and similarly if the threshold is small fewer samples are discarded. (It should be noted that the squared approximation error due to this sub-sampling scheme cannot be bounded in general for all signals, however it is selected for its computational simplicity. In particular, for the wide variety of signals considered in experimentation herein, there has not been observed squared error significantly larger than the absolute squared threshold value. Modification of this scheme to guarantee bounds on the approximation error can well be further explored.) An example of this interpolation scheme is shown in FIGS. 3A-3B.

In FIGS. 3A-3B there are shown two steps of the algorithm. In the first step, a decision is made that we can discard sample ki as it can be interpolated by samples k0 and k2. In the next step, a decision is made that we cannot discard sample k2, as it cannot be interpolated using samples k0 and k3, its neighbors. If one starts out with n2 samples need to be pruned, the complexity of this algorithm is:

ξinterp=(2ξMul+4ξSub+ξDiv)(n2−2)  (4)

Further below is a discussion of how to tune the threshold Thresh in order to obtain the desired number of {circumflex over (n)}2 samples, out of the n2 samples added by the sliding window.
In FIGS. 4A-4B, there is illustrated with respect to a stream that measures web usage, a comparison of a presently inventive intelligent sampling method against an equi-sampling technique, which samples data at a specified time interval. Preferably, a presently inventive algorithm is executed for a specific threshold and the data points within a window are reduced from M down to N (unevenly spaced). One preferably estimates the resulting periodogram (see further below) as well the periodogram derived by equi-sampling every N/M points. It is apparent from FIG. 4 that intelligent sampling provides a much higher quality reconstruction of the periodogram, because it can retain important features of the data stream. Additional examples on more datasets are provided in FIGS. 5A-5H.
The load-shedding algorithm assumes the input of a threshold value, which directly affects the resulting number of retained points within the examined window. The desirable number of final points after the thresholding is dictated by the available CPU load. An optimal threshold value would lead to sampling exactly as many points as could be processed by the currently available CPU time. However, there is no way of predicting accurately the correct threshold without having seen the complete data, or without resorting to an expensive processing phase. In FIGS. 6A-6F and FIGS. 7A-7F there are provided various examples of the spectrum approximation for different parameters of the load-shedding threshold value.
A simple estimator of the threshold value with constant complexity can be provided, derived by training on previously seen portions of the data stream. The expectation is that the training will be performed on a data subset that captures a sufficient variation of the stream characteristics. The estimator will accept as input the desired number of final samples that should remain within the examined window, along with a small subset of the current data characteristics, which—in a way—describe its ‘shape’ or ‘state’ (e.g. a subset of the data moments, its fractal dimensionality, etc.). The output of the estimator is a threshold value that will lead (with high expectation) to the desirable number of window samples.
The estimator is not expected to have zero error, but it should lead approximately to the desired compression ratio. In the majority of cases the selected threshold will lead either to higher or lower compression ratio. Intuitively, higher compression (or overestimated threshold) is preferable. This is the case, because then one does not have to resort to the additional phase of dropping randomly some of the retained samples (a sampling that is ‘blind’ and might discard crucial points, such as important local minima or maxima). In experimentation (as discussed further below), it has been verified that this desirable feature is true for the threshold estimator that is presented immediately below.
By way of a training phase, assume that F is a set of features that capture certain desirable characteristics of the examined data window w, and Pε{0, 1, . . . , |w|} describes how many points can be processed at any given time. The threshold estimator will provide a mapping F×P T, where T is a set of threshold values.
It is not difficult to visualize that data whose values change only slightly (or depict small variance of values) do not require a large threshold value. The reverse situation exists for sequences that are ‘busy’, or exhibit large variance of values. With this observation in mind, there is employed the variance within the examined window as a descriptor of the window state. Higher order moments of the data could also be used in conjunction with the variance for improving the accuracy of the predictor. However, for simplicity and for keeping the computational cost as low as possible, there is employed just the variance in our current prototype implementation.
The training phase proceeds as now described. Given the training data, a sliding window is run on them. For each data window there is computed the variance and there is executed the load-shedding algorithm for different threshold values (typically, 20, 40, . . . , 100, 120). After the algorithm execution the remaining number of data points is recorded. This process is repeated for all the extracted data windows. The result of this algorithm will be a set of triplets: [threshold,variance, number of points]. Given this, one can construct the estimator as a mapping f (numPoints,variance) Thresh, where the actual estimator is essentially stored as a 2-dimensional array for constant retrieval time. An example of this mapping is shown in FIG. 8.
It is clear that the training phase is not performed in real-time. However it happens only once (or periodically) and it allows for a very fast prediction step.
It should be pointed out that, even though we assume that the training data will provide ‘sufficient’ clues on the data stream characteristics, the estimator might come upon an input of [variance, numPoints] that has not encountered during the training phase. In this case, one can simply provide the closest match, e.g. the entry that has the closest distance (in the Euclidean sense) to the given variance and number of points. Alternatively, one could provide an extrapolation of the values, in other words, explicitly learn the mapping function. This can be achieved by constructing an RBF network based on the training triplets. Since this approach is significantly more expensive and could present overfitting problems, in the experiments below the former alternative is followed.
Further, over a period of time, the stream characteristics may gradually change, and in the end may differ completely from the training data, hence leading to inconsistent predictions. One can compensate for this by ‘readjusting’ the predictor, by also recording the observed threshold error during the algorithm execution. This will result in a more extended maintenance phase of the estimator, but this cost is bound to pay off in the long run for datasets that exhibit frequent ‘concept drifts’. This extension is not further elucidated upon herein, but it is presently noted as a potential addition for a more complex version of the threshold estimator.
Consider a signal x[ki],0≦i≦N−1, as shown in FIG. 2. Since the DFT is defined only for evenly sampled signals, there is implicitly recreated an evenly sampled signal before computing the DFT. For this, there is again used a linear interpolator (that matches the sub-sampling algorithm), thereby reconstructing a piece-wise linear evenly sampled signal. The DFT of this evenly sampled signal may be computed in terms of the sum of contributions of each of the individual line segments that constitute it. Due to the nature of the linear interpolator the contribution of each line segment to the DFT may be analytically derived in terms of only the endpoints of the segment (i.e. samples in the original unevenly sampled signal) and the distance between them. This means that one does not actually need to interpolate the unevenly sampled signal but can derive a closed form expression for the DFT under the assumption of a linear interpolation scheme. Note that while the time domain signal includes only N (uneven) samples, in order to compute the Discrete Fourier Transform (DFT) of this signal, one will preferably sample the DFT at least M=kN−1−k0 times to avoid time domain aliasing. If we denote by Xn[m] the contributions to the Fourier Transform from each of the N−1 line segments that make up the implicitly recreated evenly sampled signal, then the DFT of the whole signal can be written as:




  
    
      
        
          X
          ⁡
          
            [
            m
            ]
          
        
        =
        
          
            ∑
            
              n
              =
              1
            
            
              N
              -
              1
            
          
          ⁢
          
            
          
          ⁢
          
            X
            
              n
              ⁡
              
                [
                m
                ]
              
            
          
        
      
    
    
      
        (
        5
        )
      
    
  




where for m=1, . . . , M−1




  
    
      
        
          
            X
            n
          
          ⁡
          
            [
            m
            ]
          
        
        =
        
          
            1
            
              
                (
                
                  
                    k
                    n
                  
                  -
                  
                    k
                    
                      n
                      -
                      1
                    
                  
                
                )
              
              ⁢
              
                
                  (
                  
                    
                      2
                      ⁢
                      π
                      ⁢
                      
                        
                      
                      ⁢
                      m
                    
                    M
                  
                  )
                
                2
              
            
          
          ⁡
          
            [
            
              
                
                  (
                  
                    
                      x
                      ⁡
                      
                        [
                        
                          k
                          
                            n
                            -
                            1
                          
                        
                        ]
                      
                    
                    -
                    
                      x
                      ⁡
                      
                        [
                        
                          k
                          n
                        
                        ]
                      
                    
                  
                  )
                
                ⁢
                
                  (
                  
                    
                      e
                      
                        
                          -
                          j
                        
                        ⁢
                        
                          
                            2
                            ⁢
                            π
                            ⁢
                            
                              
                            
                            ⁢
                            
                              mk
                              
                                n
                                -
                                1
                              
                            
                          
                          M
                        
                      
                    
                    -
                    
                      e
                      
                        
                          -
                          j
                        
                        ⁢
                        
                          
                            2
                            ⁢
                            π
                            ⁢
                            
                              
                            
                            ⁢
                            
                              mk
                              n
                            
                          
                          M
                        
                      
                    
                  
                  )
                
              
              +
              
                

              
              ⁢
              
                j
                ⁢
                
                  
                    2
                    ⁢
                    π
                    ⁢
                    
                      
                    
                    ⁢
                    m
                  
                  M
                
                ⁢
                
                  (
                  
                    
                      
                        x
                        ⁡
                        
                          [
                          
                            k
                            n
                          
                          ]
                        
                      
                      ⁢
                      
                        e
                        
                          
                            -
                            j
                          
                          ⁢
                          
                            
                              2
                              ⁢
                              π
                              ⁢
                              
                                
                              
                              ⁢
                              
                                mk
                                n
                              
                            
                            M
                          
                        
                      
                    
                    -
                    
                      
                        x
                        ⁡
                        
                          [
                          
                            k
                            
                              n
                              -
                              1
                            
                          
                          ]
                        
                      
                      ⁢
                      
                        e
                        
                          
                            -
                            j
                          
                          ⁢
                          
                            
                              2
                              ⁢
                              π
                              ⁢
                              
                                
                              
                              ⁢
                              
                                mk
                                
                                  n
                                  -
                                  1
                                
                              
                            
                            M
                          
                        
                      
                    
                  
                  )
                
              
            
            ]
          
        
      
    
    
      
        (
        6
        )
      
    
  




and for m=0





  
    
      
        
          
            X
            n
          
          ⁡
          
            [
            0
            ]
          
        
        =
        
          
            1
            2
          
          ⁢
          
            (
            
              
                x
                ⁡
                
                  [
                  
                    k
                    
                      n
                      -
                      1
                    
                  
                  ]
                
              
              +
              
                x
                ⁡
                
                  [
                  
                    k
                    n
                  
                  ]
                
              
            
            )
          
          ⁢
          
            (
            
              
                k
                n
              
              -
              
                k
                
                  n
                  -
                  1
                
              
            
            )
          
        
      
    
    
      
        (
        7
        )
      
    
  




A significant benefit that equation (5) brings is that the DFT for such unevenly sampled signals can be evaluated incrementally. Hence, if the window is shifted by a fixed width such that the first n1 points are discarded, and n2 points are added at the end, then the DFT of the signal may be updated as follows:




  
    
      
        
          
            X
            new
          
          ⁡
          
            [
            m
            ]
          
        
        =
        
          
            
              X
              old
            
            ⁡
            
              [
              m
              ]
            
          
          -
          
            
              ∑
              
                n
                =
                1
              
              
                n
                1
              
            
            ⁢
            
              
            
            ⁢
            
              
                X
                n
              
              ⁡
              
                [
                m
                ]
              
            
          
          +
          
            
              ∑
              
                n
                =
                N
              
              
                N
                +
                
                  n
                  2
                
                -
                1
              
            
            ⁢
            
              
            
            ⁢
            
              
                X
                n
              
              ⁡
              
                [
                m
                ]
              
            
          
        
      
    
    
      
        (
        8
        )
      
    
  




There will now be considered the complexity of computing this update. As with several papers that analyze the complexity of the FFT, it is assumed that the complex exponentials




  e
  ⁢
  
    
      j2π
      ⁢
      
        
      
      ⁢
      
        mk
        n
      
    
    M
  




(and the intermediate value





  
    2
    ⁢
    π
    ⁢
    
      
    
    ⁢
    
      mk
      n
    
  
  M




) are considered pre-computed for all m and n. Using labels for complexity as defined in the notation, the complexity of computing one single update coefficient Xn[m] for m=1, . . . , M−1—1 may be represented as:

{circumflex over (ξ)}=6ξMul+5ξSub+ξDiv  (9)

and for m=0 as

{circumflex over (ξ)}=2ξMul+2ξSub  (10)

Finally, the complexity of updating all the M DFT coefficients in this scenario is:




  
    
      
        
          
            ξ
            update
          
          ⁡
          
            (
            
              M
              ,
              
                n
                1
              
              ,
              
                n
                2
              
            
            )
          
        
        =
        
          
            
              (
              
                
                  n
                  1
                
                +
                
                  n
                  2
                
              
              )
            
            ⁡
            
              [
              
                
                  
                    (
                    
                      M
                      -
                      1
                    
                    )
                  
                  ⁢
                  
                    (
                    
                      
                        6
                        ⁢
                        
                          ξ
                          Mul
                        
                      
                      +
                      
                        5
                        ⁢
                        
                          ξ
                          Sub
                        
                      
                      +
                      
                        ξ
                        Div
                      
                    
                    )
                  
                
                +
                
                  (
                  
                    
                      2
                      ⁢
                      
                        ξ
                        Mul
                      
                    
                    +
                    
                      2
                      ⁢
                      
                        ξ
                        Sub
                      
                    
                  
                  )
                
                +
                
                  M
                  ⁢
                  
                    
                  
                  ⁢
                  
                    ξ
                    Sub
                  
                
              
              ]
            
          
          +
          
            2
            ⁢
            M
            ⁢
            
              
            
            ⁢
            
              ξ
              Sub
            
          
        
      
    
    
      
        (
        11
        )
      
    
  




Using a presently inventive sub-sampling algorithm, one can reduce the number of samples that need to be used to update the DFT. Consider that as a result of the pruning, one can reduce n2 samples into a set of {circumflex over (n)} samples ({circumflex over (n)}2≦n2). While the reduction in the number of samples directly translates to a reduction in the complexity of the update, one should also factor in the additional cost of the sub-sampling algorithm. Comparing equations (11) and (4) it is apparent that the overall complexity of the update (including the sub-sampling) is reduced when:

ξupdate(M,n1,n2)≧ξupdate(N,n1,{circumflex over (n)}2)+ξinterp  (12)

To determine when this happens, consider a simple case when {circumflex over (n)}2=n2−1, i.e. the sub-sampling leads to a reduction of one sample. The increase in complexity for the sub-sampling is (2ξMul+4ξSub+ξDiv)(n2−2) while the corresponding decrease in the update complexity is (M−1)(6ξMul+5ξSub+ξDiv)+(2ξMul+2ξSub)+MξSub (from equation (11). Clearly, since {circumflex over (n)}2<n2≦M, one can easily realize that the reduction in complexity far outweighs the increase due to the sub-sampling algorithm. In general, equation (12) is always true when the sub-sampling algorithm reduces the number of samples (i.e., when {circumflex over (n)}2<n2).
If, at a certain time, the CPU is busy, thereby imposing a computation constraint of ξlimit, one should preferably perform a DFT update within this constraint. If ξupdate(M,n1,n2)>ξlimit one cannot use all the samples n2 for the update, and hence one needs to determine the optimal number of samples to retain {circumflex over (n)}2, such that ξupdate(M,n1,{circumflex over (n)}2)+ξinterp≦ξlimit. Specifically, one may compute this as:




  
    
      
        
          
            n
            ^
          
          2
        
        ≤
        
          
            
              
                ξ
                limit
              
              -
              
                ξ
                interp
              
              -
              
                2
                ⁢
                
                  Mξ
                  Sub
                
              
            
            
              
                
                  (
                  
                    M
                    -
                    1
                  
                  )
                
                ⁢
                
                  (
                  
                    
                      6
                      ⁢
                      
                        ξ
                        Mul
                      
                    
                    +
                    
                      5
                      ⁢
                      
                        ξ
                        Sub
                      
                    
                    +
                    
                      ξ
                      Div
                    
                  
                  )
                
              
              +
              
                (
                
                  
                    2
                    ⁢
                    
                      ξ
                      Mul
                    
                  
                  +
                  
                    2
                    ⁢
                    
                      ξ
                      Sub
                    
                  
                
                )
              
              +
              
                M
                ⁢
                
                  
                
                ⁢
                
                  ξ
                  Sub
                
              
            
          
          -
          
            n
            1
          
        
      
    
    
      
        (
        13
        )
      
    
  




Finally, one can achieve this by tuning the sub-sampling threshold Thresh based on the threshold estimator algorithm described hereinabove.
The disclosure now turns to experimentation with the algorithms and concepts discussed hereinabove.
The usefulness of the presently inventive resource-adaptive periodicity estimation depends on two factors:

    
    
        (1) The accuracy of the approach, which is indicated by the quality of the DFT approximation and its respective periodogram. If the periodogram after the load-shedding closely resembles the original one, then the provided estimate is meaningful.
        (2) The adaptiveness of the proposed methodology, which is highly dependent on the quality of the threshold estimator. An accurate estimator will lead to sampling rates that closely adapt to the current CPU loads.

The above two factors are examined separately herebelow.

    
    


The quality of the approximated Fourier coefficients is measured on a variety of periodic datasets obtained from the time-series archive at UC Riverside. These datasets only have a length of 1024, therefore it is difficult to provide a meaningful evaluation on the streaming version of the algorithm. However, by providing the whole sequence as input to the periodicity estimation unit one can evaluate the effectiveness of the load-shedding scheme in conjunction with the closed-form DFT computation on the unevenly spaced samples. One can compute the accuracy by comparing the estimated periodogram against the actual one (had we not discarded any point from the examined data window). Thus, the above experiment was run on different threshold values Thresh=20 . . . 120. For example, a value of Thresh=20 signifies that the predicted value (using the linear predictor) does not differ more than 20% from the actual sequence value.
Note that the original periodogram is evaluated on a window of M points (M=1024), while the one based on uneven sampling uses only the N remaining samples (N≦M). In order to provide a meaningful comparison between them the latter periodogram is evaluated on all M/2 frequencies—see equation (6)—, even though this is not necessary with an actual deployment of the algorithm.
The accuracy of the presently inventive methodology is compared against a naive approach that uses equi-sampling every N/M points (i.e., leading again to N remaining points within the examined window). This approach is bound to introduce aliasing and distort more the original periodogram, because (unlike the intelligent loadshedding) it does not adapt according to the signal characteristics.
FIGS. 9A-9H and FIGS. 10A-10F indicate the periodogram error introduced by the intelligent and the equi-sampling techniques. On top of each bar there is also portrayed the compression achieved using the specific threshold Thresh, computed as 100*(1−N/1024).
The results suggest that the load-shedding scheme employed by the presently inventive technique can lead to spectrum estimates of much higher quality than competing methods. In two cases (FIG. 9D, Reality Check) the equi-sampling performs better than the linear interpolator, but this occurs only for minute compression ratios (i.e., when the threshold discards less than 10 samples per 1024 points). In general the observed reduction in the estimation error compared to equi-sampling, can range from 10% to more than 90% on the 14 datasets examined herein.
To test the accuracy of the threshold estimator, longer datasets are needed, which could be used for simulating a sliding window model execution and additionally provide a training subset. Thus, for experimentation purposes, there were utilized real datasets provided by the automotive industry. As such, these are diagnostic measurements that monitor the evolution of variables of interest during the operation of a vehicle. Examples of such measurements could be the engine pressure, the torque, vibration patterns, instantaneous fuel economy, engine load at current speed, etc.
Periodic analysis is an indispensable tool in automotive industry, because predictive maintenance can be possible by monitoring the changes in the spectrum of the various rotating parts. Therefore, a change in the periodic structure of the various engine measurements can be a good indicator of machine wear and/or of an incipient failure.
The measurements used have a length of 50000 points and represent the monitoring of a variable over an extended period of time. On this data, there was employed a sliding window of 1024 points. A synthetic CPU load is generated, which is provided as input to the periodicity estimation unit. Based on the synthetic CPU trace, at any given point in time the periodicity unit is given adequate time for processing a set of points with cardinality within the range of 50 to 1024 (1024 being the length of the window). In FIGS. 11A-11C there are depicted two instances of the approximated spectrum under limited CPU resources. On the first instance the indicated available CPU of 12.41% means that only 12.41% of the total window points should remain after the load-shedding, given the available processing time.
Executing the presently inventive algorithm on the complete data stream, the accuracy of the threshold estimator is monitored. The estimator is fed with the current CPU load and provides a threshold estimate Threshest that will lead with high probability to {circumflex over (P)} remaining points (so that they could be sufficiently processed given the available CPU load). Suppose that the actual remaining points after the application of the threshold Threshest, are P. An indicator of the estimator accuracy is provided by contrasting the estimated number of points {circumflex over (P)} against the actual remaining ones P(error=|{circumflex over (P)}−P|).
The experimental results are very encouraging and indicate an average error on the estimated number of points in the range of 5% of the data window. For this experiment, if the predicted number of points for a certain threshold is 250 points, the actual value of remaining points could be (for example) 200 points. This is the case of an overestimated threshold which compressed more the flowing data stream. As mentioned before, this case is more desirable (than an underestimated threshold), because no additional points need to be subsequently dropped from the current data window (which is not bound to introduce additional aliasing problems).
A histogram of the estimator approximation error is in FIG. 12A. It is observed that for the majority of data windows the estimation error is small, while fewer instances of the algorithm execution report a large error in the threshold estimation. In FIG. 12B there is an indication of how many cases of overestimated thresholds there are and how many are underestimated. The overestimated ones (more desirable) are higher than the underestimated, which again indicates many of the attractive properties of the proposed threshold predictor.
In summary, there has been presented herein the first resource-adaptive method for periodicity estimation for streaming data. By way of brief, albeit non-restrictive recapitulation, some key aspects of a proposed method in accordance with at least one embodiment of the present invention are:
(1) An intelligent load-shedding scheme that can adapt to the CPU load using a lightweight predictor.
(2) A DFT estimation that utilizes unevenly spaced samples, provided by the previous phase.
The quality of the approximated DFT has been shown and it has also been demonstrated that the scheme can adapt closely to the available CPU resources. The intelligent load-shedding scheme has been compared against equi-sampling and improvements in the periodogram estimation ranging from 10% to 90% are shown.
Further exploration could involve an examination of whether it is possible to reduce the computational cost even further.
By way of further recapitulation, some of the important contributions set forth herein are as follows:

    
    
        An abstraction of the resource adaptation problem for periodicity estimation is provided.
        An intelligent load-shedding scheme along with a parameter estimator unit that tunes the adaptation to the current CPU load is proposed.
        A closed-form Fourier approximation is presented using unevenly spaced samples and it is shown how to update it incrementally.
        The performance of a presently inventive approach is analyzed under CPU constraints, and the complexity is measured abstractly, in terms of the number of multiplications, additions and divisions involved (making the analysis independent of the underlying processor architecture). Even though the model is very spartan in its memory utilization, no memory constraints are explicitly imposed. However, the inclusion of potential memory constraints (and other resource constraints) would be a straightforward addition to the model.
    
    


Other recent work on periodicity estimation on data streams has appeared, where the authors study sampling techniques for period estimation using sublinear space. Other work proposes sampling methods for retaining (with a given approximation error) the most significant Fourier coefficients. Still other work adapts the use of wavelet coefficients for modeling a data stream, providing also a periodicity estimator using logarithmic space complexity. However, none of these approaches addresses the issue of resource adaptation, which is one of the main contributions provided herein. It should be further noted that a presently proposed method for periodogram reconstruction based on irregularly spaced samples is significantly more lightweight than the widely used Lomb periodogram (which incurs a very high computational burden).
It is to be understood that the present invention, in accordance with at least one presently preferred embodiment, includes an arrangement for providing a spectrum estimation based on resource availability, which may be implemented on at least one general-purpose computer running suitable software programs. It may also be implemented on at least one integrated Circuit or part of at least one Integrated Circuit. Thus, it is to be understood that the invention may be implemented in hardware, software, or a combination of both.
Although illustrative embodiments of the present invention have been described herein with reference to the accompanying drawings, it is to be understood that the invention is not limited to those precise embodiments, and that various other changes and modifications may be affected therein by one skilled in the art without departing from the scope or spirit of the invention.

## Claims

1. A method comprising:
accessing one or more data streams;
utilizing one or more processors to perform:
determining a collection of data points for providing spectrum estimation for data in the one or more data streams;
determining whether there is sufficient resource availability for providing spectrum estimation with the collection of data points, wherein resource availability comprises one or more of processing capacity and memory;
removing data points from the collection of data points according to a load shedding scheme responsive to determining that there is insufficient resource availability for providing spectrum estimation with the collection of data points; and
providing a periodogram based on any remaining data points of the collection of data points;
wherein one or more data points of the collection of data points are removed if the one or more data points can be predicted from neighboring data points.

2. The method according to claim 1, determining a collection of data points comprises:
sampling data from the one or more data streams to provide said collection of data points.

3. The method according to claim 2, wherein:
said sampling comprises:
sampling data from the one or more data streams within a data window; and
wherein said load shedding scheme comprises determining whether to retain or discard sampled data points based on resource availability by applying a data predictor, the data predictor prompting retention of a data point of said collection of data points only if the value of the data point cannot be predicted by previously encountered portions of the one or more data streams.

4. An apparatus comprising:
one or more processors;
a system memory operatively connected to the one or more processors;
wherein, responsive to execution of computer readable program code accessible to the one or more processors, the one or more processors are configured to:
access one or more data streams;
determine a collection of data points to provide spectrum estimation using a collection of data points for data in the one or more data streams;
determine whether there is sufficient resource availability for providing spectrum estimation with the collection of data points;
remove data points from the collection of data points according to a load shedding scheme responsive to determining that there is insufficient resource availability for providing spectrum estimation with the collection of data points; and
provide a periodogram based on any remaining data points of the collection of data points;
wherein one or more data points of the collection of data points are removed if the one or more data points can be predicted from neighboring data points.

5. The apparatus according to claim 4, determining a collection of data points comprises:
sampling data from the one or more data streams.

6. The apparatus according to claim 5, wherein:
said sampling comprises:
sampling data from the one or more data streams within a data window; and

wherein said load shedding scheme comprises determining whether to retain or discard sampled data points based on resource availability by applying a data predictor, the data predictor prompting retention of a data point only if the value of the data point cannot be predicted by previously encountered portions of the one or more data streams.

7. A tangible, non-transitory program storage device, readable by machine, embodying a program of instructions that when executed by a processor of the machine enable the machine to:
accessing one or more data streams;
determining a collection of data points for providing spectrum estimation for data in the one or more data streams;
determining whether there is sufficient resource availability for providing spectrum estimation with the collection of data points;
removing data points from the collection of data points according to a load shedding scheme responsive to determining that there is insufficient resource availability for providing spectrum estimation with the collection of data points; and
providing a periodogram based on any remaining data points of the collection of data points;
wherein one or more data points of the collection of data points are removed if the one or more data points can be predicted from neighboring data points.