Measurement scale for non-uniform data sampling in N dimensions

[0002] This application claims the benefit of provisional application No. 60/115,752 filed Jan. 13, 1999.

[0003] This invention relates to methods for sampling objects. In particular it provides a system implementing a rigorous mathematical framework and provably correct practical algorithmic implementations for digitally recording a signal having any number of degrees of freedom.

[0004] A physical or virtual object is recorded for application specific analysis by sampling or measuring it at various points. An object might be a one dimensional electrical signal, price history or other time series, a two dimensional image, a three dimensional vector field, or an N dimensional surface in, say, configuration space, momentum space, or phase space, all useful in physics. One important example is the current use of uniformly spaced samples of signals or images to analyze these objects for subsequent processing like filtering or convolution. This analysis involves calculating approximations to the Fourier transform, wavelet transform or other related transforms. The aggregate of sample points generated by the measuring process creates a mosaic from the original physical object. The degree of congruence of the mosaic to the original depends on the number of sample points used as well as on the judicious location of the sample points.

[0005] A classical example of such is recording of an audio signal for replay by sampling the signal amplitude at uniform time intervals less than half the period of the highest frequency present in the signal. The so-called Nyquist theorem assures that such a recording enables reconstruction of the original signal. If the uniform sampling occurred at a lesser rate, in many cases an attempt to replay the sampled signal would result in a distortion termed “artifacting” in which harmonics of the signal would not correctly reproduce.

[0006] Limited resources place constraints upon the generation of sample points. Therefore sampling is often a fairly expensive operation in terms of system resources consumed. Depending on the application, the number of sample points generated might be constrained by the resolution of the measuring equipment used (cost constraint), by the amount of storage available to save the samples (space constraint) as well as by the observation time window of opportunity to obtain the samples (time constraint). Thus, a set of samples can be regarded as a somewhat scarce resource. Therefore, each sample point should have impact and provide a maximal contribution of information content to the whole mosaic.

[0007] Since there is a cost associated with the generation of each sample point, some optimization of sampling is generally required. Particularly in the case of objects having periodically recurring patterns there is much redundancy in samples taken in uniform intervals. Or, more often, an object might have mostly unchanging features with little useful information content. But, say 10% of it contains volatile features of great interest. This has led to the use of non-uniformly spaced, or irregularly patterned, data samples.

[0008] For instance, a signal might exhibit rapid change (or volatility) in one concentrated region, but otherwise be fairly smooth. Or, an image might exhibit sharp edges as well as regions of unchanging, or predictable, patterns. In both cases, many sample points are required in certain concentrated intervals or areas where much information can be gleaned from the complex changes occurring there, while few sample points are needed elsewhere where little useful information can be extracted.

[0009] The number of sample points required using the standard uniformly spaced samples rises dramatically as the number of dimensions increases. Thus if a one dimensional signal requires O(N) samples, then a two dimensional image would need O(N²) samples, a three dimensional O(N³) 3 samples and so on. (A quantity Q(x) is said to be of order N^k, written O(N^k), if Q(x)/x^kis bounded.

[0010] What is required is a method to specify the manner of sampling an object wherein each the totality of samples provides comprehensive and non-redundant information about the object, while minimizing the number of samples required, especially as the dimension of the sampled object increases.

[0011] The following methods are generally known in the prior art:

[0012] Accept the limitations of uniform sampling and utilize Fourier Transform approximation techniques like the Discrete Fourier transform and its FFT implementation;

[0013] Approximate, generally via interpolation, non-uniformly spaced samples by uniformly spaced samples of sufficiently high resolution in critical regions tailored to the particular object;

[0014] Using pseudo-random numbers, generate an irregular scale for measurements. A Monte Carlo method is typically used to implement this;

[0015] Use multi-resolution wavelets; and

[0016] Adaptively sample data.

[0017] The Fourier Transform plays a dual role in N dimensional data recording. First consider a one dimensional signal representing time variations of some variable of interest, say voltage or price. The Fourier expansion of a signal is capable of representing an analog signal by a sequence of complex numbered Fourier coefficients, the amplitudes and phases of the Fourier components, whereby the content of a signal with a bounded frequency spectrum is accurately represented by a discrete data series. This itself results in effective data compression, especially for signal having significant periodic components. Secondly, the determination of the Fourier coefficients of the signal is determined by sampling of the signal, which brings one back to the need for the present invention.

[0018] The standard technique is to employ a Discrete Fourier Transform (DFT) as a substitute for the Fourier Transform. Furthermore the FFT or Fast Fourier Transform has been developed, which is an efficient implementation of the DFT. The DFT, therefore is an approximation that inherently distorts the spectrum generated by the Fourier transform. Nevertheless, it is tolerated as the most commonly used digital and image processing technique because it is analytically tractable and readily implemented on diverse hardware/software platforms.

[0019] The approximation, generally via interpolation, of non-uniformly spaced samples by uniformly spaced samples of sufficiently high resolution also distorts the calculated spectrum. But it is also potentially wasteful of sample points since it requires uniformly spaced samples as anchors around which irregularly spaced samples are interpolated. It is also difficult to implement its mathematics. Indeed, the overwhelming general practice is still to use the first approach.

[0020] The use of pseudo-random numbers is analytically or mathematically difficult to describe, if not pragmatically intractable. That is, it is not readily describable mathematically, say as a set of equations. A mathematical description is generally needed in practice, to develop algorithms that automatically map sample points to the irregularly spaced random points. It is also computationally complex because the deltas, or differences between each sample point and neighboring random points, generally must be recalculated because the random points used might not be reproducible.

[0021] Another, relatively new, general technique is multi-resolution analysis. This involves splitting the stream of possibly non-uniformly spaced data into streams of different resolution. In some circumstances, this technique is also referred to as “decimation”. Wavelets also use a multi-resolution analysis. The use of multi-resolution wavelets tries to solve the problem of mathematically decomposing non-stationary or highly irregular signals and images, into simpler components on multi-tiered scales. But the approach is generally difficult to implement in practice. Although there have been claims for very specific applications, its effectiveness vis-a-via standard DFT based techniques for general applications is not apparent. Indeed, the 1965 Cooley-Tukey FFT algorithm immediately paved the way for the revolution in medical imaging that accelerated in the 1970s. Such a similar impact is not apparent, even more than 13 years after the framework for wavelets was articulated in a coherent manner.

[0022] There were other attempts preceding wavelets, at calculating the Fourier transform of non-stationary objects, notably the Gabor Transform. But again, these techniques were not general enough to be of use off-the-shelf, as the FFT. Indeed, a major motivation underlying wavelets was to supplant the Gabor transform method with a multi-resolution approach.

[0023] Adaptive sampling has been used in many specialized applications in connection with the FFT. Its use of specialized and sometimes customized sampling algorithms (to decide on how to adapt to data changes) limits its usefulness as a general purpose tool. It requires much a priori knowledge of the specific data set being sampled.

[0024] U.S. Pat. No. 5,708,432, appears to disclose a sampling digitizer with a sampling master clock divided by a factor k which is derived from prime integers. The '432 patent is concerned with sampling a periodic function, where the number of samples desired during each cycle exceeds the response ability of the sampling system. A technique is used in which fewer samples are taken from successive cycles of the signal and are combined in order to simulate the effect of many samples taken during one cycle. To do this, a clock is divided by an integer K which is varied according to an algorithm. N consecutive samples are taken over M cycles of the signal. Coherent signaling is achieved when the signal frequency and the divided clock frequency are relatively prime. The particular algorithm finds the relatively prime fraction by selection from a Farey series of a given order. Unlike the present invention, the '432 invention selects a set of samples separated by constant interval. The use of relatively prime ratios to guarantee that the same points are not selected on a perfectly periodic signal is different from the present technique, which would work to sample even a pulse, and does not require that the signal be periodic.

[0025] An abstract in the prior art makes reference to sampling employing prime numbers. Proceedings the European Signal Processing Conference held in Trieste on Sep. 10-13, 1996.

[0026] The abstract states:

[0027] We demonstrate that it can be advantageous from a computational point of view to use a two-stage realization instead of a single-stage realization for sample rate conversions with prime numbers. One of the stages performs a conversion by a factor of two . . . . The other stage changes the sample rate by the rational factor N/2 . . . .

[0028] U.S. Pat. No. 5,243,343 uses variable time interval data sampling. See col. 5, line 18-19. U.S. Pat. No. 5,444,459 is essentially identical with the '343 patent.

[0029] Other prior art shows

[0030] 1. non-linear sampling without reference to the prime integer algorithm of the present invention (U.S. Pat. Nos. 4,188,583; 4,142,146; 4,903,021;and 5,229,668 [FIG. 3])

[0031] 2. the use of prime numbers in Hartley transform processors (U.S. Pat. No. 4,062,060; Boussakta et al., 24 Elect. Ltrs. No. 15, pp. 926-28 (Jul. 21, 1988); 136 IEE Proc. Pt. G, No. 5 (October 1989) pp. 269-77.)

[0032] 3. Other patents of interest include U.S. Pat. No. 5,712,809 entitled “Method and apparatus for performing fast reduced coefficient discrete cosine transform”; U.S. Pat. No. 5,682,524 entitled “Databank system with method for efficiently storing non-uniform data records”; U.S. Pat. No. 5,497,152 entitled “Digital to digital conversion using non-uniform sample rates”; U.S. Pat. No. 4,271,500 entitled “Device for converting anon-uniformly sampled signal with short-time spectrum to a uniformly sampled signal”; U.S. Pat. No. 5,712,635 “Digital-to-analog conversion using non-uniform sample rates”; U.S. Pat. No. 4,999,799 entitled “Signal processing apparatus for generating a Fourier transform”; and U.S. Pat. No. 4,969,700 entitled “Computer aided holography”.

[0033] 4.

[0034] The invention relates to signal processing, primarily digital signal processing employing N dimensional data sampling required to extract information about the object for subsequent processing. In order to extract information from a signal it is sampled at discrete points and the sampled values are fed into calculation algorithms. For example the sampled values may be used to calculate the coefficients of a Fourier transform of the signal by known algorithmic techniques, thereby extracting useful information about the Fourier spectrum of the signal. Other transforms may be calculated using other kernel functions to generate Radon, Hartley, or discrete sine and cosine transforms. The discrete points at which sampling occurs may for example be discrete times for a signal that varies with time. As a specific example of the utility of the invention, it is utilized to process an N dimensional Fourier and related transform using the non-uniform samples of the D scale and FFTs. This decomposition of a Fourier transform into a series of FFTs will be detailed later.

[0035] This invention is concerned with the selection of the values of the discrete points at which sampling occurs. Since most signals are composed of a few periodic components, there is much redundancy when the sampling takes place at regular intervals. Sampling at regular intervals also leads to artifacting, i.e. false data caused by the confusion of an interval with its harmonics. On the other hand it is algorithmically simpler to sample at regular intervals and it is necessary to keep the calculations manageable.

[0036] This invention provides a method including algorithms to generate an irregularly spaced set of sampling points which still provides the practical calculation advantage of uniform sampling. This methodology can also be tuned or configured to generate an approximation of a truly (pseudo)random set of points. Further, it can be tailored such that the resolution or width between neighboring points in the set is always within pre-specified maximum and minimum bounds. That is useful when sampling different data sets, each with a different number of significant digits. For a one dimensional situation, i.e. where the signal is a function of one parameter such as time, the set of sampling points is composed as the union of n sets of points, each set being points separated by regular intervals of 1/p_nwhere p_nis a prime number. Each set is denoted a “scale” for prime number p_n. The union is therefore characterized by a set {P_n} of distinct prime numbers. The union is intrinsically non-periodic. And for even a small set of distinct primes the distribution of points in the union appears almost random. This union is used to determine where the signal is sampled. The points at which sampling takes place is not dependent on the values of the signal. If the signal cannot be accurately reconstructed from the sampled values one merely recurses and increases the set of distinct primes, which is always possible because of their infinite number. Alternately, there might be a-priori knowledge of an object's information bandwidth. Then the prime numbers used to generate the set of points can be pre-determined, such that the highest prime number just exceeds the known, effective, bandwidth. For example, audio applications generally require bandwidth ranges of 5-20 KHz. Therefore there is no need to use prime numbers of order 10⁵when this application domain only uses an effective bandwidth of order 10⁴.

[0037] The invention may be implemented in a real-time measurement operation by running a series of clocks ticking at intervals characterized by 1/p_nfor the n^thclock, and sampling the signal at each clock tick. For a more than one dimensional situation, where the signal is characterized by multiple parameters, each parameter is treated as in the one dimensional case, and the sampling space is treated as the Cartesian product of one in dimensional spaces. The invention relies upon the incommensurate nature of different primes to prevent the invention from degenerating into an inefficient sampling technique such as using regular intervals would provide.

[0038] Although the invention is characterized by the use of primes, there are other techniques that may be used to generate the sampling points. The key features are that the points are derived from the union of sets of points generated by simple calculations and the union avoids a periodic structure. Another variation of this sampling technique is to offset the points associated with each prime scale. The offsets might be random or also prime numbers. The intent would be to better capture phase information in data analysis.

[0039] An efficient method is presented for sampling and thereby enabling the recording of, or measurements of objects, which is successful regardless of the irregularity or periodicity of the object. The sampling occurs at values of an independent variable. For example, where the sampled object is an electrical signal, the independent variable could be time. Where the sampled object is a multidimensional physical object, the independent variable could be a spatial coordinate parametrizing the space in which the physical object is embedded. Collectively, the points or values of the independent variable at which samples are taken is termed a “scale”. This terminology is consistent with the usual association of a scale with the ruler used to measure quantities.

[0040] This new measurement scale is presented with special properties that makes it particularly effective for sampling non-uniform data. The new scale is contrasted with the familiar ruler scale as a point of reference. The new scale is based, fundamentally, on the existing mathematics of set theory, real analysis of the continuum, and number theory. That is, it applies theoretical knowledge from set theory, real analysis and number theory to solve practical problems related to data sampling. The theoretical framework provides the necessary rigorous foundation to guarantee that this new scale yields provably correct results in all applications.

[0041] The new scale, herein termed a D scale, is applicable to problems in signal processing, including digital signal processing (DSP), image processing, scientific and engineering computational applications, data acquisition and statistical data analysis, to mention a few application domains.

[0042] The D scale can also be configured or tailored to match the varying bandwidth, resolution or accuracy, and point distribution requirements of the widely disparate application domains mentioned above. Specifically, the minimum and maximum resolutions of any particular D scale can be guaranteed by an algorithm, discussed later, that determines which prime scales to include when a D scale is constructed. Further, the D scale can be constructed with prime scale components such that the distribution of its points is made somewhat random for certain applications. Or, the points can be distributed to accentuate the multi-frequency components of an application domain, resulting in more points towards the center of a unit interval, and tapering off at the edges. This flexibility in generating point distributions is effected by including more or fewer prime scales in a D scale.

[0043] The invention employs a method, or an algorithm, to perform nonuniform sampling. It can be embodied or implemented as a product (hardware based, software based or combined hardware/software) to implement nonuniform sampling.

[0044] Among possible commercial applications are the following:

[0045] Real-time Data Acquisition

[0046] The D scale, along with its sampling algorithm, improves the efficiency of real-time data acquisition. Data acquisition is at the heart of many processes in such diverse fields as seismology, oceanography, medical imaging and process control, to name a few. The D scale and its algorithm samples judiciously, say when volatile data is detected, and not when the data is uninteresting. In contrast with standard techniques that leads to reduced cost of data acquisition since fewer sample points are generally used. The effectiveness of the sampling process is improved since useful samples are generated, not the many extraneous samples of a uniformly spaced procedure.

[0047] Post-acquisition Data Analysis

[0048] The D scale can also be used to analyze data already generated by a data acquisition or sampling device. For example, the D scale is used to generate a closed form expression for the Fourier Transform (described later). Using the closed form expression, it is then possible to quantify the bias in not using points in regions that are not varying enough to warrant attention. Therefore, it can effectively filter the acquired data or samples to focus on the more significant data sections, thereby reducing the quantity of data that must be processed to characterize the object sampled. In this manner the D scale acts as a means for data compression. This has application in the analysis of large quantities of data produced by NMR or tomography devices, and reduces the burden of analyzing such data. The same is true of the large quantity of data produced by seismography measuring apparatus, as an example. Further, other application domains produce resonance shapes which might contain only 10% useful information content.

[0049] Transforming Irregular Data into Periodic Components

[0050] The D scale can be used to partition an original, irregularly patterned data set to create multiple components. Each component contains points which, although not periodically patterned, are all reciprocal multiples of the same prime number 1/p_k, which is what bounds them to the same prime scale. Then, because the points on each prime scale is bandlimited by p_k, the set of points on each prime scale can be decomposed via standard Fourier analysis into periodic components. Thus, the D scale enables a two stage decomposition. The first stage decomposes the original set of irregularly spaced points into sets of uniformly spaced and bandlimited points, each associated with a prime scale. Then the second stage decomposes each bandlimited decomposition into the usual periodic components of a Fourier series. Also, each partition has a different resolution (1/p_k)such that the components form a nested set of resolutions.

[0051] Characterization of Noisy data

[0052] The D scale can filter random noise into periodic components thereby providing a useful characterization of the data, which reduces the effect of random noise in the data.

[0053] FIG. 1 shows a bounded, 1D, single resolution prime number 3, D scale.

[0054] FIG. 2 shows an unbounded, 1D, single resolution prime number 3, D scale.

[0055] FIG. 3 shows a unit interval, 1D, multi-resolution D scale comprising a composite scale using the prime numbers 2, 3 and 5.

[0056] FIG. 4 shows an unbounded, 1D, multi-resolution composite D scale.

[0057] FIG. 5 depicts the steps of the process for calculating the Fourier transform spectrum of an object.

[0058] FIG. 6 depicts the 1 dimensional component of an N dimensional D Scale comprising potentially an unlimited number of prime resolution scales. The union of points on these scales creates an irregular pattern of points on the composite D Scale.

[0059] FIG. 7 depicts a one dimensional multi-resolution D scale bounded by sampling time. The spectrum's boundary is determined by the prime number of the resolution used. FIG. 7 shows this effect.

[0060] FIG. 8 depicts the reciprocity theorem.

[0061] FIG. 9 depicts a benchmark signal used to evaluate the efficiency of the invention.

[0062] FIG. 10 depicts the exact spectrum of the benchmark signal of FIG. 9.

[0063] FIG. 11adepicts the amplitude spectrum of the benchmark signal arrived at by using Mathematica's Fourier function.

[0064] FIG. 11bdepicts the amplitude spectrum as a function of ω, across its entire range.

[0065] FIG. 12 depicts the amplitude plot of the spectrum contribution from the D scale using prime 2.

[0066] FIG. 13 depicts the amplitude plot of the spectrum contribution from the D scale using prime 5.

[0067] FIG. 14 depicts the amplitude plot of the spectrum contribution from the D scale using prime 11.

[0068] FIG. 15 depicts the spectrum from D scale primes 11, 5, and 2 contributions.

[0069] FIG. 16 depicts the spectrum from D scale primes 29, 19, and 11 contributions.

[0070] The D measurement scale has three components; a set of scale points used-to select data samples, a set of pre-calculated distances between nearest neighbors and an algorithm that efficiently associates or maps sample points to the scale points. The rigorous mathematical underpinning of the D scale involves two fundamental mathematical concepts; the notion of a one-to-one correspondence between points in two sets and that of the measure of a set. The measure of a set corresponds roughly to the inherent concept of size, e.g. length or area, but is defined precisely to cover sets whose inherent size is not obvious. The Lebesgue measure of a one dimensional set is a familiar example of the use of the concept of measure, and permits assigning a measure for such oddities as the set of all irrational points, etc.

[0071] The familiar ruler uses a subset of rational numbers (i.e., numbers that can be expressed as fractions, typically multiples of some high inverse power of two, which defines the resolution of the scale) for measurements. The rulers sampling algorithm finds the multiple of the inverse power of two that is nearest to the sample point. For example, an irrational number is approximated by a rational number of the required resolution.

[0072] The D scale satisfies several criteria for the utility of a measurement scale. It has

[0073] a. A degree of regularity, or of a uniform pattern, of the scale points for ease-of implementation, that is to say the definition of the scale points is simple to calculate,

[0074] b. A degree of irregular pattern of scale points to accommodate clustering of data points, that is to say the scale points exhibit a highly non-uniform distribution,

[0075] c. Multi-resolution capability to accommodate different orders of magnitudes in sample point resolutions,

[0076] d. Minimal number of scale points needed for effective measurement.

[0077] Note that the familiar ruler meets criteria a and c, providing uniformly spaced markings and multi-resolution as powers of 10. In contrast, a scale whose points are randomly selected only meets criteria b above. The D Scale meets all four criteria indicated above, therefore making it most effective for nonuniform data sampling.

[0078] Working Description of D Scale Invention

[0079] The D measurement scale is described both analytically and geometrically. Then, rigorous mathematical proofs are presented on the correctness of its operation. After that, the scale's sampling technique is described and its utility proven. Finally, useful and commercially viable examples are provided illustrating uses of the D scale.

[0080] The N-Dimensional D scale consists of selected points on the N dimensional rational number coordinate space (or grid). This scale will be shown to meet each of the criteria set forth for the D scale. It will also be shown to be effective in sampling and/or processing non-uniformly spaced data.

[0081] The set of measurement points in a D scale is described in stages. First, the one dimension (1D) scale is described. This is then generalized to a description of the N dimensional scale.

[0082] The 1D measurement scale is multi-resolutional, i.e. it is comprised of multiple scales of different resolutions. The 1D, unit interval scale with a single resolution is first described. Then, the unbounded 1 D scale with a single resolution is described. That is followed by the multi-resolution, unbounded, 1D scale.

[0083] Consider a one dimensional bounded single resolution scale, e.g. the set of points on the 1D, half-closed interval [0, 1), that are multiples of some prime number. Or more formally, consider the set of points S_{pio, 1)}={x ε [0, 1): x=n/p, n=1..p−1, pε {prime numbers}}, which is depicted in FIG. 1 for the prime number 3.

[0084] Next, extend this definition to an unbounded set of points by replicating the set of points on [0,1) across the entire number line. More formally, consider the set of points S_p={x ε (−∞,∞):x=j+n/p, j=−∞ . . . ∞, n=1..p−1, p ε {prime numbers}}. S_pis referred to as the p^thprime resolution D Scale. This is depicted in FIG. 2.

[0085] Next consider a multi-resolution bounded D scale. Consider N 1D, unit interval, prime resolution scales, each corresponding to a different prime number. A multi-resolution D Scale is then the superposition or union of measurement points from each of these N prime resolution scales. Formally, using the previous definition of S_p[0,1), the 1D multi-resolution D Scale is: S_l[0,1)={x ε ∪S_p[0,1):p=p_l. . . p_n, p_iε {prime numbers)), or equivalently, S_l[0,1)=∪S_p[0,1)where p=p_l. . . p_n, p_iε {prime numbers)). This is shown in FIG. 3.

[0086] Next, consider N 1D, unbounded, prime resolution scales, each corresponding to a different prime number. A multi-resolution D Scale is then the superposition or union of measurement points from each of these N prime resolution scales. Formally, if S_pis as defined previously, then the 1D multi-resolution D Scale is S_l={x ε ∪S_p:p=p_l. . . p_n, p_iε (prime numbers }), or equivalently, S_l=∪S_pwhere p=p_l. . . p_n, p_iε (prime numbers }). This is shown geometrically in FIG. 4.

[0087] Summarizing, each measurement point in the multi-resolution 1D D Scale is located by its three coordinates (p,i,j), where p is the prime resolution scale in which the point is located, i is the integral interval on the prime resolution scale p, and j is the fraction of the prime number in that i^thinterval. Therefore, if a measurement point's location on the rational number based ruler is d, then the relation of these scales is d=i+j/p.

[0088] Progressing now to multidimensional examples, a single resolution, n dimensional, D scale is defined for the bounded unit hypercube. First, consider Cartesian cross product of the measurement points on the one dimensional half-closed interval [0, 1), namely, S_{pN[0, 1)}=Π S_p[0,1)=S_p[0,1)× . . . ×S_p[0,1)(N times). An N dimensional point on the unit hypercube D Scale is then described as X ε Π S_p[0,1)→X=(x_l. . . x_N)→x_jε S_p[0,1).

[0089] This has also been extended to the unbounded hypercube. First, consider the measurement points on the one dimensional half-closed interval [0, 1), previously defined as S_{p[0, 1)}. Define the Cartesian cross-product of N of these 1 dimensional unbounded scales as follows: S_pN=Π S_p=S_p× . . . ×S_p(N times). An N dimensional point on the unbounded unit hypercube D Scale is then described as X ε S_pN→X ε Π S_pX=(x_l. . . x_N)→x_jε S_p.

[0090] This in turn is expandable into a multi-resolution, N dimensional, D Scale. For example, an N-dimensional multi-resolution scale on a unit hypercube is created by the union of N dimensional, unit hypercube, prime resolution scales, for each prime number included in the scale. Formally, the N dimensional, multi-resolution, D Scale is S_N[0,1)={X ε ∪S_pN[0,1):p=p_l. . . p_n, p_iε {prime numbers}}, or equivalently, S_N[0,1)=∪S_pN[0,1)where p=p_l. . . p_n, p_iε {prime numbers}}, or equivalently, S_N[0,1)=(x_l. . . x_N):x_iε ∪S_p[0,1)where p=p_l. . . p_n, p_iε {prime numbers}}.

[0091] An N-dimensional multi-resolution scale on an unbounded hypercube is created by the union of N dimensional, unbounded, prime resolution scales, for each prime number included in the scale. Formally, the N dimensional, multi-resolution, D scale is: S_N={X ε ∪S_pN:p=p_l. . . p_n, p_iε {prime numbers}}, or equivalently S_N=∪S_pN:p=p_l. . . p_n, p_iε {prime numbers}, or equivalently, S_N={(x_l. . . x_N):x_jε ∪S_p:p=p_l. . . p_n, p_iε {prime numbers}.

[0092] I have accomplished rigorous proofs of three theorems that demonstrate properties establishing the utility of the D scale. The first property is the scale points' uniqueness. The second is the scale points' ability to completely cover a line (in 1dimension) or region (in N dimensions) without overlap. And the third property is the scale points' capability to provide any desired resolution through nesting of prime resolution scales. These properties are direct consequences of the intrinsic properties of prime numbers. Indeed, that is why only prime numbers are included in the D Scale. But more generally, any other set of numbers should have these characteristics as well. Specifically, using random or prime numbered offsets should be sure to obey these theorems. These theorems underlie the validity of the scale and provide the necessary confidence that the scale always performs its measurement function correctly.

[0093] Theorem 1 is the Fundamental Uniqueness Theorem

[0094] Each D Scale point is in one and only one prime resolution scale. That is, for any pair of N-th dimensional prime resolution scales S_{i, N}and S_j,N, in the D Scale it follows that S_{i, N}∩ S_j,N=ø where I, j ε {prime numbers}.

[0095] Theorem 2 (One Dimension)

[0096] The set of D Scale point neighborhoods covers the number line (−∞,∞) (i.e., R1). And the covering does not contain any overlapping neighborhoods: ∪Np_i,j=(−∞,∞).

[0097] Corollary 1 (1 Dimension): Any sample point can be assigned or mapped to one and only one D Scale point.

[0098] Theorem 2 (N Dimensional)

[0099] The set of D Scale N dimensional point neighborhoods covers the R^Ncontinuum. And the covering does not contain any overlapping neighborhoods. ∪_Npi,j=∞^N. (These neighborhoods are defined below.)

[0100] Corollary 1 (N Dimensions): Each sample point can be assigned to one and only one point neighborhood in the D to Scale.

[0101] Definition: The resolution (or span) of a D Scale is the maximal distance between adjacent points such that they are distinguishable. The gaps between adjacent D Scale points varies. But there are applications where the maximal gap between adjacent points must be guaranteed to within a specified resolution.

[0102] Theorem 3 (Infinitely Scalable Resolution)

[0103] The span of the D Scale decreases to 0 as the prime resolution tends to ∞. Conversely, the prime resolution tends to ∞ as the span decreases to 0. This theorem guarantees that the resolution of a D Scale can be adjusted to as small a value as required, by incorporating higher prime resolution scales.

[0104] Corollary 4: A prime resolution scale can always be found to implement a specified maximum gap on an N dimensional D scale. Additionally, a theorem is proved Theorem 4: A necessary condition for neighboring scale points to always have a minimum interval, ε, between them, is that:

[0105] |1/p_i−1/p_i+1|<ε

[0106] for all p_iprime scales.

[0107] This is readily proven using mathematical induction, since:

[0108] |1/p₀−1/p₁|<ε→|n/p₀−m/p₁|<εn<p₀, m<p₁.

[0109] If

[0110] |1/p_i−1/p_i+1|<ε→|n/p_k−m/p_l|<εn<p_k, m<_pl

[0111] Then:

[0112] |1/p_i+1−1/p_i+2|<ε→|1/p_i+1−1/p_i+2|<ε→|n/p_k−m/p_l|<εn<p_k, m<p_l

[0113] by the induction condition. And

[0114] |1/p_i+1−1/p_i+2|<ε→|n/p_i+1−m/p_i+2|n<p_i+1, m<p_i

[0115] Therefore |n/p_k+1−m/p_l|<ε n<p_k+1, m<p_ll<k+1. Therefore the induction hypothesis is true for case k+1.

[0116] Therefore a D scale can be constructed following an induction argument. Start with the lowest primes that satisfy the constraint. Then add a third prime that satisfies the constraint with the second prime, and onward. This constructs a D scale layer by layer, or prime subscale by subscale.

[0117] Non-Uniform Sampling of Data Using the D Scale sampling can be performed in real-time, as each data point is individually received. Or sampling can be performed after a collection of data points is first gathered. In either case, an algorithm is presented to “measure” or sample data points using the D Scale.

[0118] This algorithm is first presented for 1 dimensional data sets and then generalized to N dimensional samples. The algorithm is presented as a theorem, with the proof of its utility relegated to the examples.

[0119] Theorem 5 (Sample Measurement Theorem in 1D)

[0120] A sample point x can always be assigned, or mapped, to one and only one point in a D Scale according to the following algorithm:

[0121] Start from the lowest prime resolution scale in the D Scale. Determine if the sample point is in the neighborhood of any of the points in the D Scale. If it is, then assign that scale point to it. Otherwise, continue to the next prime resolution and perform the same search for nearest measurement point. The sample point must eventually find itself in the neighborhood of some point in the D Scale. Formally:

[0122] for (P_iwhere P_i=p_l. . . p_n, p_iε {prime numbers in this D Scale})

[0123] Consider a sample point x, expressed as a decimal. Then x=a.b, where a is an integral part of the decimal, and b is the fractional part. The integral part of this sample, a, is the a-th interval on the prime resolution scale P_i. The fractional part, b, is used to determine which multiple of 1/P_ito use, as follows. Consider: (.b)*P_i

[0124] The floor function extracts the greatest integer j less than or equal to f: j=└(.b)*Pi┘

[0125] Therefore: j/P_i=<.b

[0126] Also: .b=<(j+1)/P_i

[0127] Otherwise, j+1=└(.b)*Pi┘ which contradicts the above.

[0128] Thus: j/P_i=<.b=<(j+1)/P_i

[0129] The sample point's remaining distance from (j/P_i) is:

[0130] ε=b−(j/P_i)

[0131] As verification, the sample point can be located on the a^thinterval of the j^thelement of prime resolution scale P_i. Or:

[0132] x=a+(j/P_i)+ε==a.b

[0133] Next, determine if the sample point is in the neighborhood of either i th or j+1 st measurement points of prime resolution P_i; P_i,jor P_i,j+1respectively.

[0134] Let δ_j,R=distance of P_i,jto its nearest right neighbor among the D Scale measurement points. δ_j+1,L=distance of P_i,j+1to its nearest left neighbor among the D Scale measurement points.

[0135] If (ε=<δ_j,R/2)

[0136] then measurement point P_i,jis selected as that sample's nearest approximation. That D Scale point is used in subsequent calculations to represent the original sample point.

[0137] Else if (ε=<δ_j+1,L/2)

[0138] then measurement point P_i,j+1is selected.

[0139] Else continue the same logic, using the next prime resolution scale.

[0140] Theorem 4 (Sample Measurement Theorem in N Dimensions)

[0141] An N dimensional sample point can always be assigned, or mapped, to one and only one N dimensional point in a D Scale according to the following algorithm:

[0142] Start from the lowest prime resolution scale in the D Scale. Determine if the sample point is in the neighborhood of any measurement points in the D Scale by examining the 1D neighborhoods of each of its coordinates, at the same resolution scale. If it is in the neighborhood of a scale point, then assign that scale point to it. Otherwise, continue to the next prime resolution and perform the same search for nearest measurement point. The sample point must eventually find itself in the neighborhood of some point in the D Scale.

[0143] Formally:

[0144] For (P_iwhere P_i=p_l, . . . p_n, p_iε {prime numbers in this D Scale})

[0145] Consider an N dimensional sample point X, expressed as:

[0146] X=(x1, x2, . . . x_N)

[0147] Then apply the 1D algorithm to each coordinate x_ito that coordinate point's nearest 1D scale point.

[0148] Coordinates at different dimensions x_lwill find their neighborhood at different prime resolutions x_pk.

[0149] Eventually each coordinate must find its neighborhood.

[0150] The sample point will then be assigned to the scale point as:

[0151] X=(x1, x2, . . . x_N)→ (X_p1, X_p2. . . X_pN)

[0152] General Properties of the D Scale

[0153] FIGS. 6 and 7 summarize the key properties of the N dimensional D Scale, even as the 1 dimensional case is highlighted. One property is that the scale is composed of potentially an unlimited number of prime resolution scales each of which has points that are uniformly spaced. The union of points on these scales creates an irregular pattern of points on the composite D Scale . The availability of a countable infinite number of prime resolution scales leads to any desired resolution. Further, the scale can be tailored to the application domain of interest such as voice or images. Thus, the resolution between neighboring points can be made bounded, or the point pattern generated to be more or less random.

[0154] Another property is that each scale point has a neighborhood about it that does not overlap with the neighborhood of any other scale point. This property leads to a measurement algorithm that is guaranteed to identify the neighborhood within which a sample point resides, in order to map it to a point on the scale. FIG. 6 illustrates this effect. Thus, it is not possible for a sample point to “fall-through the cracks” and not map to any scale point.

[0155] Yet another property is that the spectrum of an object's (e.g., signal, image) sample points on any one prime scale is necessarily band limited. The spectrum's boundary is determined by the prime number of the resolution used. FIG. 7 shows this effect.

[0156] Contrasting this scale with the standard ruler, it requires far fewer measurement points than a ruler of the same resolution or depth. Thus, at any depth measured by the number of decimal places (and therefore resolution), there are fewer ruler points per unit interval than with the standard ruler.

[0157] Example of Implementation of the D Scale

[0158] The D Scale can be completely implemented in software, in hardware or in both software and hardware.

[0159] Consider the case of post-acquisition data analysis mentioned earlier. Then, the scale, with its associated mapping, can be implemented, along with sampling, entirely in software, on already acquired data. Such is commonly the case in medical imaging or seismology, to touch on a few general application domains. Or, consider the case of real-time data acquisition, say in radar. Then it can be implemented entirely in hardware such that the sampling is automatically performed in real-time, as data is received by a data acquisition device. Or the scale's points can be mapped to the D scale in software, while its associated sampling algorithm is implemented in hardware.

[0160] The example includes:

[0161] 1. A set of prime numbers used to comprise the scale's multi-resolutions. The set is selected with the aim to bound the range of resolutions, according to an application domain constraint, with as few points as needed.

[0162] 2. Set of points in the scale as the union of points in each of the scale's resolutions.

[0163] 3. Distances between adjacent points on scale. It is saved in-memory or elsewhere.

[0164] 4. Sampling algorithm.

[0165] Calculating the N Dimensional Fourier Transform Using the D Scale

[0166] The underlying mathematical foundation of many diverse applications includes the Fourier Transform. These applications are in such diverse domains as signal processing to include digital signal processing (DSP), image processing, scientific and engineering computational applications, data acquisition and statistical data analysis. But there are several difficulties in actually applying the Fourier Transform directly to measurement samples.

[0167] Therefore in practice, several approximations of it are routinely used instead. They include the Discrete Fourier Transform which is typically calculated using the FFT family of algorithms, the Discrete Cosine Transform and others. There have been many attempts at calculating the Fourier transform of a function using irregular samples.

[0168] The problem is then to calculate an accurate spectrum approximation of the following Fourier transform (using the time-frequency analog):

[0169] F(ω_m)=1/(2π) ∫h(t) exp (−iω_mt)dt

[0170] given a non-uniformly spaced set of time samples.

[0171] Part of the problem of calculating the transform is to provide feedback on how good the approximation is, given the set of sample points applied. The algorithm also can provide guidance on whether a redeployment of those sample points provide a more accurate spectrum. Further, it can be used to determine whether more sample points are needed to gain noticeably increased accuracy, or whether the point of diminishing returns has been reached.

[0172] The present invention is used to provide a closed form set of equations, to calculate the Fourier transform of a complex valued function given a (finite or infinite) set of non-uniformly spaced sample points. The algorithm can also be applied without change to its procedural structure, to calculate similar transforms, such as, but not inclusively, the Harley transform, Discrete Cosine Transform (DCT), Fourier series. The algorithm is fast because it expresses the Fourier transform as a series of Fast Fourier Transforms (FFT).

[0173] The Fourier transform is calculated by starting with first principles, from the definition of a Riemann Integral of a multi-dimensional function. The Riemann partition used to approximate the integral (say as an area under a curve), is determined by the set of measurement points in a time domain D scale. Also, the set of frequency components calculated is determined by a set of measurement points on a frequency domain D scale. Finally, the Shannon Sampling theorem is applied to calculate the approximate time domain ,function that corresponds to the just calculated frequency spectrum.

[0174] The Fourier Transform of a function is expressed as a sum of band limited functions. Each of these band limited functions is derived from a subset of points in the D scale associated with one of its prime scales. Because of that, each of the bandlimited functions can be calculated efficiently using the standard FFT. Therefore, the Fourier Transform can be calculated with the same order of efficiency as FFT. Further, because each function component is bandlimited, each one can be calculated using the Shannon Sampling Theorem (which relies on the Nyquist Sampling Theorem to reconstruct an original object from its bandlimited sample points). Therefore, an approximation of the original object can be reconstructed from the set of all components, providing the necessary feedback on the degree of likeness to the original object. Then, the sampling can be adjusted according to the desired accuracy, by changing the sample point deployment and possibly adding more samples.

[0175] The steps of the process are depicted in FIG. 5 illustrating the iterative steps to approximate the Fourier Transform of an object (in N dimensions), using the sample points available.

[0176] This Fourier Transform calculation can only meaningfully calculate as many points in the frequency spectrum as are available sample points. That is, it is not possible to generate more information than what was provided. This constraint is readily expressed by the following matrix equation:

[0177] Frequency vector=Fourier Decomposition Matrix*Time Sample Vector[ω₀. . . ω_N−1]=M×[t₀. . . t_N−1]

[0178] This equation shows how the set of derived points in a spectrum is structurally related to the sampled point inputs. The detailed nature of this matrix structure follows. The points in the frequency vector are determined according to the Reciprocity theorem. This theorem states that:

[0179] If A is the scale bound in one domain (typically time), Ω is the scale bound in a second domain (typically frequency), and N is the total number of samples in the first domain, then:

[0180] A*Ω=N.

[0181] An alternate formulation is:

[0182] If Δt is the sampling interval in one domain (typically time), Δf is the interval in a second domain, (typically frequency), and N is the total number of samples in the first (time) domain, then

[0183] Δt*Δf=1/N

[0184] Consider sampling a signal on prime scale p, for T_ptime intervals, where T_pis a prime number. Then the boundary is T_p. Then Δt=1/p where p is the subscale's prime sampling rate and the number of samples N=T_p*p on subscale p. On the frequency side the roles are reversed. Δf=1/T_pon subscale T_p. Let p represent the sampling rate of a D Scale's prime subscale.

[0185] This structure satisfies the reciprocity theorem.

[0186] Δt*Δf→1/p*1/T_p=1/(T_p*p)=1/N

[0187] The set of derived frequencies is then related by the reciprocity theorem to the sample points by a one-to-one correspondence between time subscales and corresponding frequency subscales.

[0188] These calculations then yield an expression for the Fourier Transform decomposition for each frequency subscale characterized by T_pifor each corresponding frequency subscale:

[0189] F(2πm/T_pi)=Real(2πm/T_pi)+iImaginary(2πm/T_pi)

[0190] or

[0191] F(2πm/T_pi)=|F(2πm/T_pi)| exp (−iθ)

[0192] Where:

[0193] |F(2πm/T_pi)|=SQRT[Imaginary(2πm/T_pi)²+Real(2πm/T_pi)²]−Amplitude

[0194] θ(2πm/T_pi)=Arctan(Imaginary(2πm/T_pi)/Real(2πm/T_pi))−Phase

[0195] IF: $F_n\left({\omega }_m\right)=\frac{1}{2}\pi {\int }_{-\infty }^{+\infty }h\left(t\right)e^{-i{\omega }_{m^t}}dt=\frac{1}{2}\pi {\int }_0^{T_{\rho \max }}h\left(t\right)e^{-{\mathrm{i\omega }}_{m^t}}dt$

[0196] THEN: $\left|F\left(2\pi m/T_{p_i}\right)\right|=\frac{1}{2}\pi \sum _{k=\left\{\min \mathrm{prime}\right\}}^{\left\{\max \mathrm{prime}\right\}}\sum _{n=0}^{n=\left(T_{\mathrm{pk}}*p_k\right)-1}{\delta }_{pk,n\%\mathrm{pk}}h\left(n/p_k\right)e^{-\mathrm{i2\pi }\left(\frac{T_{\mathrm{pk}}}{T_{\mathrm{pi}}}\right)\mathrm{nm}/T_{\mathrm{pk}}*p_k}$

[0197] where:

[0198] The prime scales associated with prime numbers p_iare components of a D Scale

[0199] T_{p i}is the observation length, in units of time, on prime scale p_i.

[0200] T_{p i}must be a prime number. T_pmaxis the maximum value.

[0201] m=0, . . . T_pi−1 for each frequency prime scale T_pi

[0202] where k=min to max prime scales.

[0203] δ_{pk,n %pk}is the interval between a point on the D Scale located by (p_k,n %p_k) and its next, right, nearest neighbor on the same D Scale. (% is an operator extracting the non-integral part of a fraction.)

[0204] {max prime} has a greater bandwidth than the function being transformed. Σ′ is a sum over the prime numbers in the scale.

[0205] If the original signal, h(t)is expressed as a complex function h(t)=h_r(t)+i h_i(t)

[0206] then:

[0207] {max prime} n=(T_pk*p_k)−1

[0208] Real(2πm/T_pi)={fraction (1/2)}πΣ Σδ_{pk,n %pk}h_r(n/p_k) cos (2π(T_pk/T_pi)nm/T_pk*p_k)+δ_{pk,n %pk}h_i(n/p_k)sin

[0209] (2π(T_pk/T_pi)nm/T_pk*p_k)

[0210] k={min prime} n=0

[0211] Imaginary(2πm/T_pi)={fraction (1/2)}πΣ Σδ_{pk,n %pk}h_r(n/p_K) cos (2π(T_pk/T_pi)nm/T_pk*p_k)+δ_{pk,n %pk}h_i(n/p_k)sin

[0212] (2π(T_pk/T_pi)nm/T_pk*p_k)

[0213] Discussion of Closed Form Solution

[0214] Consider the innermost terms in the Fourier Transform Decomposition:

[0215] n=T_pk*p_k−1

[0216] Σ δ_{pk,n %pk}h(n/p_k)e^−(pkpipk_k

[0217] n=0

[0218] Each such term has the form of a Discrete Fourier Transform. Therefore the efficient Fast Fourier Transform (FFT) algorithms can be applied to each of these terms. Note that there is an additional term in the exponent, (Tpk/Tpi), which is integral to the invention. But this is a constant for the series and does not change the FFT structure.

[0219] Another point to consider is that although the outer summation appears to be unbounded, it is bounded in any physical situation. That is, only a finite number of prime scales need to be included in the D Scale, because no physical process can have an infinite bandwidth. Practical constraints also dictate realistic bandwidths as well. Thus, a voice application does not need generally need to consider a base bandwidth of higher than 20 KHz.

[0220] Yet another point is that this equation quantifies just how much is thrown out by the standard approximation of a Fourier Transform by a Discrete Fourier transform. Indeed, all but one prime scale is disregarded in the usual approximations. Therefore, the noticeable accuracy that must result from including the multi-resolutions will have real impact on the many commercial applications that are based on the Fourier Transform.

[0221] Finally, because this decomposition is derived from first principles, it preserves the phase as well as amplitude of the original Fourier transform.

[0222] There are many practical ways to sample data using the D Scale, and then to calculate the Fourier transform. The varied approaches depend on operating circumstances, such as real-time vs. post-data acquisition analysis, or available a-priori knowledge or not, as well as other cases. An embodiment is given below to highlight one manner by which the D Scale is used to calculate the Fourier Transform of sampled data during signal processing.

[0223] Practical Example of Processing a One Dimensional Signal

[0224] An example is presented that implements the use of the D scale in signal processing. Although the invention is applicable to multidimensional signals and images, etc. from many fields, a preferred embodiment will be described for a weak, sudden, signal over a telephone wire pair having a small bandwidth.

[0225] The typical signal of this type has an exponential decay. This particular shape has a spectrum that is directly derivable from its Fourier transform. Its exact spectrum serves as a benchmark to compare the relative accuracy of the usual, textbook, Discrete Fourier Transform approximation against this invention.

[0226] Two results emerge from this example. One is that the usual, textbook, approach requires many more uniformly spaced samples, yet still does not generate an accurate replica of the benchmark spectrum of the sound signal. Second, this invention uses far fewer points and does yield an accurate spectrum, comparable to the original.

[0227] This benchmark then validates this invention's methodology and decomposition equation.

[0228] FIG. 9 shows the signal whose spectrum will be extracted. In order to demonstrate the effectiveness of the invention the benchmark signal was reproduced in a closed analytic form using Mathematica.

[0229] FIGS. 10a-bare the mathematically exact spectrum of FIG. 9 to which will be compared the methodology of the prior art and of the invention.

[0230] Result of Prior Art Approximation of Signal

[0231] The Fourier Transform cannot be generally applied directly to data sets, uniform or not. One barrier is that data sets cannot, in the vast number of cases, be expressed in a mathematically tractable form. More important, the exponential term in the Fourier Transform causes partial to complete cancellation of terms generated from uniform sampling, leading to inaccurate spectrum calculations.

[0232] These limitations notwithstanding, the usual approach is to approximate the Fourier Transform by a Discrete Fourier Transform (DFT) on uniformly spaced data. However, if the uniform rate is lower than the maximum possible bandwidth, then vital information could be lost in the calculated spectrum. Therefore the uniform sampling rate is set to the Nyquist bandwidth. This employs only one term of many that the invention will utilize: $F\left(2\pi m/T_{\mathrm{pi}}\right)\approx \frac{1}{2}\pi \sum _{n=0}^{T_{\mathrm{pk}}*p_k}h\left(n/p_k\right)e^{-\mathrm{i2\pi }nm/T_{\mathrm{pk}}·p_k}$

[0233] The calculated spectrum is arrived at by using Mathematica's Fourier function. As may be seen from FIG. 11, this differs markedly from the exact spectrum in FIG. 10. Note, in particular, that the amplitude at 0 is much greater than the expected value of 2. Also note that the curvature in the lower left hand corner differs from the exact spectrum.

[0234] Description of the Non-Uniform Sampling Device

[0235] A measuring device is used to generate real-time, non-uniformly spaced samples, of a sudden sound over a telephone wire pair. The basic wire-pair telephony audio range has a 4 KHz baseband bandwidth. But the sound is weak enough that the effective base bandwidth is about 6 Hz. There is no loss in generalization, using this low frequency.

[0236] Non-uniformly spaced samples -are generated by a measuring device that dynamically changes its sampling rate based on the volatility of its most recent samples. That is, it generates more samples when recent samples suggest higher volatility, and far fewer samples where the time series is relatively constant. It does this by superimposing initial uniform sampling rates with new sampling rates that it activates when volatility increases. This provides a multi-resolution sampling effect. FIG. 4 shows how this device can be constructed using off-the-shelf components.

[0237] Invention's Spectral Determination of Non-Uniformly Sampled Input Signal

[0238] Consider the original signal as shown in FIG. 9. Its most volatile region is between 0-2 followed by 2-5 time units. Lastly, the region between 5-11 is relatively constant, without much information. therefore the non-uniform sampling device constructed above is programmed to partition the incoming data into three D scale subscales, at rates 1/11, 1/7, 1/3.

[0239] The initial rate is 1/3 for the entire 11 time units. The 1/7 rate is applied only until the volatility reaches a low enough threshold, or 5 time units. Finally, the highest sampling rate of 1/11 is only applied for the first 2 time units, where volatility is highest.

[0240] The D scale is constructed from the selected prime sampling rates, noted in the previous section. The first step creates the D scale's prime subscale vectors, in this case three. The Mathematica ‘N’ operation calculates decimal values from the fractions. These subscales are then joined and then most importantly, sorted from low to high numeric values. This step effectively mixes the values from the prime subscales. Then identical entries are removed. This occurs because each unit interval is a part of each scale, like 1,2, etc. Finally, as a Mathematica artifact, the last entry in the temporary, intermediate, scale is appended.

[0241] Do[sp[i]=N[table[(n−1)/p[[i]], {n, 1, p[[i]]pt[[i]]}]], {i, Length[p]}];

[0242] scaleTmp=Sort[Flatten[Join[Table[sp[i], {i, 1, length[p]}]]]];

[0243] scale=Table[Null, {0}];

[0244] Do[If[UnsameQ[scaleTmp[[i+1]], scaleTmp[[i]]]scale=Append[scale, scaleTmp[[i]]], Null], {i, 1, Length[scaleTmp]−1}];

[0245] scale=Append]scale, scaleTmp[[Length[scaleTmp]]]];

[0246] creates the D scale. The intervals or widths between adjacent D scale points is calculated:

[0247] delta=Table[scale[[i+1]]−scale[[i]], {i, 1, Length[scale]−1}];

[0248] The values of the original signal are calculated at the D scale points. Non uniform sampling at the scale points is determined by

[0249] h=Table[Exp [−.5 scale[[i]]] UnitStep[[i]]], {i, 1, Length[scale]−1}];

[0250] The D scale is indexed into its prime subscales. This is a programming artifact, not intrinsic to implementing the algorithm. The aim here is to effectively unmix the values in the D scale (due to the sort operation above), and place their D scale locations, or indexes, into their associated prime subscale. Each subscale vector is generated by first locating its members in the D scale. A member is identified if the product of this subscale's prime value and a scale value is an integer. Then, that member's D scale index is placed in the subscale array by:

[0251] Clear[b]

[0252] Do[

[0253] Do[

[0254] If[

[0255] IntegerQ[

[0256] Rationalize[

[0257] scale[[n]] p[[m]]

[0258] ]

[0259] ],

[0260] b[m]=Append[b[m], n],

[0261] Null

[0262] ], {n, Length[scale]}

[0263] ], {m, Length[p]}

[0264] ];

[0265] This step adjusts the index vectors created above. It removes those boundary scale entries like 1,2,etc. that were placed in each subscale. These scale points are assigned arbitrarily to one of the subscales only, in this case the first. Then it is removed from the other subscales. Again, this is a programming artifact and would not be necessary if another implementation approach were taken.

[0266] Finally, the implementation of Σ S is:

[0267] Do[f[1]=

[0268] Table[

[0269] {2Pi m/pt[[1]],

[0270] Abs[

[0271] Sum[

[0272] Sum[

[0273] h[[ c[k][[n]] ]]

[0274] delta[[ c[k][[n]] ]]

[0275] Exp[

[0276] −2Pi (pt[[k]]/pt[[l]])

[0277] (nm)/(p[[k]]pt[[k]]) I

[0278] ]

[0279] , {n, 1, Length[c]]−1}

[0280] ]

[0281] ]

[0282] }

[0283] , {m, 0, Length[c[1]]}

[0284] ]

[0285] , {l, 1, Length[p]}

[0286] ];

[0287] the decomposition equation is implemented, now that the data structures are in place. Following the structure of the implementation, a set of frequency amplitude vectors is generated, from the sets of sampled values in the D scale's prime subscales, again in this case three. The frequency vector is denoted f[l], where l is indexed from 1 to Length[p], or three in this case. The D scale values are fetched indirectly via the c[ ] vector created in a previous step. Thus, h[ c[k][[n]] ] is the value of the original signal ‘h’ at the scale point located by its location in prime subscale vector c[k], at position n in that vector.

[0288] The amplitudes of these newly generated frequency vectors f[l] are plotted, to highlight the underlying approach of generating each subscale of the spectrum's D scale. The Mathematica ‘N’ operation calculates decimal values from the complex numbers whose real and imaginary components are expressed as fractions.

[0289] Finally, the spectrum is plotted by juxtaposing the contributions from each of the spectrum D scale's prime subscales (each of which is shown in FIGS. 12-14).

[0290] Note that it accurately replicates the original spectrum from 0-10. Specifically, it has a value of 2 at 0, and its curvature is very similar to the original. Thus, the value of 0.5 is reached at about 1.75, as in the original. And the value of 1 is reached at about 1 on the x axis, again as in the original spectrum. See FIG. 15. Thus the result of the present invention is excellent compared to the prior art.

[0291] Finally, the spectrum was recalculated (see FIG. 16) using more data samples, at higher rates. The primes used were:

[0292] p={29, 19, 11} instead of p={11, 7, 3}

[0293] The sampling time intervals remained the same as in the prior case. Note that it also generates an authentic replica of the original spectrum. But, as expected, it is smoother than the previously generated one because more sample points are applied. The point at value 0 is not indicated, although it would be the correct value of 2, because the code used inadvertently started from the next to first index, not the first scale index. In practice the difference between a physical signal and the sampled transformed result may be calculated and from the difference it may be determined whether to recurse the sampling process while increasing the set of distinct primes.

[0294] Advantages of the D Scale

[0295] The advantages of using the Druck Scale are apparent when comparing the benchmark, textbook and invention calculation results.

[0296] e. The number of sample points required is dramatically smaller than using the textbook approach. For it required a high data sampling rate, at twice the maximum bandwidth, throughout the measurement period.

[0297] f. The results are much more accurate because it does not rely on the coarse Discrete Fourier Transform approximation. In fact, as the closed form expression indicates, the DFT approximation is really one of many DFT components of the Fourier Transform above. That is, the invention shows that the standard approximation throws out most of the relevant terms.

[0298] g. The invention is faster because fewer points are used in the calculation. And the invention shows how the calculation is comprised of decoupled terms, each with the FFT form. Therefore the FFT efficiency can be applied, as well as parallel processing optimizations.

[0299] h. The invention also enables the phase of the spectrum to be accurately calculated, because it relies on first principles. This is quite important for many applications, where phase is the determining factor in designed or observed behavior.

[0300] Although the invention has been described in view of particular embodiments, the scope of the invention should be determined by the following claims: