Fourier transform for timestamped network data

Abstract

Spectrum estimates of unevenly spaced timestamped data collected over a network. Unevenly spaced data is not suitable for processing by traditional methods such as Fast Fourier Transforms. Spectrum estimates of timestamped data collected over a network are calculated using a Continuous Fourier Transform. If the time between samples is long, spectrum estimates may be computed incrementally.

Description

FIELD OF THE INVENTION

The present invention pertains to dealing with measurements taken at unevenly-spaced intervals in a computer network, and deriving a frequency spectrum from these measurements.

ART BACKGROUND

Many companies, such as Agilent Technologies, Hewlett-Packard, and others, support sensor nodes based on IEEE 1451 standards. Such standards allow the deployment of sensors on networks ranging from dedicated Ethernet local area networks to corporate internets to the World Wide Internet.

Current TCP/IP networks currently offer no quality-of-service guarantees. This causes a problem for many measurement applications, such as difficulty taking measurements at accurate and pre-determined time intervals.

U.S. Pat. No. 5,566,180 to Eidson et. al. teaches a method and apparatus for providing precise control of timing in distributed nodes in a network. U.S. Pat. No. 6,278,710 to Eidson provides further enhancement to time synchronization protocols for a distributed system.

These protocols provide the ability to accurately timestamp measurements taken by nodes on a network. Nevertheless, the times that measurements are taken may not be evenly spaced.

Many applications require the processing of data comprising sequences of such timestamped measurements. But many commonly used applications, such as signal processing algorithms, operate on data collected at evenly spaced times. One such common signal processing task is providing spectral data through the calculation of a Fourier Transform. What is needed is a method of calculating Fourier Transforms on unevenly spaced timestamped network data.

SUMMARY OF THE INVENTION

Spectrum estimates of unevenly spaced timestamped network data are calculated using estimates of a Continuous Fourier Transform. If the time between samples is long, the estimates may be computed incrementally.

BRIEF DESCRIPTION OF THE DRAWINGS

The present invention is described with respect to particular exemplary embodiments thereof and reference is made to the drawings in which:

FIG. 1 shows a block diagram of network data collection,

FIG. 2 is a graph of data collected over a network,

FIG. 3 is a graph of time between data samples collected over a network, and

FIG. 4 is a spectrum calculated from the collected data.

DETAILED DESCRIPTION

FIG. 1 shows a block diagram of network data collection. Signal source 100 provides a signal to smart transducer interface module (STIM) 110 . STIM 110 communicates 120 with Ethernet Controller 130 . Ethernet Controller 130 communicates 140 to the Internet 150 . Workstation 170 also communicates 160 with the Internet 150 . In a data collection application such as that contemplated by the present invention, Data Acquisition software module 180 collects data from source 100 through STIM 110 through the communications paths shown. Fourier Transform software module 190 uses the data collected by Data Acquisition module 180 and calculates a frequency spectrum of the data. In the preferred embodiment, Data Acquisition module 180 is written in Java to maintain portability. Fourier Transform module 190 is written in Matlab. These two modules may be implemented in many different computer languages, and may be combined.

The Internet is a wild beast. When protocols such as TCP/IP are used, packet delivery time is not guaranteed. The routing and therefore transit times of individual packets between end nodes such as Ethernet Controller 130 and workstation 170 may vary on a packet-to-packet basis. To this may be added delays interposed by proxy servers, routers, switches, other Internet traffic, and the like.

In the example of FIG. 1 , Signal Source 100 is an Agilent Technologies 33120A Function Generator producing an amplitude modulated signal with a center frequency of 0.01 Hz modulated at 0.01 Hz. This produces a signal with strong peaks at 0.01 Hz and 0.02 Hz.

Workstation 170 has a Java-written data acquisition program 180 acquiring 1000 data samples of Signal Source 100 through STIM 110 and Ethernet Controller 130 . Data acquisition program 180 makes a separate HTTP connection for each data sample. One connection with multiple data requests could also be performed. Each data sample includes an accurate time stamp as well as the value of Signal Source 100 .

FIG. 2 is a graph of the collected data. Note that the distribution of points in the horizontal, or time axis, is not uniform. The inter-measurement delay is better seen in FIG. 3 , showing an almost four-to-one variation in sample acquisition times.

The variability in sample acquisition time makes calculating spectral data using traditional approaches difficult.

The common approach to deriving spectral data from a series of measurements is to calculate a Discrete Fourier Transform or a Fast Fourier Transform. These are derived from the Continuous Fourier Transform.

The Fourier Transform F(ω) of a real or complex function x(t) is $$ F ⁢   ⁢ ( ω ) = 1 2 ⁢   ⁢ π ⁢   ⁢ ∫ - ∞ ∞ ⁢ ⅇ 2 ⁢   ⁢ π ⁢   ⁢ j ⁢   ⁢ ω ⁢   ⁢ t ⁢   ⁢ x ⁢   ⁢ ( t ) ⁢   ⁢ ⅆ t ( 1 )

For any real frequency ω and where j={square root over (−1)}. For the present invention, we will deal with t measured in seconds and ω to be frequency measured in Hertz. In most signal processing applications, the signal is sampled at a constant time interval Δt yielding N samples x(0),x(Δt), . . . ,x((N−1)Δt). The Discrete Fourier Transform (DFT) arises as the discrete (summation) form of Equation (1) above: $$ F ⁢   ⁢ ( ω i ) = 1 N ⁢   ⁢ ∑ k = 0 N - 1 ⁢   ⁢ ⅇ 2 ⁢   ⁢ π ⁢   ⁢ j ⁢   ⁢ ω i ⁢   ⁢ k ⁢   ⁢ Δ ⁢   ⁢ t ⁢   ⁢ x ⁢   ⁢ ( k ⁢   ⁢ Δ ⁢   ⁢ t ) ( 2 )

Computing the sums in Equation (2) in the straightforward way requires O(N 2 ) computations. The Fast Fourier Transform (FFT) is a means for calculating these sums in O(N log N) computations. Unfortunately, the correctness of the derivations of both the DFT and FFT depend on the sampling rate Δt being constant. Thus the Discrete and Fast Fourier Transforms are not suited to handling unevenly sampled timestamped data.

In the present invention, let

t 0 ,t 1 , . . . t N−1   (3)

be the times the signal was sampled, and let

x 0 ,x 1 , . . . x N−1   (4)

be the measured values of the signal x(t) at those times, i.e. x i =x(t i ).

It is also useful to assume that the time stamps are precise, having many significant digits provided, and accurate. This condition is met if the accuracy of the clock being used to timestamp ticks in fine time quanta is small in comparison to the time constants of the physical system being measured. When this condition is met, time may be considered to be continuous.

In order to evaluate Equation (1), it is assumed that the signal x(t) is identically zero before and after the time during which it is measured. Other reasonable assumptions on the signal x(t) in between the times it is measured include:

Zero Order Hold (ZOH) The signal maintains a constant value up until the time of the next sample. That is, if t i ≦t<t i+1 , then x(t)=x(t i ). The Zero Order Hold assumption is reasonable if the signal is oversampled.

First Order Hold (FOH) The signal is continuous and piecewise linear, with first derivative discontinuities at the sampling times. That is, if t i ≦t<t i+1 , then

x ( t )=λ x ( t i )+(1−λ) x ( t i+1 ) where λ=( t−t i )/( t i+1 −t i ).

Locally Polynomial The signal is locally given by a quadratic, cubic, or other polynomial function that interpolates or fits several adjacent samples.

Bandlimited The signal's spectrum must have no power outside of a particular frequency band.

One approach to dealing with unevenly sampled timestamped data is to resample the data evenly, interpolating between the data in order to yield evenly spaced samples, and then perform an FFT on the resampled data. In resampling, any of the assumptions (zero order hold, first order hold, polynomial, etc.) can be used to determine resampled values at the desired evenly spaced time points.

The Fast Fourier Transform (FFT) works most efficiently if its input has a length which is an integer power of two. Choose an integer k>log 2 N (for example, k=┌log 2 N┐) and let N′=2 k be the number of points to be in the re-sampled data. The times of the re-sampled data will then be t′ 0 , . . . ,t′ N−1 , where $$ t i ′ = i N - 1 ⁢   ⁢ ( t N - 1 - t 0 ) + t 0 .

Resampling is performed according to the assumptions made on the nature of the data. In the simplest case, assuming zero-order-hold, the resampled datum x′ i at time t′ i is x′ i =x j where j is the largest integer such that t j ≦t′ i . If first-order-hold is assumed, then the resampled datum x′ i at time t′ i is a linear interpolation of the two closest points x′ i =λx j +(1−λ)x j+1 , where again j is the largest integer such that $$ t j ≤ t i ′ ⁢   ⁢ and ⁢   ⁢ λ = t j + 1 - t i ′ t j + 1 - t j .

Higher order scaling functions may also be applied, such as quadratic, cubic, or other polynomial methods.

According to the present invention, assuming zero-order-hold on the sampled signal, the Fourier Transform F(ω) of Equation (1) becomes $$ F ⁢   ⁢ ( w ) = 1 2 ⁢   ⁢ π ⁢   ⁢ ∑ i = 0 N - 1 ⁢   ⁢ ∫ t i t i + 1 ⁢ ⅇ 2 ⁢   ⁢ π ⁢   ⁢ j ⁢   ⁢ ω ⁢   ⁢ t ⁢   ⁢ x ⁢   ⁢ ( t ) ⁢   ⁢ ⅆ t ( 5 )   ⁢ = 1 2 ⁢   ⁢ π ⁢   ⁢ ∑ i = 0 N - 1 ⁢   ⁢ ∫ t i t i + 1 ⁢ ⅇ 2 ⁢   ⁢ π ⁢   ⁢ j ⁢   ⁢ ω ⁢   ⁢ t ⁢   ⁢ x ⁢   ⁢ ( t i ) ⁢   ⁢ ⅆ t ( 6 )   ⁢ = 1 ( 2 ⁢   ⁢ π ) 3 / 2 ⁢   ⁢ ω ⁢   ⁢ j ⁢   ⁢ ∑ i = 0 N - 1 ⁢   ⁢ ( ⅇ 2 ⁢   ⁢ π ⁢   ⁢ j ⁢   ⁢ ω ⁢   ⁢ t i + 1 - ⅇ 2 ⁢   ⁢ π ⁢   ⁢ j ⁢   ⁢ ω ⁢   ⁢ t i ) ⁢   ⁢ x i ( 7 )

Normally, F(ω) is computed at about N/2 frequencies, requiring on the order of O(N 2 ) computations to calculate Equation (7). For moderate values of N this may be acceptable.

If the time between samples is long, the Fourier Transform may be computed incrementally. Initially set F(ω)←0 for each frequency ω for which a spectral estimate is desired. Then after each sample x i measured at time t i is received, set $$ F ⁢   ⁢ ( ω ) ← F ⁢   ⁢ ( ω ) + 1 ( 2 ⁢   ⁢ π ) 3 / 2 ⁢   ⁢ ω ⁢   ⁢ j ⁢   ⁢ ∑ i = 0 N - 1 ⁢   ⁢ ( ⅇ 2 ⁢   ⁢ π ⁢   ⁢ j ⁢   ⁢ ω ⁢   ⁢ t i + 1 - ⅇ 2 ⁢   ⁢ π ⁢   ⁢ j ⁢   ⁢ ω ⁢   ⁢ t i ) ⁢   ⁢ x i ( 8 )

This incremental process can be used to provide an estimate of the spectrum which is updated after the arrival of each new measurement.

Equation (7) may be implemented for example as the Matlab source program shown as Computer Listing 1 at the end of this application. In this embodiment, time is rescaled to the interval 0≦t≦2π. This resealing is useful as the values for the timestamps may be numerically quite large, for example the number of seconds from an epoch date such as 1 Jan. 1970. Computing sines, cosines or exponentials of large arguments may produce inaccurate results. The Continuous Fourier Transform is calculated on scaled frequencies, and the resulting frequencies are rescaled to correspond to the original time scale.

Applying the program of Computer Listing 1 to the data shown in FIG. 2 produces the spectrum shown in FIG. 4 . This spectrum shows strong peaks at 0.01 Hz and 0.02 Hz.

The foregoing detailed description of the present invention is provided for the purpose of illustration and is not intended to be exhaustive or to limit the invention to the precise embodiments disclosed. Accordingly the scope of the present invention is defined by the appended claims.

Computer Listing 1
function [freqs, y] = FT_timestamped(t, x, algorithm)
% [freqs, y] = FT_timestamped(t, x)
%
% t      Real vector: times that the data were measured.
% x      Vector: the data.
% freqs  Real vector: frequencies in Hertz.
% y      Complex vector: y(i) is the complex amplitude of the data at
         frequency freqs(i).
% Notes: length(t) must equal length(x).
%
if( nargin < 3)
         algorithm = 2;
end
len_t = length(t);
i = 1:len_t−1; % Indexes the whole of t & x EXCEPT for the last
element
% Check parameter validity
mbrealvector(t);
mbvector(x);
if( len_t ˜= length(x) )
         error(‘t and x do not have the same length’);
end
t=t(:);
x=x(:);
% Max freq. per Nyquist criterion
max_Hz = len_t / (2*(t(end)−t(1)))
nfreqs = ceil(len_t / 2)
freqs = [1:nfreqs−1];
if algorithm==1
         % Algorithm 1: Continuous Fourier Transform
         % Shift and scale time to lie in the interval [0, 2*pi].
         % This scaling is important: algorithms for computing
         % sines, cosines, or comlex exponentials generally become
         % inaccurate for large arguments.
         t = t − t (1);
         t_end = t(end)
         scale = 2*pi / t(end);
         t = scale * t;
         y = zeros(size(freqs)); % Preallocate space for y
         for omega = freqs
%        y(omega) = sum( (exp(j*omega*t(i))−exp(j*omega*t(i+1))).*x(i) )
. . .
%        * (j/(omega*sqrt(2*pi)));
         y(omega) = sum( (exp(j*omega*t(i))−exp(j*omega*t(i+1))).*x(i) )
./ omega;
         end
         y = y * (j/(2*pi));
         % Scale the frequencies to correct time scale
         freqs = (max_Hz / nfreqs) * freqs;
         dc = sum( x(1:end−1) .* (t(2:end)−t(1:end−1)) ) ./ (t(end)−t(1));
         freqs = [−freqs(end:−1:1) 0 freqs];
         y = [conj(y(end:−1:1)) dc y];
elseif (algorithm==2) | (algorithm==3)
         % Algorithm 2 & 3: Resampling followed by FFT
         len_newt = 2{circumflex over ( )}ceil(log2(len_t)) % Smallest power of 2 that is >=
len_t
         newt = 0:len_newt−1;
         newt = newt*((t(end)−t(1))/(len_newt−1)) + t(1);
         new_max_Hz = len_newt / (2*(t(end)−t(1)))
         newt1_minus_t1 = newt(1)-t(1)
         newtend_minus_tend = newt(end)−t(end)
         newx = zeros(size(newt));
         size_newx = size(newx)
         nextSampleIndex = 1;
         % ZOH resampling in O(N) time
         %t = [t ; inf];
         for i=1: length (newt)−1
         while( ˜( newt(i) < t(nextSampleIndex+1) ) )
         nextSampleIndex = nextSampleIndex + 1;
         end
         switch algorithm
         case 2,
         % Zero Order Hold
         newx(i) = x(nextSampleIndex);
         case 3,
         % First Order Hold
         lambda = (t(nextSampleIndex+1)−
newt(i))/(t(nextSampleIndex+1)−t(nextSampleIndex));
         newx(i) = lambda*x(nextSampleIndex) + (1−
lambda)*x(nextSampleIndex+1);
         end
         end
         newx(end) = x(end);
         %diffx = max(abs(newx(:)−x(:)))
         y = fftshift(fft(newx)) ./ length(newx);
         %freqs = ([−len_newt/2:len_newt/2−1]) * (new_max_Hz /
(len_newt/2));
         freqs = ([−len_newt/2:len_newt/2−1]) ./ (t(end)−t(1));
         size_y = size(y)
         size_f = size(freqs)
end

Claims

1. A method of calculating spectrum estimates comprising:collecting timestamped data over a computer network, and calculating spectrum estimates using a Continuous Fourier Transform where the spectrum estimates F(ω) of data samples having values xi timestamped at time ti are calculated as F⁢ ⁢(ω)=1(2⁢ ⁢π)3/2⁢ ⁢ω⁢ ⁢j⁢ ⁢∑i=0N-1⁢ ⁢(ⅇ2⁢ ⁢π⁢ ⁢j⁢ ⁢ωi+1-ⅇ2⁢ ⁢π⁢ ⁢j⁢ ⁢ωi)⁢ ⁢xi.

2. A method of calculating spectrum estimates comprising:collecting timestamped data over a computer network, and calculating spectrum estimates using a Continuous Fourier Transform where the spectrum estimates F(ω) of data samples having values xi timestamped at time ti are calculated incrementally as F⁢ ⁢(ω)←F⁢ ⁢(ω)+1(2⁢ ⁢π)3/2⁢ ⁢ω⁢ ⁢j⁢ ⁢∑i=0N-1⁢ ⁢(ⅇ2⁢ ⁢π⁢ ⁢j⁢ ⁢ω⁢ ⁢ti+1-ⅇ2⁢ ⁢π⁢ ⁢j⁢ ⁢ω⁢ ⁢ti)⁢ ⁢xi.

3. A computer system for computing spectrum estimates comprising:means for collecting timestamped data from a device connected to a network, and means for computing a Continuous Fourier Transform on the collected timestamped network data producing spectrum estimates where the spectrum estimates F(ω) of data samples having values xi timestamped at time ti are calculated as F⁢ ⁢(ω)=1(2⁢ ⁢π)3/2⁢ ⁢ω⁢ ⁢j⁢ ⁢∑i=0N-1⁢ ⁢(ⅇ2⁢ ⁢π⁢ ⁢j⁢ ⁢ωi+1-ⅇ2⁢ ⁢π⁢ ⁢j⁢ ⁢ωi)⁢ ⁢xi.

4. A computer system for computing spectrum estimates comprising:means for collecting timestamped data from a device connected to a network, and means for computing a Continuous Fourier Transform on the collected timestamped network data producing spectrum estimates where the spectrum estimates F(ω) of data samples having values xi timestamped at time ti are calculated incrementally as F⁢ ⁢(ω)←F⁢ ⁢(ω)+1(2⁢ ⁢π)3/2⁢ ⁢ω⁢ ⁢j⁢ ⁢∑i=0N-1⁢ ⁢(ⅇ2⁢ ⁢π⁢ ⁢j⁢ ⁢ω⁢ ⁢ti+1-ⅇ2⁢ ⁢π⁢ ⁢j⁢ ⁢ω⁢ ⁢ti)⁢ ⁢xi.

5. A computer readable medium carrying one or more sequences of instructions from a user of a computer system for computing spectrum estimates of timestamped network data, wherein execution of the one or more sequences of instructions by one or more processors causes the one or more processors to perform the steps of:collecting timestamped data from a network device, and calculating a Continuous Fourier Transform on the collected timestamped data producing spectrum estimates where the spectrum estimates F(ω) of data samples having values xi timestamped at time ti are calculated as F⁢ ⁢(ω)=1(2⁢ ⁢π)3/2⁢ ⁢ω⁢ ⁢j⁢ ⁢∑i=0N-1⁢ ⁢(ⅇ2⁢ ⁢π⁢ ⁢j⁢ ⁢ωi+1-ⅇ2⁢ ⁢π⁢ ⁢j⁢ ⁢ωi)⁢ ⁢xi.

6. The computer readable medium of claim 5 where the spectrum estimates F(ω) of data samples having values xi timestamped at time ti are calculated incrementally as F⁢ ⁢(ω)←F⁢ ⁢(ω)+1(2⁢ ⁢π)3/2⁢ ⁢ω⁢ ⁢j⁢ ⁢∑i=0N-1⁢ ⁢(ⅇ2⁢ ⁢π⁢ ⁢j⁢ ⁢ω⁢ ⁢ti+1-ⅇ2⁢ ⁢π⁢ ⁢j⁢ ⁢ω⁢ ⁢ti)⁢ ⁢xi.