System and methods of nonuniform data sampling and data reconstruction in shift invariant and wavelet spaces

Abstract

A system and methods for converting data for an object of interest, which is characterized by a function f, between a digital form and an analog form. In one embodiment, the method includes the steps of locally selecting a plurality of data sample points in the order of N, N being an integer, wherein the data sample points are in a first form of data type, performing a transformation in a shift invariant space to the locally selected data sample points to obtain a new data set that is in a second form of data type different from the first form of data type, and reconstructing the object of interest f from the new data set. In one embodiment, the first form of data type is one of the digital form and the analog form, and the second form of data type is one of the digital form and the analog form that is different from the first form of data type. In other words, the first form of data type can be in digital form, and the corresponding second form of data type is in analog form, and vice versa.

Description

[0003] Some references, which may include patents, patent applications and various publications, are cited in a reference list and discussed in the description of this invention. The citation and/or discussion of such references is provided merely to clarify the description of the present invention and is not an admission that any such reference is “prior art” to the invention described herein. All references cited and discussed in this specification are incorporated herein by reference in their entireties and to the same extent as if each reference was individually incorporated by reference. In terms of notation, hereinafter, “[n]” represents the nth reference cited in the reference list. For example, [5] represents the 5th reference cited in the reference list, namely, Aldroubi and Grochenig, Nonuniform sampling and reconstruction in shift-invariant spaces, SIAM Rev. 43 (2001), No. 4, 585-620.

FIELD OF THE INVENTION

[0004] The present invention generally relates to a system and methods for data representation of an object of interest. More particularly, the present invention relates to a system and methods for data representation of an object of interest through nonuniform data sampling and data reconstruction in shift-invariant and wavelet spaces.

BACKGROUND OF THE INVENTION

[0005] Any object of interest in the world can be represented by corresponding information, or data. Information or data can then be processed so that people may have access to it. As known to people skilled in the art and as used in the specification throughout, among other things, the term data generally refers to numerical, graphical or alphabetical information representing one or more objects of interest. Data may be obtained or reduced through empirical observations or from instrument readings. The term data processing generally refers to various operations on data. In particular, among other functions, data processing includes data recording, transmission of data, manipulation of data, data conversion, data reporting, data reconstruction, storage of data, retrieval of data, and usage of data.

[0006] One aspect of data processing relates to digital data processing. Modern digital data processing of functions (or signals or images) always uses a discretized version of the original signal f that is obtained by sampling f on a discrete set. The question then arises whether and how f can be recovered from its samples, which often is termed as “sampling problem”. Therefore, among other things, one objective of research on the sampling problem is to quantify the conditions under which it is possible to recover particular classes of functions from different sets of discrete samples, and a related objective is to use these analytical results to develop explicit reconstruction schemes for the analysis and processing of digital data. Specifically, the sampling problem often includes two main parts as follows:

[0007] (a) First, given a class of functions V on d, find conditions on sampling sets X={xjεd:jεJ}, where J is a countable index set, under which a function fεV can be reconstructed uniquely and stably from its samples {f(xj):xjεX)}; and

[0008] (b) Second, find efficient and fast numerical algorithms that recover any function fεV from its samples on X.

[0009] In some applications, it may be justified to assume that the sampling set X={xj:jεJ} is uniform, i.e., that X forms a regular n-dimensional Cartesian grid, where n is an integer greater than zero. For examples, as shown in FIG. 1, curve 101 representing a function f is sampled on a one-dimensional uniform grid 103. Moreover, as shown in FIG. 2, grid 203 represents a two-dimensional Cartesian grid. These grids can find some uses in real world and can be used to practice the present invention. For example, a digital image is often acquired by sampling light intensities on a uniform grid. Data acquisition requirements and the ability to process and reconstruct the data simply and efficiently often justify this type of uniform data collection.

[0010] However, in many realistic situations the data are known only on a nonuniformly spaced sampling set. This nonuniformity is a fact of life and prevents the use of the standard methods from Fourier analysis. Correspondingly, several nonuniform grids, such as a polar sampling grid 213, a spiral sampling grid 223, and a nonuniform grid 233, all schematically shown in FIG. 2, can find some more uses in real world and can be used to practice the present invention. The following examples are typical and indicate that nonuniform sampling problems are more pervasive in science and engineering.

[0011] Communication theory and practice: When data from a uniformly sampled signal (function) are lost, the result is generally a sequence of nonuniform samples. This scenario is usually referred to as a missing data problem. Often, missing samples are due to the partial destruction of storage devices, e.g., scratches on a CD. As shown in FIG. 3, a missing data problem is simulated by randomly removing samples from a slice of a three-dimensional magnetic resonance (MR) digital image 301. The simulated image 303 with 50% random missing samples displays a rather incomplete construction of image 301 and loses important details of image 301. Image 303 can be considered as a sampling of image 301, which is the object of interest at this example, on the nonuniform grid 233.

[0012] Astronomical measurements: The measurement of star luminosity gives rise to extremely nonuniformly sampled time series. Daylight periods and adverse nighttime weather conditions prevent regular data collection as further discussed in [111].

[0013] Medical imaging: Computerized tomography (CT) and magnetic resonance imaging (MRI) frequently use the nonuniform polar sampling grid 213 and spiral sampling grid 223 as further discussed in [21, 90], respectively.

[0014] Other applications using nonuniform sampling sets occur in geophysics as further discussed in [92], spectroscopy as further discussed in [101], general signal/image processing as further discussed in [13, 22, 103, 106], and biomedical imaging such as MRI as discussed above in connection with images 301 and 303 and ultrasonic images 401 and 403 as shown in FIG. 4 and as further discussed in [20, 59, 90, 101], respectively. For examples, as shown in FIG. 4, detected edge points 402 of the left ventricle of a heart from a two-dimensional ultrasound image 401 constitute a nonuniform sampling of the left ventricle's contour, while boundary 404 of the left ventricle shown is reconstructed from the detected edge sample points 402. More information about modern techniques for nonuniform sampling and applications can be found in [16].

[0015] Efforts have been made accordingly to find better, fast and efficient techniques for nonuniform data sampling and data reconstruction such that one can reconstruct an object of interest from nonuniform data sampling with limit data points, and/or can convert from digital type (i.e., some data points) to analog type (i.e., a corresponding image), or vice versa.

[0016] However, among other things, standard methods of signal processing and data sampling encounter several problems as follows:

[0017] First, these methods often are not local due to the infinite range of the band limited function model utilized by these methods.

[0018] Second, these methods are not multi-dimensional, i.e., they do not work well for signals related to time series, images, 3-D images, spectral images as well as samples of any high dimensional images with dimension larger than 3.

[0019] Third, these methods are normally not that fast. In fact, standard methods based on Fast Fourier Transformation (“FFT”) are currently fastest methods one can get. Still, in an ideal situation, it has an order of a computational complexity NlogN that is slow when N is large, where N is the order of the samples.

[0020] Fourth, these methods may not be easily utilized in association with parallel processing because the non-localization due to the infinite range of the band limited function model utilized by these methods, which further limits the ability of these methods to take advantages of computing power offered by parallel processing.

[0021] Therefore, there is a need to develop system and methods for nonuniform sampling and applications of data representing an object to address the aforementioned deficiencies and inadequacies.

SUMMARY OF THE INVENTION

[0022] The present invention overcomes the disadvantages of the prior art and relates to new approaches to sampling, converting and processing data. Specifically, the present invention relates to system and methods of sampling and reconstruction of data from a data set that contains information regarding one or more objects of interest. Moreover, the present invention relates to system and methods of converting data between a digital form and an analog form for a data set that contains information regarding one or more objects of interest.

[0023] In one aspect, the present invention relates to a method for sampling data for an object of interest, which is characterized by a functions, from a data set {f(xj)}obtained from data sample points xjε[M0, M1], j being an indexing number, [M0, M1] representing an interval where the function f is sampled, comprising the steps of:

[0024] a. constructing a first matrix =M0M1, with matrix elements in the form of jk=φxj−k), wherein φ is a continuous function with compact support of size S, S being an integer, j being an integer from 1 to J, J being an integer ≧M1−M0+2S−1, and k being an integer from M0−S+1 to M1+S−1;

[0025] b. obtaining a second matrix =M0, M1=*, with matrix elements in the form of 1${\mathrm{kl}}_{}=\sum _{j=1}^J\stackrel{\_}{\phi \left(x_j-k\right)}\phi \left(x_j-l\right),$

[0026] wherein each of k, l is an integer from M0−S+1 to M1+S−1;

[0027] c. obtaining a first vector b=*y with vector elements in the form of 2$b_k=\sum _{j=1}^J\stackrel{\_}{\phi \left(x_j-k\right)}{y}_j$

[0028] for k=M0−S+1, . . . ,M1+S−1, wherein y=(y1, . . . , yj) is a sampling vector;

[0029] d. solving a plurality of equations in the form of c=lb to obtain a second vector c=(ck), wherein k is an integer from M0−S+1 to M1+S−1;

[0030] e. computing the restriction of the function f to [M0,M1] to obtain a new data set, f(x), in the form of 3$f\left(x\right)=\sum _{k=M_0-S+1}^{M_1+S-1}{c}_k\phi \left(x-k\right)\mathrm{for}x\in \left[M_0,M_1\right];\mathrm{and}$

[0031] f. reconstructing the object of interest f from the new data set.

[0032] In one embodiment, the step of reconstructing the object of interest f from the new data set comprises the step of obtaining the object of interest f as the least square approximation of the sampling vector y from 4$\sum _{j=1}^J{\&LeftBracketingBar;y_j-f\left(x_j\right)\&RightBracketingBar;}^2=\underset{h\in V\left(\phi \right)}{\min }\sum _{j=1}^J{\&LeftBracketingBar;y_j-h\left(x_j\right)\&RightBracketingBar;}^2,$

[0033] wherein V(φ) represents a shift-invariant space. When y arises as the sampling vector, y={f(xj)}, an exact reconstruction of f is obtained from the new data set.

[0034] In another aspect, the present invention relates to a system for sampling data for an object of interest, which is characterized by a function f, from a data set {f(xj)} obtained from data sample points xjε[M0, M1], j being an indexing number, [M0, M1] representing an interval where the function f is sampled, comprising a processor, wherein the processor is capable of performing the steps of:

[0035] a. constructing a first matrix =M0,M1 with matrix elements in the form of jk=φ(xj—k), wherein φ is a continuous function with compact support of size S, S being an integer, j being an integer from 1 to J, J being an integer ≧M1−M0+2S−1, and k being an integer from M0−S+1 to M1+S−1;

[0036] b. obtaining a second matrix =M0,M1=* with matrix elements in the form of 5${\mathrm{kl}}_{}=\sum _{j=1}^J\stackrel{\_}{\phi \left(x_j-k\right)}\phi \left(x_j-l\right),$

[0037] wherein each of k, l is an integer from M0−S+1 to M1+S−1;

[0038] c. obtaining a first vector b=*y with vector elements in the form of 6$b_k=\sum _{j=1}^J\stackrel{\_}{\phi \left(x_j-k\right)}{y}_j$

[0039] for k=M0−S+1, . . . ,M1+S−1, wherein y=(y1, . . . , yj) is a sampling vector;

[0040] d. solving a plurality of equations in the form of c=1b to obtain a second vector c=(ck), wherein k is an integer from M0−S+1 to M1+S−1;

[0041] e. computing the restriction of the function f to [M0,M1] to obtain a new data set, f(x), in the form of 7$f\left(x\right)=\sum _{k=M_0-S+1}^{M_1+S-1}{c}_k\phi \left(x-k\right)\mathrm{for}x\in \left[M_0,M_1\right];\mathrm{and}$

[0042] f. reconstructing the object of interest f from the new data set.

[0043] In yet another aspect, the present invention relates to a method for converting data for an object of interest, which is characterized by a functions between a digital form and an analog form, comprising the steps of:

[0044] a. locally selecting a plurality of data sample points in the order of N, N being an integer, wherein the data sample points are in a first form of data type;

[0045] b. performing a transformation in a shift-invariant space to the locally selected data sample points to obtain a new data set that is in a second form of data type different from the first form of data type; and

[0046] c. reconstructing the object of interest f from the new data set.

[0047] In one embodiment, the first form of data type is one of the digital form and the analog form, and the second form of data type is one of the digital form and the analog form that is different from the first form of data type. In other words, the first form of data type can be in digital form, and the corresponding second form of data type is in analog form, and vice versa. Alternatively, the present invention can be practiced to convert data from digital form to digital form, or from analog form to analog form.

[0048] Furthermore, the step of locally selecting a plurality of data sample points comprises the step of obtaining a data set {f(xj)} from data sample points xjε[M0, M1], j being an indexing number up to the order of N, [M0, M1] representing an interval where the function f is sampled.

[0049] Moreover, the step of performing a transformation in a shift-invariant space comprises the steps of:

[0050] a. constructing a first matrix =M0,M2 with matrix elements in the form of jk=φ(xj−k), wherein φ is a continuous function with compact support of size S, S being an integer, j being an integer from 1 to J, J being an integer ≧M1−M0+2S−1, and k being an integer from M0−S+1 to M1+S−1;

[0051] b. obtaining a second matrix =M0,M1=* with matrix elements in the form of 8${\mathrm{kl}}_{}=\sum _{j=1}^J\stackrel{\_}{\phi \left(x_j-k\right)}\phi \left(x_j-l\right),$

[0052] wherein each of k, l is an integer from M0−S+1 to M1+S−1;

[0053] c. obtaining a first vector b=* y with vector elements in the form of 9$b_k=\sum _{j=1}^J\stackrel{\_}{\phi \left(x_j-k\right)}{y}_j$

[0054] for k=M0−S+1, . . . ,M1+S−1, wherein y=(y1, . . . , yj) is a sampling vector;

[0055] d. solving a plurality of equations in the form of c=b to obtain a second vector c=(ck), wherein k is an integer from M0−S+1 to M1+S−1; and

[0056] e. computing the restriction of the function f to [M0, M1] to obtain the new data set, f(x), in the form of 10$f\left(x\right)=\sum _{k=M_0-S+1}^{M_1+S-1}{c}_k\phi \left(x-k\right)\mathrm{for}x\in \left[M_0,M_1\right].$

[0057] Additionally, the step of reconstructing the object of interest f from the new data set comprises the step of obtaining the object of interest f as the least square approximation of the sampling vector y from 11$\sum _{j=1}^J{\&LeftBracketingBar;y_j-f\left(x_j\right)\&RightBracketingBar;}^2=\underset{h\in V\left(\phi \right)}{\min }\sum _{j=1}^J{\&LeftBracketingBar;y_j-h\left(x_j\right)\&RightBracketingBar;}^2,$

[0058] wherein V(φ) represents the shift-invariant space.

[0059] In a further aspect, the present invention relates to a system for converting data for an object of interest, which is characterized by a functions between a digital form and an analog form, comprising a processor, wherein the processor is capable of performing the steps of:

[0060] a. locally selecting a plurality of data sample points in the order of N, N being an integer, wherein the data sample points are in a first form of data type;

[0061] b. performing a transformation in a shift-invariant space to the locally selected data sample points to obtain a new data set that is in a second type of data type different from the first form of data type; and

[0062] c. reconstructing the object of interest f from the new data set.

[0063] In one aspect, the present invention relates to a method for recovering information about an object of interest f from a data set {f(xj)} obtained from nonuniform data sample points xjεX, j being an indexing number, X representing the set of points where the object of interest f has been sampled, comprising the steps of:

[0064] a. selecting data f(xj) from the data set {f(xj)};

[0065] b. constructing an approximation operator QX;

[0066] c. applying the approximation operator QX to the data f(xj) to obtain an approximation QXf(xj);

[0067] d. constructing a projection operator P;

[0068] e. applying the projection operator P to the approximation QXf(xj) to obtain a first approximation f1=P QXf(xj);

[0069] f. obtaining an error e=f(xj)−f1;

[0070] g. applying the projection operator P to the error e to obtain a first approximation of error e1=P e;

[0071] h. obtaining a second approximation f2=f1+e1; and

[0072] i. returning to step (e) until a sequence fn=f1+e1+e2+e3+ . . . +en−1 is obtained, wherein function fn converges to the object of interest f.

[0073] Another aspect of the present invention relates to a system for recovering information about an object of interest f from a data set {f(xj)} obtained from nonuniform data sample points xjεX, j being an indexing number, X representing the set of points where the object of interest f has been sampled. In one embodiment, the system has an input device for receiving the information in the data set {f(xj)}, and a processing unit communicating to the input device and performing the steps of:

[0074] a. selecting data f(xj) from the data set {f(xj)};

[0075] b. constructing an approximation operator QX;

[0076] c. applying the approximation operator QX to the data f(xj) to obtain an approximation QXf(xj);

[0077] d. constructing a projection operator P;

[0078] e. applying the projection operator P to the approximation QXf(xj) to obtain a first approximation f1=P QXf(xj);

[0079] f. obtaining an error e=f(xj)−f1;

[0080] g. applying the projection operator P to the error e to obtain a first approximation of error e1=P e;

[0081] h. obtaining a second approximations f2=f1+e1; and

[0082] i. returning to step (e) until a sequence fn=f1+e1+e2+e3+ . . . +en−1 is obtained, wherein function fn converges to the object of interest f.

[0083] A further aspect of the present invention relates to a method for recovering information about an object of interest f from a data set {f(xj))} obtained from nonuniform data sample points xjεX, j being an indexing number, X representing the set of points where the object of interest f has been sampled in presence of noise. Practicing the present invention in one embodiment includes the steps of:

[0084] a. selecting data fj′ from the data set {fj′=f(xj)+ηj}, wherein ηj represents a corresponding noise component;

[0085] b. constructing an initialization functions as the summation of function fj′: f′=ΣjεJfj′βj, wherein βj represents the jth component of a partition of unity;

[0086] C. constructing an approximation operator QX;

[0087] d. applying the approximation operator QX to the function f′ to obtain an approximation QXf′;

[0088] e. constructing a projection operator P;

[0089] f. applying the projection operator P to the approximation QXf′ to obtain a first approximation f1=P QXf′;

[0090] g. obtaining an error e=f−f1;

[0091] h. applying the projection operator P to the error e to obtain a first approximation of error e1=P e;

[0092] i. obtaining a second approximation f2=f1+e1; and

[0093] j. returning to step (f) until a sequence fn=f1+e1+e2+e3+ . . . +en−1 is obtained, wherein function fn converges to a function f∞that adequately describes the object of interests

[0094] Another further aspect of the present invention relates to system for recovering information about an object of interest f from a data set {f(xj)} obtained from nonuniform data sample points xjε:X, j being an indexing number, X representing the set of points where the object of interest f has been sampled in presence of noise. In one embodiment, the system has an input device for receiving the information in the data set {f(xj)}, and a processing unit coupled to and communicating to the input device and performing the steps of:

[0095] a. selecting data fj′ from the data set {fj′=f(xj)+σj}, wherein ηj represents a corresponding noise component;

[0096] b. constructing an initialization function f′ as the summation of function fj′: f′=ΣjεJfj′βj, wherein βj represents the jth component of a partition of unity;

[0097] c. constructing an approximation operator QX;

[0098] d. applying the approximation operator QX to the function f′ to obtain an approximation QXf′;

[0099] e. constructing a projection operator P;

[0100] f. applying the projection operator P to the approximation QXf′ to obtain a first approximation f1=P QXf′;

[0101] g. obtaining an error e=f−f1;

[0102] h. applying the projection operator P to the error e to obtain a first approximation of error e1=P e;

[0103] i. obtaining a second approximation f2=f1+e1; and

[0104] j. returning to step (f) until a sequence fn=f1+e1+e2+e3+ . . . +en−1 is obtained, wherein function fn converges to a function f∞ that adequately describes the object of interest f.

[0105] In yet another aspect, the present invention relates to a method for recovering information about an object of interest f from a data set {f(xj)} obtained from nonuniform data sample points xjεX, j being an indexing number, X representing the set of points where the object of interest f has been sampled in presence of noise, comprising the steps of sampling a set of data samples from the data set {f(xj)}, iteratively framing the sampled set of data samples until the sampled set of data samples converges to a new set of data samples that has a size significantly smaller than the data set {f(xj)}, and reconstructing the object of interest f from the new set of data samples.

[0106] In another aspect, the present invention relates to a system for recovering information about an object of interest f from a data set {f(xj)}obtained from nonuniform data sample points xjεX, j being an indexing number, X representing the set of points where the object of interest f has been sampled in presence of noise, comprising means for sampling a set of data samples from the data set {f(xj)}, means for iteratively framing the sampled set of data samples until the sampled set of data samples converges to a new set of data samples that has a size significantly smaller than the data set {f(xj)}, means for reconstructing the object of interest f from the new set of data samples. In one embodiment, the sampling means comprises an input device, the framing means comprises a processor, and the reconstructing means comprises an output device.

[0107] In yet another aspect, the present invention relates to a method for recovering information about an object of interest f from a data set of averages {(f, ψxj} obtained from nonuniform data sample points xjεX, j being an indexing number, X representing the set of points where the object of interest f has been sampled, comprising the steps of:

[0108] a. selecting data f, ψxj from the data set {f, ψxj};

[0109] b. constructing a quasi-reconstruction operator Ax for obtaining a quasi-reconstruction AXf, wherein AXf=Σjεj f, ψxjβj, wherein βj represents the jth component of a partition of unity, and {ψxj:xjεX} is a set of functionals that act on f;

[0110] c. constructing a projection operator P;

[0111] d. applying the projection operator P to the quasi-reconstruction Axf to obtain a first approximation f1=PAXf;

[0112] e. obtaining an error e=(f−f1);

[0113] f. applying the projection operator P and the quasi-reconstruction operator AX to the error e to obtain a first approximation of error e1=P AXe;

[0114] g. obtaining a second approximation f2=f1+e1; and

[0115] h. returning to step (e) until a sequence fn=f1+e1+e2+e3+ . . . +en−1=P AX(f−fn−1)+fn−1 is obtained, wherein function fn converges to the object of interest f.

[0116] Another aspect of the present invention relates to a system for recovering information about an object of interest f from a data set of averages {f,ψxj} obtained from nonuniform data sample points xjεX, j being an indexing number, X representing the set of points where the object of interest f has been sampled. In one embodiment, the system has an input device for receiving the information in the data set {f,ψxj}, and a processor coupled to and communicating to the input device and performing the steps of:

[0117] a. selecting data f,ψxj from the data set {f,ψxj};

[0118] b. constructing a quasi-reconstruction operator AX for obtaining a quasi-reconstruction AXf, wherein AXf=ΣjεJf,ψxjβj, wherein βj represents the jth component of a partition of unity, and {ψxj:xjεX} is a set of functionals that act on f,

[0119] c. constructing a projection operator P;

[0120] d. applying the projection operator P to the quasi-reconstruction AXf to obtain a first approximation f1=PAXf;

[0121] e. obtaining an error e=f−f1;

[0122] f. applying the projection operator P and the quasi-reconstruction operator AX to the error e to obtain a first approximation of error e1=P AX e;

[0123] g. obtaining a second approximation f2=f1+e1; and

[0124] h. returning to step (e) until a sequence fn=f1+e1+e2+e3+ . . . +en−1=P AX(f−fn−1)+fn−1 is obtained, wherein function fn converges to the object of interest f.

[0125] A further aspect of the present invention relates to a method for recovering information about an object of interest f from a data set {f,ψxj} obtained from nonuniform data sample points xjεX, j being an indexing number, X representing the set of points where the object of interest f has been sampled in presence of noise. Practicing the present invention in one embodiment includes the steps of:

[0126] a. selecting data fj′ from the data set {fj′=f(xj)+ηj}, wherein ηj represents a corresponding noise component;

[0127] b. constructing an initialization function f′ as the summation of function fj′: f′=ΣjεJfj′βj, wherein βj represents the jth component of a partition of unity;

[0128] c. constructing a quasi-interpolant operator QX;

[0129] d. applying the quasi-interpolant operator QX to the function f′ to obtain an approximation QXf′;

[0130] e. constructing a projection operator P;

[0131] f. applying the projection operator P to the approximation QXf′ to obtain a first approximation f1=P QXf′;

[0132] g. obtaining an error e=f−f1;

[0133] h. constructing a quasi-reconstruction operator AX for obtaining a quasi-reconstruction AXf, wherein AXf=ΣjεJf,ψxjβj, and {ψxj:xjεX} is a set of functionals that act on f;

[0134] i. applying the projection operator P and the quasi-reconstruction operator AX to the error e to obtain a first approximation of error e1=PAX e;

[0135] j. obtaining a second approximation f2=f1+e1; and

[0136] k. returning to step (g) until a sequence fn=f1+e1+e2+e3+ . . . +en−1=f1+(I−P AX)fn−1 is obtained, wherein I is an unit operator and function fn converges to a function f∞ that adequately describes the object of interest f.

[0137] A further aspect of the present invention relates to system for recovering information about an object of interest f from a data set {f,ψxj} obtained from nonuniform data sample points xjεX, j being an indexing number, X representing the set of points where the object of interest f has been sampled in presence of noise. In one embodiment, the system has an input device for receiving the information in the data set {f(xj)}, and a processor coupled to and communicating to the input device and performing the steps of:

[0138] a. selecting data fj′ from the data set {fj′=f,ψxj+ηj}, wherein ηj represents a corresponding noise component;

[0139] b. constructing an initialization function f′ as the summation of function fj′: f′=ΣjεJfj′βj, wherein βj represents the jth component of a partition of unity;

[0140] c. constructing a quasi-interpolant operator QX;

[0141] d. applying the quasi-interpolant operator QX to the functions to obtain an approximation QXf′;

[0142] e. constructing a projection operator P;

[0143] f. applying the projection operator P to the approximation QXf′ to obtain a first approximation f1=P QXf′;

[0144] g. obtaining an error e=f−f1;

[0145] h. constructing a quasi-reconstruction operator AX for obtaining a quasi-reconstruction AXf wherein AXf=ΣjεJ f,ψxjβj, and {ψxj:xjεX} is a set of functionals that act on f;

[0146] i. applying the projection operator P and the quasi-reconstruction operator AX to the error e to obtain a first approximation of error e1=PAX e;

[0147] j. obtaining a second approximation f2=f1+e1; and

[0148] k. returning to step (g) until a sequence fn=f1+e1+e2+e3+ . . . +en−1=f1+(I−P AX)fn−1 is obtained, wherein I is an unit operator and function fn converges to a function f∞ that adequately describes the object of interest f.

[0149] In one embodiment, the approximation operator QX comprises an interpolation operator, wherein the interpolation operator provides piecewise linear interpolation. Alternatively, the approximation operator QX comprises a quasi-interpolation operator, wherein the quasi-interpolation operator provides step-wise approximation. The projection operator P comprises a bounded projection. Alternatively, the projection operator P comprises a universal projection.

[0150] The data set, being in the form of {f(xj)}, {fj′=f(xj)+ηj}, {f,ψxj}, or {fj′=f,ψxj+ηj}, can comprise a variety of types of data in a first form of analog data type or in a second form of digital data type, or vice versa. Each type of data can comprise a variety of kinds of data in terms of substance. For examples, the data set can comprise a plurality of digital data related to sound including communications between person to person, person to machine, machine to machine including, but not limited to, data collected for commercial purpose, educational purpose, intelligence purpose, national security purpose, etc. The data set can also comprise a plurality of digital data related to at least one image, wherein the at least one image is one of a plurality of magnetic resonance image, a plurality of computerized tomography images, a plurality of optical image, a plurality of ultra sound images, a plurality of electronic images that are transmittable over a network, a plurality of satellite images, a plurality of three-dimensional images, a plurality of spectral images, a plurality of n-dimensional images with n being an integer greater than 1, any combination of them, or the like. Moreover, the data set can comprise a plurality of biological data, wherein the plurality of biological data contain a plurality of data from the group of a plurality of biological data related to genes, a plurality of biological data related to proteins, a plurality of biological data related to cells, a plurality of biological data related to bacteria, or a plurality of biological data related to tissues, any combination of them, or the like. Furthermore, the data set can comprise a plurality of data from at least one measurement of the object of the interest, wherein the at least one measurement is one of at least one astronomical measurement, at least one seismic measurement, at least one marine measurement, at least one geophysical measurement, at least one atmospheric measurement, at least one engineering measurement, at least one physical measurement, or at least one chemical measurement, any combination of them, or the like. Additionally, the data set can comprise a plurality of data that are static, a plurality of data related to time dependent signals, any combination of them, or the like.

[0151] In one embodiment, the processor comprises a microprocessor. The system further has an input device coupled to and in communication with the processor, wherein the input device comprises at least one device selected from the group of a processor interface, a GUI, a scanner, a CD-ROM, a diskette, a computer coupled to a network, and a networking device. Moreover, the system further has an output device coupled to and in communication with the processor, wherein the output device comprises at least one device selected from the group of a GUI, a printer, a CD-ROM, a diskette, a memory device, a computer coupled to a network, and a networking device.

[0152] For the purpose of practicing the present invention, an object of interest can be any subject under the sun in which human being is interested, which includes human being himself or herself and any and all human being activities.

[0153] Other advantages and uses for the present invention will be more clearly understood by reference to the remainder of this document.

[0154] The present invention overcomes the disadvantages of the prior art and relates to new approaches to sampling, converting and processing data. Specifically, the present invention relates to system and methods of sampling and reconstruction of data from a data set that contains information regarding one or more objects of interest. Moreover, the present invention relates to system and methods of converting data between a digital form and an analog form for a data set that contains information regarding one or more objects of interest.

[0155] Many more applications can be found for these and other aspects of the present invention. For examples, among other things and advantages, the methods and related algorithms of the present invention are local, i.e., the reconstruction of the data or signal of an object of interest depends only on the given data in the neighborhood of each sampling point as opposed to the infinite range of the band limited function model of the prior art.

[0156] Second, the methods and related algorithms of the present invention are multi-dimensional, which make them suitable for signals, time series, images, 3-D images, spectral images as well as samples of any high dimensional images with dimension larger than n, n being an integer greater than 1.

[0157] Third, the methods and related algorithms of the present invention are fast in the order of N with a computational complexity of order N, where N is about the number of samples. In other words, the number of operation needed for reconstruction is about and proportional to the number of samples. By contrast, as set forth above, standard methods based on FFT are of complexity NlogN. Consequently, the methods and related algorithms of the present invention are well adapted for processing large data sets in the fields of, for examples, medical image processing, geophysics, seismology, and oil exploration.

[0158] Fourth, localization of the methods and related algorithms of the present invention makes it as a great candidate in combination with data segmentation to facilitate parallel processing, thereby increasing speed of reconstruction.

[0159] Fifth, the methods and related algorithms of the present invention can be used to improve data conversion between digital form and analog form in many applications such as improving recovery of faults in a fault tolerant system rather than building heavily redundant systems, improving analysis/correction of outlier data points for reduction of noise, application in medical systems including computer x-ray tomography, magnetic resonance imaging, ultra-sound images, telemedicine such as remote specialist interpretation of medical images and the like, improving registration/synthesis of color separated images in printing processes and systems, enhancement of image or sound files recorded with older or inferior technologies including movies, early TV, early photographs, or the like, improving satellite imaging, applications in geophysics, seismology, oil exploration and the like.

[0160] Indeed, the methods and related algorithms of the present invention can be utilized to provide better, efficient, and fast data processing, data conversion, and data reconstruction for any and all fields where data for an object of interest are collected and need to be processed.

[0161] These and other aspects of the present invention will become apparent from the following description of the preferred embodiment taken in conjunction with the following drawings, although variations and modifications therein may be affected without departing from the spirit and scope of the novel concepts of the disclosure.

BRIEF DESCRIPTION OF THE DRAWINGS

[0162] FIG. 1 schematically shows a typical sampling problem that can be solved by practicing the present invention. (a). A function f defined on , represented by curve 101, is sampled on a uniform grid 103. (b). The same function f is sampled on a nonuniformly spaced grid 113. The sampling locations xj are marked by the symbol x, and the sampled values f(xj) by a circle o.

[0163] FIG. 2 schematically shows several types of sampling grids that can be utilized to practice the present invention. Top left: A Cartesian sampling grid 203, which is used in signal and image processing whenever possible because its simplicity. Top right: A polar sampling grid 213, which is used, for example, in computerized tomography. In this case, the two-dimensional Fourier transform {circumflex over (f)} (not shown) is sampled with the goal of reconstructing f. Bottom left: A spiral sampling grid 223, which is used, for example, for fast MRI by direct signal reconstruction from spectral data on spirals. Bottom right: A typical nonuniform sampling grid 233, which can be used, for examples, in spectroscopy, astronomy, geophysics, and other signal and image processing applications.

[0164] FIG. 3 schematically shows a missing data problem that can be solved by practicing the present invention. Left: Original digital MRI image 301 with 128×128 samples. Right: MRI image 303 with 50% randomly missing samples.

[0165] FIG. 4 schematically shows sampling and boundary reconstruction from ultrasonic images that can be obtained by practicing the present invention. Left: Detected edge points 402 of the left ventricle of a heart from a two-dimensional ultrasound image 401 constituting a nonuniform sampling of the left ventricle's contour. Right: Boundary 404 of the left ventricle reconstructed from the detected edge sample points 402.

[0166] FIG. 5 schematically shows the iterative reconstruction algorithm according to one embodiment of the present invention.

[0167] FIG. 6 schematically shows reconstruction of a function f with ∥f∥2=3.5 using the iteration algorithm as set forth in FIG. 5 according to one embodiment of the present invention. (a). Curve 601 representing function f belonging to the shift-invariant space generated by the Gaussian function exp(−x2/σ2), σ≈0.81, and its sample values {f(xj), xjεX} marked by + (density ε≈0.8). (b). Curve 611 representing error ∥f−fn|L2 against the number of iterations. (c). Curve 621 representing final error f−fn after 10 iterations. (d). Curve 631 representing reconstructed function f10 (continuous line) and original samples {f(xj), xjεX}.

[0168] FIG. 7 displays a missing data reconstruction performed according to one embodiment of the present invention. (a). Original digital MRI image 701 with 128×128 samples. (b). MRI image 703 with 50% randomly missing samples. (c). Reconstruction or image 711 using the iterative reconstruction algorithm according to one embodiment of the present invention. The corresponding shift invariant space is generated by φ(x, y)=β3(x)×β3(y), where β3=χ[0, 1]*χ[0, 1]*χ[0, 1]*, χ[0, 1] is the B-spline function of degree 3.

[0169] FIG. 8 schematically shows reconstruction of a function f with additive noise using the iteration algorithm as set forth in FIG. 5 according to one embodiment of the present invention. (a). Curve 801 representing function f belonging to the shift-invariant space generated by the Gaussian function exp(−x2/σ2), σ=0.81. (b). Curve 811 representing function f with an additive white noise (SNR≈0 db). (c). Curve 821 representing Noisy signal sampled on a nonuniform grid with maximal gap 0.51. (d). Curve 831 representing reconstructed function f10 after 10 iterations (continuous line) and original signal f(dotted line) from curve 801.

[0170] FIG. 9 schematically shows a block-diagram of a system that can be utilized to practice the present invention.

[0171] FIG. 10 shows the shift-invariant spline spaces with the B-spline of order 3 with φ[−2,2] and S=2 that can be utilized to practice the present invention.

[0172] FIG. 11 shows a curve 1101 as a plot of the operation count as a function of the number of sampling points.

[0173] FIG. 12 schematically shows reconstruction of a function f with additive noise according to one embodiment of the present invention with a signal on the interval [0,128]. Since S=2, one needs at least M1−M0+2S−1=131 samples. The actual sampling set of FIG. 12 has 200 sample points as calculated and satisfies the maximum gap condition maxj(xj+1−xj)≈0.67<1. (a). Curve 1201 representing original signal. (b). Curve 1211 representing original signal with noise. (c). Curve 1221 representing noisy signal sampled on a nonuniform grid. (d). Curve 1231 representing reconstructed original signal after 10 iterations (continuous line) in comparison with original signal dotted line from curve 1201.

[0174] FIG. 13A schematically shows the dependence of the reconstruction on the generator φ, where the generator φ is a B-spline of order N, i.e. φN=X[−1/2,1/2]* . . . *[1/2,1/2] (N+1-fold convolution). The data set (xj, yj)j=1, . . . ,J is generated by sampling a function fεV (φ5) and then noise is added. The top left picture shows the original signal and the noisy sampled data. Then each subplot depicts the optimal approximation of these data in the spline space V (φN), N=0, . . . , 6, starting with an approximation by a step function frecεV (φ0) via an approximation by a piecewise linear function frecεV(φ1) and ending with a smooth approximation frecεV (φ6). (a). Curve 1301 representing the original signal and curve 1303 representing the corresponding noisy sampled data. (b). Curve 1311 representing the optimal approximation of the data as given in (a) in the spline space V(φN), N=0, starting with an approximation by a step function frecεV (φ0), i.e., N=0, in comparison with curve 1301 representing original signal with noise. (c). Same as (b) but with N=1. (d). Same as (b) but with N=2.

[0175] FIG. 13B schematically shows the dependence of the reconstruction on the generator φ, where the generator φ is a B-spline of order N, i.e. φN=X[−1/2,1/2]* . . . *[1/2,1/2] (N+1-fold convolution) same as FIG. 13A. (e). Same as (b) of FIG. 13A but with N=3. (f). Same as (e) but with N=4. (e). Same as (b) but with N=5. (e). Same as (e) but with N=6.

[0176] FIG. 14 schematically shows the associated banded matrix of the linear system according to one embodiment of the present invention.

[0177] FIG. 15 schematically shows the power of the method of data segmentation according to one embodiment of the present invention.

DETAILED DESCRIPTION OF THE INVENTION

[0178] Various embodiments of the invention are now described in detail. Referring to the drawings, like numbers indicate like parts throughout the views. As used in the description herein and throughout the claims that follow, the meaning of “a,” “an,” and “the” includes plural reference unless the context clearly dictates otherwise. Also, as used in the description herein and throughout the claims that follow, the meaning of “in” includes “in” and “on” unless the context clearly dictates otherwise.

[0179] The terms used in this specification generally have their ordinary meanings in the art, within the context of the invention, and in the specific context where each term is used. Certain terms that are used to describe the invention are discussed below, or elsewhere in the specification, to provide additional guidance to the practitioner in describing the compositions and methods of the invention and how to make and use them. For convenience, certain terms may be highlighted, for example using italics and/or quotation marks. The use of highlighting has no influence on the scope and meaning of a term; the scope and meaning of a term is the same, in the same context, whether or not it is highlighted. It will be appreciated that the same thing can be said in more than one way. Consequently, alternative language and synonyms may be used for any one or more of the terms discussed herein, nor is any special significance to be placed upon whether or not a term is elaborated or discussed herein. Synonyms for certain terms are provided. A recital of one or more synonyms does not exclude the use of other synonyms. The use of examples anywhere in this specification, including examples of any terms discussed herein, is illustrative only, and in no way limits the scope and meaning of the invention or of any exemplified term. Likewise, the invention is not limited to various embodiments given in this specification.

[0180] As used herein, “about” or “approximately” shall generally mean within 20 percent, preferably within 10 percent, and more preferably within 5 percent of a given value or range. Numerical quantities given herein are approximate, meaning that the term “about” or “approximately” can be inferred if not expressly stated.

[0181] Several embodiments are now described with reference to the FIGS. 1-15, in which like numbers indicate like parts throughout the FIGS. 1-15.

Overview of the Inventive System

[0182] In general, the present invention can be practiced in association with an inventive system 100 as shown in FIG. 9. In one embodiment, the system 100 includes a computer 8 having a host processor 12 with a display 14, such as a monitor, having a graphic-user interface (GUI) 20 displaying data for at least one peripheral device 10. The host processor 12 includes at least one microprocessor. The peripheral device 10 can be an input device to receive data or out put device for output data, shown here as a printer, in operative communication with the host processor 12. The printer 10 and host processor 12 can be in communication through any media, such as a direct wire connection 18 or through a network or the Internet 16. Even wireless communication between the printer 10 and host processor 12 can be used with the present invention as long as data can be exchanged between the printer 10 and host processor 12. The GUI 20 is generated by a GUI code as part of the operating system (O/S) of the host processor 12. As those skilled in the art will appreciate, the data to be processed in the system 100 can be stored in a memory device (not shown), can be imported from a device (not shown) in communication with the system 100, or can be provided in the form of a computer readable medium. Examples of computer readable media include: recordable type media such as floppy disks and CD-ROMs and transmission type media such as digital and analog communication links. Likewise, the result of processed data can be displayed in GUI 20, printed by printer 10, saved in a memory device (not shown), processed in another device, transmitted over network 16 or can be, provided in the form of the computer readable medium. Furthermore, the system 100 can be operated to practice the present invention with other data processing software such as Matlab© and Labview©.

[0183] Computer 8 can be a mainframe computer, a workstation, a desktop computer, a portable computer, any combination of them, or the like. Indeed, several computers 8 can be coupled together to form a local network for parallel processing. Or, computers 8 can communicate to each other through network 16 such as Internet for processing data according to one embodiment of the present invention.

[0184] The processing unit of each computer may comprise a processor such as host processor 12. The input device may include a microprocessor interface. The input device further may include at least one device selected from the group of a GUI, a scanner, a CD-ROM, a diskette, a computer coupled to a network, a networking device, any combination of them, or the like. Additionally, the system may further include an output device coupled to the processing unit, wherein the output device comprises at least one device selected from the group of a GUI, a printer, a CD-ROM, a diskette, a computer coupled to a network, a networking device, any combination of them, or the like.

[0185] Without intent to limit the scope of the invention, exemplary methods and their related results according to the embodiments of the present invention are given below. Note that titles or subtitles may be used in the examples for convenience of a reader, which in no way should limit the scope of the invention. Moreover, certain theories are proposed and disclosed herein; however, in no way they, whether they are right or wrong, should limit the scope of the invention so long as data are processed, sampled, converted, or the like according to the invention without regard for any particular theory or scheme of action.

EXAMPLES

Example 1

A Unified Framework for Uniform and Nonuniform Data Sampling in Shift Invariant Spaces

[0186] Approach Overview. Modern digital data processing of functions (or signals or images) always uses a discretized version of the original signal f that is obtained by sampling f on a discrete set. The question then arises whether and how f can be recovered from its samples.

[0187] In one aspect, the present invention relates to a method for recovering information about an object of interest f from a data set {f(xj)} obtained from uniform or nonuniform data sample points xjεX, j being an indexing number, X representing the set of points where the object of interest f has been sampled, wherein X may be perturbed by a noise, comprising the steps of sampling a set of data samples from the data set {f(xj)}, iteratively framing the sampled set of data samples until the sampled set of data samples converges to a new set of data samples that has a size significantly smaller than the data set {f(xj)}, and reconstructing the object of interest f from the new set of data samples.

1. Relevant Concepts

[0188] Sampling in Paley-Wiener Spaces: Bandlimited Functions. Since infinitely many functions can have the same sampled values on X={xj}jεJ⊂d, the sampling problem becomes meaningful only after imposing some a priori conditions on f. The standard assumption is that the function f on d belongs to the space of bandlimited functions BΩ; i.e., the Fourier transform 12$\hat{f}\left(\xi \right)={\int }_{R^d}f\left(x\right)e^{-2\pi i⟨\xi ,x⟩}dx$

[0189] of f is such that {circumflex over (f)} (ξ)=0 for all ξ∉Ω=[−ω, ω]d for some ω<∞ (see, e.g., [15, 44, 47, 55, 62, 72, 78, 88, 51, 112] and the review papers [27, 61, 65]). The reason for this assumption is a classical result of Whittaker [114] in complex analysis which states that, for dimension d=1, a function fεL2 ()∩B[−1/2,1/2] can be recovered exactly from its samples {f (k):kεZ} by the interpolation formula

f (x)=Σkεf(k)sinc(x−k)  (1.1)

[0190] where sinc (x)=sin πx/πx. This series gave rise to the uniform sampling theory of Shannon [96], which is fundamental in engineering and digital signal processing because it gives a framework for converting analog signals into sequences of numbers. These sequences can then be processed digitally and converted back to analog signals via (1.1).

[0191] Taking the Fourier transform of (1.1) and using the fact that the Fourier transform of the sinc function is the characteristic function χ[−1/2,1/2] shows that for any ξε[−½, ½] 13$\hat{f}\left(\xi \right)=\sum _k^{}f\left(k\right)e^{2\pi ik\xi }=\sum _k^{}{⟨\hat{f},e^{\mathrm{i2}\pi k}⟩}_{L^2\left(-1/2,1/2\right)}{e}^{i2\pi k\xi }.$

[0192] Thus, reconstruction by means of the formula (1.1) is equivalent to the fact that the set {ei2πkξ, kεZ} forms an orthonormal basis of L2 (−½, ½) called the harmonic Fourier basis. This equivalence between the harmonic Fourier basis and the reconstruction of a uniformly sampled bandlimited function has been extended to treat some special cases of nonuniformly sampled data. In particular, the results by Paley and Wiener [87], Kadec [71], and others on the nonharmonic Fourier bases {ei2πxkξ, kεZ} can be translated into results about nonuniform sampling and reconstruction of bandlimited functions [15, 62, 89, 94]. For example, Kadec's theorem [71] states that if X={xε:|xk−k|≦L<¼} for all kεZ, then the set {ei2πxkξ, kεZ} is a Riesz basis of L2 (−½, ½); i.e., {ei2πxkξ, kεZ} is the image of an orthonormal basis of L2 (−½, ½) under a bounded and invertible operator from L2 (−½, ½) onto L2 (−½, ½). Using Fourier transform methods, this result implies that any bandlimited function fεL2∩B[−1/2, 1/2] can be completely recovered from its samples {f(xk), kεZ}, as long as the sampling set is of the form X {xkε:|xk−k|<¼}kεZ.

[0193] The sampling set X={xkε:|xk−k|<¼}kεZ in Kadec's theorem is just a perturbation of Z. For more general sampling sets, the work of Beurling [23, 24], Landau [74], and others [18, 58] provides a deep understanding of the one-dimensional theory of nonuniform sampling of bandlimited functions. Specifically, for the exact and stable reconstruction of a bandlimited function f from its samples {f (xj):xjεX}, it is sufficient that the Beurling density 14$\begin{array}{cc}D\left(X\right)=\underset{r\to \infty }{\lim }\underset{y\in R}{\inf }\frac{\#X\bigcap \left(y+\left[0,r\right]\right)}{r}& \left(1.2\right)\end{array}$

[0194] satisfies D (X)>1. Conversely, if f is uniquely and stably determined by its samples on X⊂, then D (X)≧1 [74]. The marginal case D (X)=1 is very complicated and is treated in [79, 89, 94].

[0195] It should be emphasized that these results deal with stable reconstructions. This means that an inequality of the form 15${\&LeftDoubleBracketingBar;f\&RightDoubleBracketingBar;}_p\le {C\left(\sum _{x_j\in X}^{}{\&LeftBracketingBar;f\left(x_j\right)\&RightBracketingBar;}^p\right)}^{1/p}$

[0196] holds for all bandlimited functions fεLP∩BΩ. A sampling set for which the reconstruction is stable in this sense is called a (stable) set of sampling. This terminology is used to contrast a set of sampling with the weaker notion of a set of uniqueness. X is a set of uniqueness for BΩif f|x=0 implies that f=0. Whereas a set of sampling for B[1/2, 1/2] has a density D≧1, there are sets of uniqueness with arbitrarily small density. See [73, 25] for examples and characterizations of sets of uniqueness.

[0197] While the theorems of Paley and Wiener and Kadec about Riesz bases consisting of complex exponentials ei2πxkξ are equivalent to statements about sampling sets that are perturbations of Z, the results about arbitrary sets of sampling are connected to the more general notion of frames introduced by Duffin and Schaeffer [40]. The concept of frames generalizes the notion of orthogonal bases and Riesz bases in Hilbert spaces and of unconditional bases in some Banach spaces [2, 5, 6, 12, 14, 15, 20, 28, 29, 46, 66, 97].

[0198] Sampling in Shift-Invariant Spaces. The series (1.1) shows that the space of bandlimited functions B[1/2, 1/2] is identical with the space 16$\begin{array}{cc}{V}^2\left(\sin c\right)=\left\{\sum _{k\in Z}^{}{c}_k\sin c\left(x-k\right):\left(c_k\right)\in l^2\right\}.& \left(1.3\right)\end{array}$

[0199] Since the sinc function has infinite support and slow decay, the space of bandlimited functions is often unsuitable for numerical implementations. For instance, the pointwise evaluation 17$f\mapsto f\left(x_0\right)=\sum _{k\in Z}^{}{c}_k\sin c\left(x_0-k\right)$

[0200] is a nonlocal operation, because, as a consequence of the long-range behavior of sinc, many coefficients ck will contribute to the value f(x0). In fact, all bandlimited functions have infinite support since they are analytic. Moreover, functions that are measured in applications tend to have frequency components that decay for high frequencies, but these functions are not bandlimited in the strict sense. Thus, it has been advantageous to use non-bandlimited models that retain some of the simplicity and structure of bandlimited models but are more amenable to numerical implementation and are more flexible for approximating real data [13, 63, 64, 86, 103, 104]. One such example are the shift-invariant spaces that are discussed in this specification.

[0201] A shift-invariant space is a space of functions on d of the form 18$V\left({\phi }_1,\dots ,{\phi }_r\right)=\left\{\sum _{i=1}^r\sum _{j\in Z^d}^{}{c}_j^i{\phi }_i\left(x-j\right)\right\}.$

[0202] Such spaces have been used in finite elements and approximation theory [34, 35, 67, 68, 69, 98] and for the construction of multiresolution approximations and wavelets [32, 33, 39, 53, 60, 70, 82, 83, 95, 98, 99, 100]. They have been extensively studied in recent years (see, for instance, [6, 19, 52, 67, 68, 69]).

[0203] Sampling in shift-invariant spaces that are not bandlimited is a suitable and realistic model for many applications, e.g., for taking into account real acquisition and reconstruction devices, for modeling signals with smoother spectrum than is the case with bandlimited functions, or for numerical implementation [9, 13, 22, 26, 85, 86, 103, 104, 107, 110, 115, 116]. These requirements can often be met by choosing “appropriate” functions φi. This may mean that the functions φi have a shape corresponding to a particular “impulse response” of a device, or that they are compactly supported, or that they have a Fourier transform |{circumflex over (φ)}i(ξ)| that decays smoothly to zero as |ξ|→∞.

[0204] Uniform Sampling in Shift-Invariant Spaces. Early results on sampling in shift-invariant spaces concentrated on the problem of uniform sampling [7, 9, 10, 11, 37,64, 105, 108, 113, 116] or interlaced uniform sampling [110]. The problem of uniform sampling in shift-invariant spaces shares some similarities with Shannon's sampling theorem in that it requires only the Poisson summation formula and a few facts about Riesz bases [7, 9]. The connection between interpolation in spline spaces, filtering of signals, and Shannon's sampling theory was established in [11, 109]. These results imply that Shannon's sampling theory can be viewed as a limiting case of polynomial spline interpolation when the order of the spline tends to infinity [11, 109]. Furthermore, Shannon's sampling theory is a special case of interpolation in shift-invariant spaces [7, 9, 113, 116] and a limiting case for the interpolation in certain families of shift-invariant spaces V(φn) that are obtained by a generator φn=φ* . . . * φ consisting of the n-fold convolution of a single generator φ[9].

[0205] In applications, signals do not in general belong to a prescribed shift-invariant space. Thus, when using the bandlimited theory, the common practice in engineering is to force the function f to become bandlimited before sampling. Mathematically, this corresponds to multiplication of the Fourier transform {circumflex over (f)} of f by a characteristic function χΩ. The new function fa with Fourier transform {circumflex over (f)}={circumflex over (f)}χΩ is then sampled and stored digitally for later processing or reconstruction. The multiplication by χΩ before sampling is called prefiltering with an ideal filter and is used to reduce the errors in reconstructions called aliasing errors. It has been shown that the three steps of the traditional uniform sampling procedure, namely prefiltering, sampling, and postfiltering for reconstruction, are equivalent to finding the best L2-approximation of a function in L2∩BΩ[9, 105]. This procedure generalizes to sampling in general shift-invariant spaces [7, 9, 10, 85, 105, 108]. In fact, the reconstruction from the samples of a function should be considered as an approximation in the shift-invariant space generated by the impulse response of the sampling device. This allows a reconstruction that optimally fits the available samples and can be done using fast algorithms [106, 107].

[0206] Nonuniform Sampling in Shift-Invariant Spaces. The problem of nonuniform sampling in general shift-invariant spaces is more recent [4, 5, 30, 66, 75, 76, 77, 102, 119]. The earliest results [31, 77] concentrate on perturbation of regular sampling in shaft-invariant spaces and are therefore similar in spirit to Kadec's result for bandlimited functions. For the L2 case in dimension d=1, and under some restrictions on the shift-invariant spaces, several theorems on nonuniform sampling can be found in [76, 102]. Moreover, a lower bound on the maximal distance between two sampling points needed for reconstructing a function from its samples was given for the case of polynomial splines and other special cases of shift-invariant spaces in [76]. For the general multivariate case in Lp, the theory was developed in [4], and for the case of polynomial spline shift-invariant spaces, the maximal allowable gap between samples was obtained in [5]. For general shift-invariant spaces, a Beurling density D≧1 is necessary for stable reconstruction [5]. As in the case of bandlimited functions, the theory of frames is central in nonuniform sampling of shift-invariant spaces, and there is an equivalence between a certain type of frame and the problem of sampling in shift-invariant spaces [5, 66, 75].

[0207] Among other things, one aim of the present invention is to provide a unified framework for uniform and nonuniform sampling in shift-invariant spaces. This is accomplished by bringing together wavelet theory, frame theory, reproducing kernel Hilbert spaces, approximation theory, amalgam spaces, and sampling. This combination simplifies some parts of the literature on sampling.

[0208] The description below contains several new results. In particular, all the well-known results are developed in weighted Lp-spaces. Extensions of frame theory and reproducing kernel Hilbert spaces to Banach spaces are discussed, and the connections between reproducing kernels in weighted Lp-spaces, Banach frames, and sampling are described.

[0209] The rest of this Example is organized as follows. Next section, section 2, introduces the relevant spaces for sampling theory and presents some of their properties. Weighted Lp-spaces and sequence spaces are defined. Wiener amalgam spaces are discussed and we also derive some convolution relations in the style of Young's inequalities. The weighted Lp-shift-invariant spaces are then introduced in, and some of their main properties are established. The sampling problem in weighted shift-invariant spaces is subsequently disclosed in section 3. Next, in section 4, some aspects of reproducing kernel Hilbert spaces and frame theory are reviewed. The discussion includes an extension of frame theory and reproducing kernel Hilbert spaces to Banach spaces. The connections between reproducing kernels in weighted Lp-spaces, Banach frames, and sampling are also discussed. Frame algorithms for the reconstruction of a function from its samples are discussed in section 5 and so are iterative reconstructions in section 6. In applications, a function f does not belong to a particular prescribed space V, in general. Moreover, even if the assumption that a function f belongs to a particular space V is valid, the samples off are not exact due to digital inaccuracy, or the samples are corrupted by noise when they are obtained by a real measuring device. For this reason, the results of the various reconstruction algorithms when the samples are corrupted by noise are presented in section 7 as well, which is an important issue in practical applications. Theorems and lemmas given in this disclosure can be and have been proved by the inventors. Interested reader can find detailed proof and related discussions in reference [f], which is hereby incorporated herein by reference in its entirety.

2. Relevant Spaces

[0210] Function Spaces. This section, among other things, provides the basic framework for treating nonuniform sampling in weighted shift-invariant spaces. The shift-invariant spaces under consideration are of the form 19$\begin{array}{cc}V\left(\phi \right)=\left\{\sum _{k\in Z^d}^{}{c}_k\phi \left(·-k\right)\right\}& \left(2.1\right)\end{array}$

[0211] where c=(ck)kεZd is taken from some sequence space and φ is the so-called generator of V(φ). Before it is possible to give a precise definition of shift-invariant spaces, one needs to study the convergence properties of the series 20$\left\{\sum _{k\in Z^d}^{}{c}_k\phi \left(·-k\right)\right\}.$

[0212] In the context of the sampling problem the function s in V(φ) must also be continuous. In addition, one may want to control the growth or decay at infinity of the functions in V(φ). Thus the generator φ and the associated sequence space cannot be chosen arbitrarily. To settle these questions, we first discuss weighted Lp-spaces with specific classes of weight functions, and we then develop the main properties of amalgam spaces. Only then we give a rigorous definition of a shift-invariant space and derive its main properties in next section. Shift-invariant spaces figure prominently in other areas of applied mathematics, notably in wavelet theory and approximately theory [33, 34]. Our presentation will be adapted to the requirements of sampling theory.

[0213] Weighted Lpv-Spaces. To model decay or growth of functions, we use weighted Lp-spaces [41]. A function f belongs to Lpv (d) with weight function v if vf belongs to Lp (d). Equipped with the norm ∥f∥Lpv=∥vf∥Lp, the space Lpv is a Banach space. If the weight function v grow rapidly as |x|→∞, then the functions in Lpv decay roughly at a corresponding rate. Conversely, if the weight function v decays rapidly, then the functions in Lpv may grow as |x|→∞.

[0214] In general, a weight function is just a nonnegative function v. We will use two special types of weight functions as examples. The weight functions denoted by X are always assumed to be continuous, symmetric, i.e., ω(x)=ω(−x), positive, and submultiplicative:

0<ω(x+y)≦ω(x)ω(y) ∀x, yεd  (2.2)

[0215] This submultiplicativity condition implies that ω(x)≧ω(0)≧1 for all xεd. For a technical reason, we impose the growth condition 21$\sum _{n=1}^{\infty }\frac{\log \omega \left(\mathrm{nk}\right)}{n^2}<\infty \forall k\in Z^d.$

[0216] Although most of the results do not require this extra condition on ω, we use it in Lemma 2.11. This condition, however, does not limit in any way the scope of the invention. For simplicity we refer to ω as a submultiplicative weight. A prototypical example is the Sobolev weight ω(x)=(1+|x|)α, with α>0. When ω=1, we obtain the usual Lp-spaces.

[0217] In addition, a weight function v is called moderate with respect to the submultiplicative weight ω or simply ω-moderate, if it is continuous, symmetric, and positive and satisfies C ω(−x)v(y)≧v(x+y) for all x, yεd. For instance, the weights v(x)=(1+|x|)β are moderate with respect to ω(x)=(1+|x|)α if and only if |β|≦α. If v is ω-moderate, then C ω(−x)v(x+y)≧v(y)=v(x+y−x), and it follows that: 22$\frac{1}{v\left(x+y\right)}\le C\omega \left(x\right)\frac{1}{v\left(y\right)}.$

[0218] Thus, the weight 1/v is also ω-moderate.

[0219] If v is ω-moderate, then a simple computation shows that

∥f(·−y)∥Lvp≦Cω(y)∥f∥Lpv

[0220] and in particular, ∥f(·−y)∥Lpω≦ω(y)∥f∥Lpω. Conversely, if Lpv is translation-invariant, then 23$\omega \left(x\right)=\underset{{\&LeftDoubleBracketingBar;f\&RightDoubleBracketingBar;}_{L_v^p}\le 1}{\sup }{\&LeftDoubleBracketingBar;f\left(·-x\right)\&RightDoubleBracketingBar;}_{L_v^p}$

[0221] is submultiplicative and v is ω-moderate. To see this, we note that 24$\omega \left(x\right)=\underset{{\&LeftDoubleBracketingBar;f\&RightDoubleBracketingBar;}_{L_v^p}\le 1}{\sup }{\&LeftDoubleBracketingBar;f\left(·-x\right)\&RightDoubleBracketingBar;}_{L_v^p}$

[0222] is the operator norm of the translation operator f→f (·−x). Since operator norms are submultiplicative, it follows that ω(x))ω(y)≧ω(x+y). Moreover, 25${\int }_{R^d}{\&LeftBracketingBar;f\left(t-x\right)\&RightBracketingBar;}^p{v\left(t\right)}^{p}dt={\int }_{R^d}{\&LeftBracketingBar;f\left(t\right)\&RightBracketingBar;}^p{v\left(t+x\right)}^{p}dt\le {\omega \left(x\right)}^p{\int }_{R^d}{\&LeftBracketingBar;f\left(t\right)\&RightBracketingBar;}^p{v\left(t\right)}^{p}dt$

[0223] Thus, ω(x) v(y)≧v(x+y), and therefore the weighted Lp-spaces with a moderate weight are exactly the translation-invariant spaces.

[0224] We also consider the weighted sequence spaces lpv(Zd) with weight v: a sequence {(ck):kεZd} belongs to lpv if ((c v)k)=(ckvk) belongs to lp with norm ∥c∥lPv=∥vc∥lp, where (vk) is the restriction of v to Zd.

[0225] Wiener Amalgam Spaces. For the sampling problem we also need to control the local behavior of functions so that the sampling operation f→(f(xj))jεJ is at least well defined. This is done conveniently with the help of the Wiener amalgam spaces W(Lpv). These consist of functions that are “locally in L∞ and globally in Lpv.” To be precise, a measurable function f belongs to W(Lpv), 1≦p<∞, if it satisfies 26$\begin{array}{cc}{\&LeftDoubleBracketingBar;f\&RightDoubleBracketingBar;}_{W\left(L_v^p\right)}^p=\sum _{k\in Z^d}^{}\mathrm{ess}\sup \left\{{\&LeftBracketingBar;f\left(x+k\right)\&RightBracketingBar;}^p{v\left(k\right)}^p;x\in {\left[0,1\right]}^d\right\}<\infty & \left(2.3\right)\end{array}$

[0226] If p=∞, a measurable function f belongs to W(L∞v) if it satisfies 27$\begin{array}{cc}{\&LeftDoubleBracketingBar;f\&RightDoubleBracketingBar;}_{W\left(L_v^{\infty }\right)}=\underset{k\in Z^d}{\sup }\left\{\mathrm{ess}\sup \left\{\&LeftBracketingBar;f\left(x+k\right)\&RightBracketingBar;v\left(k\right);x\in {\left[0,1\right]}^d\right\}\right\}<\infty .& \left(2.4\right)\end{array}$

[0227] Note that W(L∞v) coincides with L∞v.

[0228] Endowed with this norm, W(Lpv) becomes a Banach space [43, 45]. Moreover, it is translation invariant; i.e., if fεW(Lpv), then f(·−y)εW(Lpv) and

∥f(·−y)∥W(Lpv)≦Cω(y)∥f∥W(Lpv).

[0229] The subspace of continuous functions W0(Lpv)=W(C, Lpv)⊂W(Lpv) is a closed subspace of W(Lpv) and thus also a Banach space [43, 45]. We have the following inclusions between the various spaces.

[0230] THEOREM 2.1. Let v be ω-moderate and 1≦p≦q≦∞. Then the following inclusions hold:

[0231] (i) W0(Lvp)⊂W0(Lvq) and W(Lvp)⊂W(Lvq)⊂Lvp

[0232] (ii) W0(Lvp)⊂W0(Lvp), W(Lωp)⊂W(Lvp), and Lωp⊂Lvp

[0233] The following convolution relations in the style of Young's theorem [118] are useful.

[0234] THEOREM 2.2. Let v be ω-moderate.

[0235] (i) If fεLvp and gεL1ω, then f*gεLvp and ∥f*g∥Lvp≦C∥f∥Lvp∥g∥Lω1

[0236] (ii) If fεLvp and gεW(Lω1), then f*gεW(Lvp) and ∥f*g∥W(Lvp)≦C∥f∥Lvp∥g∥W(Lφ7 1)

[0237] (iii) If cεlvp and dεlω1, then c*dεlcp and ∥c*d∥lvp≦C∥c∥lvp∥d∥lω1

[0238] Remark 2.1. Amalgam spaces and their generalizations have been investigated by Feichtinger, and the results of Theorem 2.1 can be found in [42, 43, 45, 44]. The results and methods developed by Feichtinger can also be used to deduce Theorem 2.2. However, for the sake of completeness, in section 2.4 we present direct proofs of Theorems 2.1 and 2.2 that do not rely on the deep results of amalgam spaces

[0239] Shift-Invariant Spaces. Among other things, this section discusses shift-invariant spaces and their basic properties. Although some of the following observations are known in wavelet and approximately theory, they have received little attention in connection with sampling.

[0240] Given a so-called generator φ, we consider shift-invariant spaces of the form 28$\begin{array}{cc}{V}_v^p\left(\phi \right)=\left\{\sum _{k\in Z^d}^{}{c}_k\phi \left(·-k\right):c\in l_v^p\right\}& \left(2.5\right)\end{array}$

[0241] If v=1, we simply write Vp(φ). The weight function v controls the decay or growth rate of the functions in Vpv(φ). To some extent, the parameter p also controls the growth of the functions in Vpv(φ) but more importantly, p controls the norm we wish to use for measuring the size of our functions. For some applications in image processing, the choice p=1 is appropriate [36]; p=2 corresponds to the energy norm, and p=∞ is used as a measure in some quality control applications. Moreover, the smoothness of a function and its appropriate value of p, 1≦p<∞, for a given class of signals or images can be estimated using wavelet decomposition techniques [36]. The determination of p and the signal smoothness are used for optimal compression and coding of signals and images.

[0242] For the spaces Vpv(φ) to be well defined, some additional conditions on the generator φ must be imposed. For v=1 and p=2, the standard condition in wavelet theory is often stated in the Fourier domain as 29$\begin{array}{cc}0<m\le {\hat{a}}_{\phi }\left(\xi \right)=\sum _{j\in Z^d}^{}{\&LeftBracketingBar;\hat{\phi }\left(\xi +j\right)\&RightBracketingBar;}^2\le M<\infty \mathrm{for}\mathrm{almost}\mathrm{every}\xi ,\mathrm{for}& \left(2.6\right)\end{array}$

[0243] some constants m>0 and M>0 [80, 81]. This condition implies that V2(φ) is a closed subspace of L2 and that {φ(·−k):kεZd} is a Riesz basis of V2(φ), i.e., the image of an orthonormal basis under an invertible linear transformation [33].

[0244] The theory of Riesz bases asserts the existence of a dual basis. Specifically, for any Riesz basis for V2(φ) of the form {φ(·−k):kεZd}, there exists a unique function {tilde over (φ)}εV2(φ) such that {{tilde over (φ)}(·−k):kεZd} is also a Riesz basis for V2(φ) and such that {tilde over (φ)} satisfies the biorthogonality relation

{tilde over (φ)}(·), φ(·−k)=δ(k)

[0245] Where δ(0)=1 and δ(k)=0 for k≠0. Since the dual generator {tilde over (φ)} belongs to V2(φ), it can be expressed in the form 30$\begin{array}{cc}\tilde{\phi }(·)=\sum _{k\in Z^d}^{}{b}_k\phi \left(·-k\right).& \left(2.7\right)\end{array}$

[0246] The coefficients (bk) are determined explicitly by the Fourier series 31$\sum _{k\in Z^d}{b}_k{e}^{2\pi ik\xi }={\left(\sum _{k\in Z^d}{\&LeftBracketingBar;\hat{\phi }\left(\xi +k\right)\&RightBracketingBar;}^2\right)}^{-1};$

[0247] i.e., (bk) is the inverse Fourier transform of 32${\left(\sum _{k\in Z^d}{\&LeftBracketingBar;\hat{\phi }\left(\xi +k\right)\&RightBracketingBar;}^2\right)}^{-1}$

[0248] (see, for example, [8, 9]. Since aφ(ξ)−1≦1/m by (2.6), the sequence (bk) exists and belongs to l2(Zd).

[0249] In order to handle general shift-invariant spaces Vpv(φ) instead of V2(φ), we need more information about the dual generator. The following result is one of the central results in this disclosure and is useful for the treatment of general shift-invariant spaces.

[0250] THEOREM 2.3. Assume that φεW(L1ω) and that {φ(·−k):kεZd} is a Riesz basis for V2(φ). Then the dual generator {tilde over (φ)} is in W(L1ω).

[0251] As a corollary, we obtain the following theorem.

[0252] THEOREM 2.4. Assume that φεW(L1ω) and that v is ω-moderate.

[0253] (i) The space Vpv(φ) is a subspace (not necessarily closed) of Lpv and W(Lpv) for any p with 1≦p≦∞.

[0254] (ii) If {φ(·−k):kεZd} is a Riesz basis of V2(φ), then there exist constants mp>0, Mp>0 such that 33$\begin{array}{cc}{m}_p{\&LeftDoubleBracketingBar;c\&RightDoubleBracketingBar;}_{l_v^p}\le {\&LeftDoubleBracketingBar;\sum _{k\in Z^d}{c}_k\phi \left(·-k\right)\&RightDoubleBracketingBar;}_{L_v^p}\le M_p{\&LeftDoubleBracketingBar;c\&RightDoubleBracketingBar;}_{l_v^p}\forall c\in l_v^p\left(Z^d\right)& \left(2.8\right)\end{array}$

[0255] is satisfied for all 1≦p≦∞ and all ω-moderate weights v. Consequently, {φ(·−k):kεZd} is an unconditional basis for Vpv(φ) for 1≦p≦∞, and Vpv(φ) is a closed subspace of Lpv and W(Lpv) for 1≦p≦∞.

[0256] The theorem says that the inclusion in Theorem 2.4(i) and the norm equivalence (2.8) hold simultaneously for all p and all ω-moderate weights, provided that they hold for the Hilbert space V2(φ). But in V2(φ), the Riesz basis property (2.8) is much easier to check. In fact, it is equivalent to inequalities (2.6). In equalities (2.8) imply that Lpv and Vpv(φ) are isomorphic Banach spaces and that the set {φ(·−k):kεZd} is an unconditional basis of Vpv(φ). In approximation theory we say that φ has stable integer translates and is a stable generator [67, 68, 69]. When v=1 the conclusion (2.8) of Theorem 2.4 is well known and can be found in [67, 68, 69].

[0257] As a corollary of Theorem 2.4, we obtain the following inclusions among shift-invariant spaces.

[0258] COROLLARY 2.5. Assume that φεW(L1ω) and that v is ω-moderate. Then

V ω 1 (φ)⊂Vωp(φ)⊂Vωq(φ) for 1≦p≦q≦∞

[0259] and

V ω q (φ)⊂Vωp(φ) for 1≦q≦∞.

[0260] 3. The Sampling Problem in Weighted Shift-Invariant Spaces Vpv(φ)

[0261] For a reasonable formulation of the sampling problem in Vpv(φ) the point evaluations f→f(x) must be well defined. Furthermore, a small variation in the sampling point should produce only a small variation in the sampling value. As a minimal requirement, we need the functions in Vpv(φ) to be continuous, which is guaranteed by the following statement.

[0262] THEOREM 3.1. Assume that φεW0(L1ω), that φ satisfies (2.6), and that v is −ω-moderate.

[0263] (i) Vvp(φ)⊂W0(Lvp) for all p,1≦p≦∞

[0264] (ii) If fεVvp(φ), then we have the norm equivalences ∥f∥Lvp≈∥c∥lvp≈∥f∥W(Lvp)

[0265] (iii) If X={Xj:jεJ} is such that in fj,l|xj−xl∥>0, then 34$\begin{array}{cc}{\left(\sum _{x_k\in X}{\&LeftBracketingBar;f\left(x_k\right)\&RightBracketingBar;}^p{\&LeftBracketingBar;v\left(x_k\right)\&RightBracketingBar;}^p\right)}^{\frac{1}{p}}\le C_p{\&LeftDoubleBracketingBar;f\&RightDoubleBracketingBar;}_{L_v^p}\forall f\in V_v^p\left(\phi \right).& \left(3.1\right)\end{array}$

[0266] In particular, if φ is continuous and has compact support, then the conclusions (i)-(iii) hold.

[0267] A set X={xj:jεJ} satisfying in fj,l|xj−xl|>0 is called separated.

[0268] Inequality (3.1) has two interpretations. It implies that the sampling operator SX:f→f|x is a bounded operator from Vpv(φ) into the corresponding sequence space 35$l_v^p\left(X\right)=\left\{\left(c_j\right)\text{:}{\left(\sum _{j\in J}{\&LeftBracketingBar;\left(c_j\right)\&RightBracketingBar;}^p{v\left(x_j\right)}^p\right)}^{\frac{1}{p}}<\infty \right\}$

[0269] Equivalently, the weighted sampling operator SX:fv|X is a bounded operator from Vpv(φ) into Lp.

[0270] To recover a function fεVpv(φ) from its samples, we need a converse of inequality (3.1). Following Landau [74], we say that X is a set of sampling for Vpv(φ) if 36$\begin{array}{cc}{c}_p{\&LeftDoubleBracketingBar;f\&RightDoubleBracketingBar;}_{L_v^p}\le {\left(\sum _{x_j\in X}{\&LeftBracketingBar;f\left(x_j\right)\&RightBracketingBar;}^p{\&LeftBracketingBar;v\left(x_j\right)\&RightBracketingBar;}^p\right)}^{\frac{1}{p}}\le C_p{\&LeftDoubleBracketingBar;f\&RightDoubleBracketingBar;}_{L_v^p}& \left(3.2\right)\end{array}$

[0271] where cp and Cp are positive constants independent of f.

[0272] The left-hand inequality implies that if f(xj)=0 for all xjεX, then f=0. Thus X is a set of uniqueness. Moreover, the sampling operator SX can be inverted on its range and SX−1 is a bounded operator from Range (SX)⊂lpv(X) to Vpv(φ). Thus (3.2) says that a small change of a sampled value f(xj) causes only a small change of f. This implies that the sampling is stable or, equivalently, that the reconstruction of f from its samples is continuous. As pointed out in section 1, every set of sampling is a set of uniqueness, but the converse is not true. For practical considerations and numerical implementations, only sets of sampling are of interest, because only these can, lead to robust algorithms.

[0273] A solution to the sampling problem may include, among other things, the following two parts:

[0274] (a) Given a generator φ, we need to find conditions on X, usually in the form of a density, such that the norm equivalence (3.2) holds. Then, at least in principle, fεVpv(φ) is uniquely and stably determined by f|X.

[0275] (b) We need to design reconstruction procedures that are useful and efficient in practical applications. The objective is to find efficient and fast numerical algorithms that recover f from its samples f|X, when (3.2) is satisfied.

[0276] Remark 3.1.

[0277] (i) The hypothesis that X be separated is for convenience only and is not essential. For arbitrary sampling sets, we can use adaptive weights to compensate for the local variations of the sampling density [48, 49]. Let Vj={xεd:|x−xj|≦|x−xk | for all k≈j} be the Voronoi region at xj, and let γj=λ(Vj) be the size of Vj. Then X is a set of sampling for Vpv(φ) if 37$\begin{array}{cc}{c}_p{\&LeftDoubleBracketingBar;f\&RightDoubleBracketingBar;}_{L_v^p}\le {\left(\sum _{x_j\in X}{\&LeftBracketingBar;f\left(x_j\right)\&RightBracketingBar;}^p{\gamma }_j{\&LeftBracketingBar;v\left(x_j\right)\&RightBracketingBar;}^p\right)}^{\frac{1}{p}}\le C_p{\&LeftDoubleBracketingBar;f\&RightDoubleBracketingBar;}_{L_v^p}& \end{array}$

[0278]  In numerical applications the adaptive weights γj are used as a cheap device for preconditioning and for improving the ratio Cp/cp, the condition number of the set of sampling [49, 101].

[0279] (ii) The assumption that the samples {f(xj):jεJ} can be measured exactly is not realistic. A better assumption is that the sampled data is of the form 38$\begin{array}{cc}{g}_{x_j}={\int }_{R^d}f\left(x\right){\psi }_{x_j}\left(x\right)dx& \left(3.3\right)\end{array}$

[0280]  where {ψxj:xjεX} is a set of functionals that act on the function f to produce the data {gxj:xjεX}. The functionals {ψxj:xjεX} may reflect the characteristics of the sampling devices. For this case, the well-posedness condition (3.2) must be replaced by 39$\begin{array}{cc}{c}_p{\&LeftDoubleBracketingBar;f\&RightDoubleBracketingBar;}_{L^p}\le {\left(\sum _{x_j\in X}{\&LeftBracketingBar;g_{x_j}\left(f\right)\&RightBracketingBar;}^p\right)}^{\frac{1}{p}}\le C_p{\&LeftDoubleBracketingBar;f\&RightDoubleBracketingBar;}_{L^p}& \left(3.4\right)\end{array}$

[0281]  where gxj are defined by (3.3) and where cp and Cp are positive constants independent off [1].

[0282] 4. Reproducing Kernel Hilbert Spaces, Frames, and Nonuniform Sampling

[0283] As disclosed above, results of Paley and Wiener and Kadec relate Riesz bases consisting of complex exponentials to sampling sets that are perturbations of Z. More generally, the appropriate concept for arbitrary sets of sampling in shift-invariant spaces is the concept of frames discussed below. Frame theory generalizes and encompasses the theory of Riesz bases and enables us to translate the sampling problem into a problem of functional analysis. The connection between frames and sets of sampling is established by means of reproducing kernel Hilbert spaces (“RKHSs”), discussed below. Frames are introduced in next section, and the relation between RKHSs frames, and sets of sampling is developed subsequently.

[0284] RKHSs. Theorem 3.1 (iii) holds for arbitrary separated sampling sets, so in particular Theorem 3.1 (iii) shows that all point evaluations f→f(x) are continuous linear functionals on Vpv(φ) for all xεd. Since Vpv(φ)⊂Lpv and the dual space of Lpv is Lp′1/v, where 1/p+1/p′=1, there exists a function KxεLpv such that 40$f\left(x\right)=⟨f,K_x⟩={\int }_{R^d}f\left(t\right)\stackrel{\_}{K_x\left(t\right)}dt$

[0285] for all fεVpv(φ). In addition, it will be shown that KxεVp1/v(φ)

[0286] In the case of a Hilbert space H of continuous functions on d, such as V2(φ), the following terminology is used. A Hilbert space is an RKHS [117] if, for any xεd, the pointwise evaluation f→f(x) is a bounded linear functional on H. The unique functions Kx→H satisfying f(x)=(f, Kx) are called the reproducing kernels of H.

[0287] With this terminology we have the following consequence of Theorem 3.1.

[0288] THEOREM 4.1. Let v be ω-moderate. If φεW0(L1ω), then the evaluations f→f(x) are continuous functionals, and there exist functions KxεV1ω(φ) such that f(x)=f, Kx. The kernel functions are given explicitly by 41$\begin{array}{cc}{K}_x\left(y\right)=\sum _{k\in Z^d}^{}\stackrel{\_}{\phi \left(x-k\right)}\tilde{\phi }\left(y-k\right)& \left(4.1\right)\end{array}$

[0289] In particular, V2(φ) is an RKHS.

[0290] The above theorem is a reformulation of Theorem 3.1. We only need to prove the formula for the reproducing kernel. Note that Kx in (4.1) is well defined: since {tilde over (φ)}εW0(L1ω), Theorem 2.3 combined with Theorem 3.1 (iii) implies that the sequence {{tilde over (φ)}(x−k):kεZd} belongs to l1ω. Thus, by the definition of V1ω(φ), we have KxεV1ω(φ), and so KxεVpv(φ) for any p with 1≦p≦∞ and any ω-moderate weight v. Furthermore, Kx is clearly the reproducing kernel, because if f(x)=ΣkCKφ(x−k), then 42$⟨f,K_x⟩=\sum _{j,k}^{}{c}_j\phi \left(x-k\right)⟨\phi \left(·-j\right),\tilde{\phi }\left(·-k\right)⟩=\sum _k^{}{c}_k\phi \left(x-k\right)=f\left(x\right)$

[0291] Frames. In order to reconstruct a function fεVpv(φ) from its samples f(xj), it is sufficient to solve the (infinite) system of equations 43$\begin{array}{cc}\sum _{k\in Z^d}^{}{c}_k\phi \left(x_j-k\right)=f\left(x_j\right)& \left(4.2\right)\end{array}$

[0292] for the coefficients (ck). If we introduce the infinite matrix U with entries

U jk =φ(xj−k)  (4.3)

[0293] indexed by X×Zd, then the relation between the coefficient sequence c and the samples is given by

Uc=f|x.

[0294] Theorem 3.1(ii) and (iii) imply that f|Xεlpv(X). Thus U maps lpv(Zd) into lpv(X).

[0295] Since f(x)=f,Kx, the sampling inequality (3.2) implies that the set of reproducing kernels {Kxj, xjεX} spans Vp′1/V. This observation leads to the following abstract concepts.

[0296] A Hilbert frame (or simply a frame) {ej:jεJ} of a Hilbert space H is a collection of vectors in H indexed by a countable set J such that 44$\begin{array}{cc}A{\&LeftDoubleBracketingBar;f\&RightDoubleBracketingBar;}_H^2\le \sum _j^{}{\&LeftBracketingBar;⟨f,e_j⟩\&RightBracketingBar;}^2\le B{\&LeftDoubleBracketingBar;f\&RightDoubleBracketingBar;}_H^2& \left(4.4\right)\end{array}$

[0297] for two constants A, B>0 independent of fεH[40].

[0298] More generally, a Banach frame for a Banach space B is a collection of functionals {ej:jεJ}⊂B* with the following properties [54].

[0299] (a) There exists an associated sequence space Bd on the index set J, such that 45$A{\&LeftDoubleBracketingBar;f\&RightDoubleBracketingBar;}_B\le {\&LeftDoubleBracketingBar;{\left(⟨f,e_j⟩\right)}_{j\in J}\&RightDoubleBracketingBar;}_{B_d}\le B{\&LeftDoubleBracketingBar;f\&RightDoubleBracketingBar;}_B$

[0300]  for two constants A,B>0 independent of fεB.

[0301] (b) There exists a so-called reconstruction operator R from Bd into B such that 46$R\left({\left(⟨f,e_j⟩\right)}_{j\in J}\right)=f$

[0302] Relations between RKHSs, Frames, and Nonuniform Sampling. The following theorem translates the different terminologies that arise in the context of sampling theory [2, 40, 74].

[0303] THEOREM 4.2 The following are equivalent:

[0304] (i) X={xj:jεJ} is a set of sampling for Vpv(φ).

[0305] (ii) For the matrix U in (4.3), there exist a, b>0 such that 47$a{\&LeftDoubleBracketingBar;c\&RightDoubleBracketingBar;}_{l_v^p}\le {\&LeftDoubleBracketingBar;\mathrm{Uc}\&RightDoubleBracketingBar;}_{l_v^p\left(X\right)}\le b{\&LeftDoubleBracketingBar;c\&RightDoubleBracketingBar;}_{l_v^p}\forall c\in l_v^p$

[0306] (iii) There exist positive constants a>0 and b>0 such that 48$a{\&LeftDoubleBracketingBar;f\&RightDoubleBracketingBar;}_{L_v^p}\le {\left(\sum _{x_j\in X}^{}{\&LeftBracketingBar;f\left(x_j\right)\&RightBracketingBar;}^p{\&LeftBracketingBar;v\left(x_j\right)\&RightBracketingBar;}^p\right)}^{1/p}\le b{\&LeftDoubleBracketingBar;f\&RightDoubleBracketingBar;}_{L_v^p}\forall f\in V_v^p\left(\phi \right)$

[0307] (iv) For p=2, the set of reproducing kernels {Kxj:xjεX} is a (Hilbert) frame for V2(φ).

[0308] Remark 4.1.

[0309] (i) The relation between RKHSs and uniform sampling of bandlimited functions was first reported by Yao [170] and used to derive interpolating series similar to (1.1). For the case of shift-invariant spaces, this connection was established in (9). Sampling for functions in RKHSs was studied in [84]. For the general case of nonuniform sampling in shift-invariant spaces, the connection was established in [5].

[0310] (ii) The relation between Hilbert frames and sampling of bandlimited functions is well known [14, 48]. Sampling in shift-invariant spaces is more recent, and the relation between frames and sampling in shift-invariant spaces (with p=2 and v=1) can be found in [5, 30, 75, 77, 102].

[0311] (iii) The relation between Hilbert frames and the weighted average sampling mentioned in Remark 3.1 can be found in [1]. This relation is obtained via kernels that generalize the RKHS.

5. Frame Algorithms for Lpv-Spaces

[0312] Theorem 4.2 states that a separated set X={xj:jεJ} is a set of sampling for V2(φ) if and only if the set of reproducing kernels {{tilde over (K)}xj:xjεX} is a frame for V2(φ). It is well known from frame theory that there exists a dual frame {{tilde over (K)}xj:xjεX} that allows us to reconstruct the function f⊂V2(φ) explicitly as 49$\begin{array}{cc}f\left(x\right)=\sum _{j\in J}^{}⟨f,K_{x_j}⟩{\tilde{K}}_{x_j}\left(x\right)=\sum _{j\in J}^{}f\left(x_j\right){\tilde{K}}_{x_j}\left(x\right)& \left(5.1\right)\end{array}$

[0313] However, a dual frame {{tilde over (K)}xj:xjεX} is difficult to find in general, and this method for recovering a function fεV2(φ) from its samples {f(xj):xjεX} is often not practical.

[0314] Instead, the frame operator 50$\begin{array}{cc}\mathrm{Tf}\left(x\right)=\sum _{j\in J}^{}⟨f,K_{x_j}⟩K_{x_j}\left(x\right)=\sum _{j\in J}^{}f\left(x_j\right)K_{x_j}\left(x\right)& \left(5.2\right)\end{array}$

[0315] can be inverted via an iterative that we now describe. The operator 51$I-\frac{2}{A+B}T$

[0316] is contractive, i.e., the operator norm on L2(d) satisfies the estimate 52${\&LeftDoubleBracketingBar;I-\frac{2}{A+B}T\&RightDoubleBracketingBar;}_{\mathrm{op}}\le \frac{B-A}{A+B}<1$

[0317] where A, B are frame bounds for {Kxj:xjεX}. Thus, 53$\frac{2}{A+B}T$

[0318] can be inverted by the Neumann series 54$\frac{A+B}{2}{T}^{-1}=\sum _{n=0}^{\infty }{\left(I-\frac{2}{A+B}T\right)}^n$

[0319] This analysis gives the iterative frame reconstruction algorithm, which is made up of an initialization step 55$f_1=\sum _{j\in J}^{}f\left(x_j\right)K_{x_j}$

[0320] and iteration 56$\begin{array}{cc}{f}_n=\frac{2}{A+B}{f}_1+\left(I-\frac{2}{A+B}T\right)f_{n-1}& \left(5.3\right)\end{array}$

[0321] As n→∞, the iterative frame algorithm (5.3) converges to f∞=T−1f1=T1Tf=f.

[0322] Remark 5.1.

[0323] (i) The computation of T requires the computation of the reproducing frame functions {Kxj:xjεX}, which is a difficult task. Moreover, for each sampling set X we need to compute a new set of reproducing frame functions {Kxj:xjεX}.

[0324] (ii) Even if the frame functions {Kxj:xjεX} are known, the performance of the frame algorithm depends sensitively on estimates for the frame bounds.

[0325] Since accurate and explicit frame bounds, let alone optimal ones, are hardly ever known for nonuniform sampling problems, the frame algorithm converges very slowly in general. For efficient numerical computations involving frames, the primitive iteration (5.3) should therefore be replaced almost always by acceleration methods, such as Chebyshev or conjugate gradient acceleration. In particular, conjugate gradient methods converge at the optimal rate, even without any knowledge of the frame bounds [49, 56].

[0326] (iii) The convergence of the frame algorithm is guaranteed only in L2, even if the function belongs to other spaces Lpv. It is a remarkable fact that in Hilbert space the norm equivalence (4.4) alone guarantees that the frame operator is invertible. In Banach spaces the situation is much more complicated and the existence of a reconstruction procedure must be postulated in the definition of a Banach frame. In the special case of sampling in shift-invariant spaces, the frame operator T is invertible on all Vpv(φ) whenever T is invertible on V2(φ) and φ possesses a suitable polynomial decay [57].

6. Iterative Reconstruction Algorithms

[0327] Since the iterative frame algorithm is often slow to converge and its convergence is not even guaranteed beyond V2(φ), alternative reconstruction procedures have been designed [4, 76]. These procedures are also iterative and based on a Neumann series. New iterative reconstruction algorithms developed according to the present invention are disclosed and discussed below.

[0328] In one embodiment, as shown in FIG. 5, at step 503, samples {f(xj):xjεX} are obtained for function f of an object of interest. At step 505, the function f is approximated from its samples {f(xj):xjεX} using an interpolation or a quasi-interpolation Qxf. For examples, Qxf could be a piecewise linear interpolation of the samples f|X or even an approximation by step functions, the so-called sample-and-hold interpolant.

[0329] At step 507, the approximation Qxf is then projected in the space Vpv(φ) to obtain the first approximation f1=PQXfεVpv(φ), which then can be fed back, at step 501, to sampling operator Sx and produce an f1(X). The error e=f−f1 between the functions f and f1 belongs to the space Vpv(φ). Moreover, the values of e on the sampling set X can be calculated, at step 505, from {f(xj):xjεX} and (PQXf)(xj). Then we repeat the interpolation projection procedure on e and obtain a correction e1. The updated estimate is now, at step 509, f2=f1+e1, which can then be fed back to step 501 for a new round of calculation. By repeating this procedure, we obtain a sequence fn=f1+e1+e2+e3+ . . . +en−1 that converges to the function f.

[0330] In order to prove convergence results for this type of algorithm, we need the sampling set to be dense enough. The appropriate definition for the sampling density of X is again due to Beurling.

[0331] DEFINITION 6.1. A set X={xj:jεJ} is ε0-dense in d if 57$\begin{array}{cc}{R}^d=\bigcup _j{B}_{\gamma }\left(x_j\right)\forall \gamma >{\gamma }_0.& \left(6.1\right)\end{array}$

[0332] This definition implies that the distance of any sampling point to its nearest neighbor is at most 2γ0. Thus, strictly speaking, ε0 is the inverse of a density; i.e., if ε0 increases, the number of points per unit cube decreases. In fact, if a set X is ε0-dense, then its Beurling density defined by (1.2) satisfies D(X)≧γ0−1. This last relation states that Δ0-density imposes more constraints on a sampling set X than the Beurling density D(X).

[0333] To create suitable quasi-interpolants, we proceed as follows: Let {βj}jεJ be a partition of unity such that

[0334] (1)0<βj≦1 for all jεJ;

[0335] (2) supp βj⊂Bγ(xj); and

[0336] (3) ΣjεJβj=1.

[0337] A partition of unity that satisfies these conditions is sometimes called a bounded partition of unity. Then the operator QX defined by

Q X f=Σ jεJ f (xj)βj

[0338] is a quasi-interpolant of the sampled values f|X.

[0339] In this situation we have the following qualitative statement.

[0340] THEOREM 6.1. Let φ in W0(L1ω) and let P be a bounded projection from Lpv onto Vpv(φ). Then there exists a density γ>0 (γ=γ(v, p, P)) such that any fεVpv(φ) can be recovered from its samples {f(xj):xjεX} on any γ-dense set X={xj:jεJ} by the iterative algorithm 58$\begin{array}{cc}\left\{\begin{array}{c}{f}_1={\mathrm{PQ}}_{X}f\\ f_{n+1}={\mathrm{PQ}}_X\left(f-f_n\right)+f_n\end{array}\right\}& \left(6.2\right)\end{array}$

[0341] Then iterates fn converge to f uniformly and in the W(Lpv)- and Lpv-norms. The convergence is geometric, that is,

∥f−fn∥Lvp≦C∥f−fn∥W(Lvp)≦C′∥f−f1∥W)Lvp)αn

[0342] For some α=α(γ)<1.

[0343] The algorithm based on this iteration and as illustrated in FIG. 5 has been successfully utilized in practice. A reconstruction of a function f with ∥f∥2=3.5 using the iteration algorithm as set forth in FIG. 5 according to one embodiment of the present invention is shown in FIG. 6. Referring now to FIG. 6(a), curve 601 represents function f belonging to the shift-invariant space generated by the Gaussian function exp(−x2/σ2), where σ is about 0.81, and its sample values 603, {f(xj), xjεX}, marked by + (with density γ about 0.8). As shown in FIG. 6(b), curve 611 represents error ∥f−fn∥L2 against the number of iterations, indicating good convergence of the reconstruction. As shown in FIG. 6(c), curve 621 represents final error f−fn after 10 iterations. And as shown in FIG. 6(d), curve 631 represents reconstructed function f10 (continuous line) and original samples 603 {f(xj), xjεX}, showing great match between the two.

[0344] Additionally, a missing data reconstruction of an MRI image performed according to one embodiment of the present invention is shown in FIG. 7. Referring now to FIG. 7(a), original digital MRI image 701 is shown with 128×128 samples. As shown in FIG. 7(b), MRI image 703 is given with 50% randomly missing samples. As shown in FIG. 7(c), reconstruction image 711 with 50% randomly missing samples using the iterative reconstruction algorithm according to one embodiment of the present invention is given and shows great improvement(s) over image 703. The corresponding shift invariant space is generated by φ(x, y)=β3(x)×β3(y), where β3=χ[0,1]*χ[0,1]* χ[0,1]*, and χ[0,1] is the B-spline function of degree 3.

[0345] Remark 6.1. For v=1, a special case, Theorem 6.1 was proved in [4]. For a special case of the weighted average sampling mentioned in Remark 3.1, a modified iterative algorithm and a theorem similar to Theorem 6.1 can be found in [1].

[0346] Universal Projections in Weighted Shift-Invariant Spaces. Theorem 6.1 requires bounded projections from Lpv onto Vpv(φ). In contrast to the situation in Hilbert space, the existence of bounded projections in Banach spaces is a difficult problem. In the context of nonuniform sampling in shift-invariant spaces, we would like the projections to satisfy additional requirements. In particular, we would like projectors that can be implemented with fast algorithms. Further, it would be useful to find a universal projection, i.e., a projection that works simultaneously for all Lpv, 1≦p≦∞, and all weights v. In shift-invariant spaces such universal projections do indeed exist.

[0347] THEOREM 6.2. Assume φεW0(L1ω). Then the operator 59$P:f\to \sum _{k\in Z^d}^{}⟨f,\tilde{\phi }\left(·-k\right)⟩\phi \left(·-k\right)$

[0348] is a bounded projection from Lpv onto Vpv(φ) for all p, 1≦p≦∞, and all ω-moderate weights v.

[0349] Remark 6.2. The operator P can be implemented using convolutions and sampling. Thus the universal projector P can be implemented with fast “filtering” algorithms [3].

7. Data Reconstruction in Presence of Noise

[0350] In practical applications the given data are rarely the exact samples of a function fεVpv(φ). In other words, we generally do not have “perfect” data for an object of interest that can be represented by a function f but data with noise. We assume more generally that f belongs to W0(Lpv); then the sampling operator f→{f(xj):xjεX} still makes sense and yields a sequence in lpv(X). Alternatively, we may assume that fεVpv(φ), but that the sampled sequence is a noisy version of {f(xj):xjεX}, e.g., that the sampling sequence has the form {f′xj=f(xj)+ηj}εlpv(X). If a reconstruction algorithm is applied to noisy data, then the question arises whether the algorithm still converges, and if it does, to which limit it converges.

[0351] To see what is involved, we first consider sampling in the Hilbert space V2(φ). Assume that X={Xj:jεJ} is a set of sampling for V2(φ). Then the set of reproducing kernels {Kxj:xjεX} forms a frame for V2(φ), and so fεV2(φ) can be reconstructed from the samples f(xj)=(f, Kxj) with the help of the dual frame {{tilde over (K)}xj:xjεX}⊂V2(φ) in the form of the expansion 60$\begin{array}{cc}f=\sum _{j\in J}^{}⟨f,K_{x_j}⟩{\tilde{K}}_{x_j}=\sum _{j\in J}^{}f\left(x_j\right){\tilde{K}}_{x_j}& \left(7.1\right)\end{array}$

[0352] If f∉V2(φ), but fεW0(L2), say, then f(xj)≠f,Kxj in general. However, the coefficients (f,Kxj) still make sense for fεL2 and the frame expansion (7.1) still converges. The following result describes the limit of this expansion when f∉V2(φ).

[0353] THEOREM 7.1. Assume that X⊂d is a set of sampling for V2(φ) and let P be the orthogonal projection from L2 onto V2(φ). Then 61$Pf=\sum _{j\in J}⟨f,K_{x_j}⟩{\tilde{K}}_{x_j}$

[0354] for all fεL2.

[0355] The previous theorem suggests a procedure for sampling: the function f is first “prefiltered” with the reproducing kernel Kx to obtain the function fa defined by f(x)=f,Kx for all xεd. Sampling fa on X then gives a sequence of inner products fa(xj)=f,Kxj. The reconstruction (7.1) off, is then the least square approximation off by a function faεV2(φ). In the case of bandlimited functions, we have φ(x)=sinπx/(πx) and Kx(t)=(sinπt(t−x))/(π(t−x)). Then the inner product fa(x)=f,Kx=f*φ(x) is just a convolution. The filtering operation corresponds to a restriction of the bandwidth to [−½, ½], because (f*φ){circumflex over ( )}={circumflex over (f)}·χ[−1/2,1/2], and is usually called prefiltering to reduce aliasing.

[0356] In practical situations, any sampling sequence may be perturbed by noise. This perturbation can be modeled in several equivalent ways: (a) The function fεV2(φ) is sampled on X, and then noise ηjεl2 is added, resulting in a sequence f′j=f(xj)+ηj. (b) One can start with an arbitrary sequence f′jεl2(X). (c) One can sample a function fεW0(L2) which is not necessarily in V2(φ).

[0357] In this situation, according to one embodiment of the present invention, we wish to know what happens if we run the frame algorithm with the input sequence {f′j:jεJ}. If {f′jjεJ}εl2(X), we can still initialize the iterative frame algorithm by 62$\begin{array}{cc}{g}_1=\sum _{j\in J}{f}_j^{\prime }{K}_{x_j}& \left(7.2\right)\end{array}$

[0358] This corresponds exactly to the first step in the iterative frame algorithm (5.3). Then we set 63$\begin{array}{cc}{g}_n=\frac{2}{A+B}{g}_1+\left(I-\frac{2}{A+B}T\right)g_{n-1}& \left(7.3\right)\end{array}$

[0359] Since {Kxj} is a frame for V2(φ) by assumption, this iterative algorithm still converges in L2, and its limit is 64$\begin{array}{cc}{g}_{\infty }=\underset{n\to \infty }{\lim }{g}_n=\sum _{j\in J}{f}_j^{\prime }{\tilde{K}}_{x_j}& \left(7.4\right)\end{array}$

[0360] THEOREM 7.2. Let X be a set of sampling for V2(φ). Then for any {f′j:jεJ}εl2(X), the modified frame algorithm (7.3) with the initialization (7.2) converges to 65$g_{\infty }=\sum _{j\in Z^d}{f}_j^{\prime }{\tilde{K}}_{x_j}\in V^2\left(\phi \right).$

[0361] We have that 66$\sum _{j\in J}\&LeftBracketingBar;f_j^{\prime }-{g_{\infty }\left(x_j\right)}^2\&RightBracketingBar;<\sum _{j\in J}\&LeftBracketingBar;f_j^{\prime }-{g\left(x_j\right)}^2\&RightBracketingBar;$

[0362] For all gεV2(φ) with equality if and only if g=g∞. Thus g∞ fits the given data optimally in the least squares sense.

[0363] Next we investigate the iterative algorithm (6.2) in the case of noisy samples {f′j:jεJ}εlpv(X). We use initialization 67$\begin{array}{cc}{f}_1=P{Q}_X{f}^{\prime }=P\left(\sum _{j\in J}{f}_j^{\prime }{\beta }_j\right)& \left(7.5\right)\end{array}$

[0364] and define the recursion as in (6.2) by

f n =f 1 +(I−PQX)fn−1  (7.6)

[0365] The convergence of this algorithm is clarified in the following theorem (see FIG. 8).

[0366] THEOREM 7.3. Under the same assumptions as in Theorem 6.1, the algorithm (7.6) converges to a function f∞εVpv(φ), which satisfies PQXf∞=PQX{f′j}.

[0367] A reconstruction of a function f with additive noise using the iteration algorithm as set forth in FIG. 5 and further disclosed above according to one embodiment of the present invention is shown in FIG. 8. Referring now to FIG. 8(a), curve 801 represents a function f belonging to the shift-invariant space generated by the Gaussian function exp(−x2/σ2) with σ≈0.81. As shown in FIG. 8(b), curve 811 represents the function f with an additive white noise, where SNR≈0 db). As shown in FIG. 8(c), curve or distribution 821 represents noisy signal sampled on a nonuniform grid with maximal gap 0.51. As shown in FIG. 8(d), curve 831 represents reconstructed function f10 only after 10 iterations (continuous line) in comparison (and agreement) with original signal f (dotted line) from curve 801.

EXAMPLE 2

Fast Local Reconstruction Methods for Uniform and Nonuniform Data Sampling in Shift Invariant Spaces

[0368] Approach Overview. In one aspect, the present invention relates to a method for fast reconstruction of a function f from its samples f(xj), where f belongs to a shift invariant space V(φ) and φ is a generator. If the generator φ has compact support, then the reconstruction is local quite in contrast to methods based on band-limited functions. In one embodiment, it is shown that the matrix of the corresponding linear system of equations is a positive-definite banded matrix. This special structure makes possible the fast local reconstruction algorithm in O(S2J) operations, where J is the number of samples and S is the support length of the generator φ. Further optimization can be achieved by means of data segmentation. Several numerical simulations according to the present invention are provided.

[0369] In another aspect, the present invention relates to a method for converting data for an object of interest, which is characterized by a functions, between a digital form and an analog form. In one embodiment, the method includes the steps of locally selecting a plurality of data sample points in the order of N, N being an integer, wherein the data sample points are in a first form of data type, performing a transformation in a shift-invariant space to the locally selected data sample points to obtain a new data set that is in a second form of data type different from the first form of data type, and reconstructing the object of interest f from the new data set. In one embodiment, the first form of data type is one of the digital form and the analog form, and the second form of data type is one of the digital form and the analog form that is different from the first form of data type. In other words, the first form of data type can be in digital form, and the corresponding second form of data type is in analog form, and vice versa.

1. Introduction

[0370] As set forth above, shift-invariant spaces serve as a universal model for uniform and nonuniform sampling of functions. Normally, the objective of the so-called sampling problem is either to recover a signal (function) f from its samples {f(xj) Z} or to approximate a data set (xj,yj) by a suitable function f satisfying f(xj)≠Yj. Obviously this problem is ill-posed, and so a successful reconstruction requires some a priori information about the signal. Usually it is assumed that f is contained in the span of integer translates of a given generator φ. In technical terms, the original function f has the form 68$f\left(x\right)=\sum _{k\in Z^{c_k}}\phi \left(x-k\right)$

[0371] and belongs to the shift-invariant space V(φ).

[0372] Until recently the only choice for p was the cardinal sine function 69$\phi \left(x\right)=\frac{\sin \pi \alpha x}{\pi \alpha x},$

[0373] since in this case V(φ) coincides with the band-limited functions of band-width 2α. Despite the existence of fast numerical methods [49], this model has some drawbacks because it is non-local, and the behavior of f at a point x also depends on samples far away from x. For this reason, one works with truncated versions of the cardinal sine. This idea leads naturally to work in shift-invariant spaces with a generator p of compact support. The concept of shift-invariant spaces first arose in approximation theory and wavelet theory [33,34,68]. Its potential for the systematic treatment of sampling problems was recognized much later.

[0374] Among other things, one goal of this Example is the investigation of specific numerical issues of nonuniform sampling in shift-invariant spaces. Usually the reconstruction algorithms are based on general frame methods or they follow simple iterative schemes. Of course, these can be applied successfully in the context of shift-invariant spaces as well as set forth above. We adopt the philosophy that the optimal numerical algorithms always have to use the special structure of the problem. Hence general purpose algorithms have to be fine-tuned to achieve their optimal performance. Here we use the peculiar structure of shift-invariant spaces with compactly supported generator to solve the sampling problem.

[0375] Among other things, we want to reconstruct a function fεV(φ) from a finite number of samples f(xj) taken from an interval xjε[M0, M1]. In our derivation of the special structure we use a frame theoretic argument and combine it with the fact that the generator φ has compact support. The resulting algorithm is local in the sense that the complete reconstruction of a function in V(φ) on the interval [M0, M1] only requires samples from [M0,M1] (quite in contrast to bandlimited functions). By comparing reconstructions in different spline spaces we find that, among other things, the algorithm can be used as an efficient denoising procedure for noisy samples.

[0376] The rest of the disclosure for Example 2 is organized as follows: In Section 2 we present the precise technical details of shift-invariant spaces, state the well-known equivalence of sampling problems with the construction of certain frames, and discuss some general reconstruction techniques associated to frames. In Section 3 we exploit the special structure of shift-invariant spaces to derive a local reconstruction algorithm of order O(J) and discuss the numerical issues involved. Section 4 explains the results of the numerical simulations and provides, as an example, a code for the main algorithm.

2. Relevant Concepts

[0377] 2.1. Shift-Invariant Spaces, Frames, and Sampling. Let φ be a continuous function with compact support of size S so that

supp φ[−S,S].  (1)

[0378] For convenience we assume that S is a positive integer and thus φ(±S)=0. Then the shift-invariant space V(φ) is defined as

V (φ)={fεL2():f=Σckφ(x−k) for (ck)εl2(Z)}.  (2)

[0379] To guarantee the stability of these representations, we assume that the generator φ is stable, which means that there exists a constant C0 such that (3) 70$C^{-1}{\&LeftDoubleBracketingBar;c\&RightDoubleBracketingBar;}_{l^2}\le {\&LeftDoubleBracketingBar;\sum _{k\in Z}{c}_k\phi \left(·-k\right)\&RightDoubleBracketingBar;}_2\le C{\&LeftDoubleBracketingBar;c\&RightDoubleBracketingBar;}_{l^2}$

[0380] for all finite sequences c=(ck)kεZ, or equivalently, the translates φ(.−k),kεZ, form a Riesz basis for V(φ). As a consequence, V(φ) is a closed subspace of L2(R) and inherits the inner product .,. of L2(R).

[0381] The sampling problem in V(φ) may be interpreted as the following problem: Given a set of sampling points X={xj:jεZ} arranged in increasing order xj<xj+1 and a sequence of samples {f(xj):jεZ}, of a function fεV(φ), we would like to recover the original function f in a stable and numerically efficient way. Here stability means that there exist (possibly unspecified) constants A, B>0 such that 71$\begin{array}{cc}A{\&LeftDoubleBracketingBar;f\&RightDoubleBracketingBar;}_2\le {\left(\sum _{j\in Z}{\&LeftBracketingBar;f\left(x_j\right)\&RightBracketingBar;}^2\right)}^{1/2}\le B{\&LeftDoubleBracketingBar;f\&RightDoubleBracketingBar;}_2\forall f\in V\left(\phi \right).& \left(4\right)\end{array}$

[0382] A sampling set satisfying (4) is called a set of stable sampling.

[0383] Thus, for (4) to be valid, we need point evaluations in V(φ) to be well-defined. This is guaranteed by the following lemma given in [5].

[0384] Lemma 1. If φ is continuous and satisfies the condition ΣkεZ maxxε[0,1]|f(x+k)|<∞, in particular, if φ is continuous with compact support, then for all xεR there exists a function KxεV(φ) such that f(x)=f,Kx for fεV(φ). We say that V(φ) is a reproducing kernel Hilbert space.

[0385] Explicit formulas for Kx, are known, see [5], but we do not need them here. Note that with the kernels Kx, the sampling inequality (4) can be formulated as 72$A{\&LeftDoubleBracketingBar;f\&RightDoubleBracketingBar;}_2\le {\left(\sum _{j\in Z}{\&LeftBracketingBar;⟨f,k_{x_j}⟩\&RightBracketingBar;}^2\right)}^{1/2}\le B{\&LeftDoubleBracketingBar;f\&RightDoubleBracketingBar;}_2,$

[0386] which is equivalent to saying that the set {Kxj:jεZ} is a frame for V(φ).

[0387] Let U be the infinite matrix with entries

U jk =φ(xj−k) j,kεZ.  (5)

[0388] Then the sampling problem in V(φ) can be formulated in several distinct ways [120].

[0389] Lemma 2. If φ satisfies the condition of Lemma 1, then the following are equivalent:

[0390] (i) X={xj:jεZ} is a set of sampling for V(φ).

[0391] (ii) There exist A, B>0 such that

A∥c∥ l 2 ≦∥Uc∥ l 2 ≦B∥c∥ l 2 ∀cεl 2 (Z)

[0392] (iii) The set of reproducing kernels {Kxj:jεZ} is a frame for V(φ).

[0393] REMARK. It is difficult to characterize sets of sampling for V(φ). If φ is a B-spline of order N, i.e., φ=X[0,1]* . . . * X[0,1](N+1 convolutions), then the main result of [120] implies that the maximum gap condition supjεZ(xj+1−xj)=δ<1 is sufficient for the conditions of Lemma 2 to hold.

[0394] 2.2. General Purpose Reconstructions. Lemma 2 leads to some general reconstruction techniques that are always applicable.

[0395] A. Linear Algebra Solution. One could simply try to solve the (infinite) system of linear equations 73$\begin{array}{cc}\sum _{k\in Z}{c}_k\phi \left(x_j-k\right)=f\left(x_j\right)\forall j\in Z,& \left(6\right)\end{array}$

[0396] for the coefficients (ck), or in the notation of (5) with f|x=(f(xj))jεZ

UC=f|x   (7)

[0397] B. The Normal Equations. Frequently it is better to consider the associated system of normal equations [124]

U*Uc=U*f|x   (8)

[0398] This approach has the advantage that the matrix T:=U*U is a positive operator on l2(z). Furthermore, if the input y=(yj)jεZ does not consist of a sequence of exact samples of fεV(φ), then the function 74$f=\sum _{k\in Z}{c}_k\phi \left(·-k\right)$

[0399] corresponding to the solution c=(U*U)−1U*y solves tile least squares problem 75$\begin{array}{cc}\sum _{j\in Z}{\&LeftBracketingBar;y_j-f\left(x_j\right)\&RightBracketingBar;}^2=\underset{h\in v\left(\phi \right)}{\min }\sum _{j\in Z}{\&LeftBracketingBar;y_j-h\left(x_j\right)\&RightBracketingBar;}^2.& \left(9\right)\end{array}$

[0400] C. Frame Approach. Lemma 2(iii) suggests to use versions of the frame algorithm to find a reconstruction of f. By (iii) the frame operator which is defined as 76$\begin{array}{cc}Sf\left(x\right)=\sum _{j\in Z}⟨f,K_{x_j}⟩K_{x_j}\left(x\right)=\sum _{j\in Z}f\left(x_j\right)K_{x_j}\left(x\right)& \left(10\right)\end{array}$

[0401] is invertible and its inverse defines the dual frame xj=S−1KxjεZ. Then the reconstruction is given by 77$\begin{array}{cc}f\left(x\right)=\sum _{j\in J}⟨f,K_{x_j}⟩•K_{x_j}=\sum _{j\in J}f\left(x_j\right)•K_{x_j}\left(x\right).& \left(11\right)\end{array}$

[0402] We observe that the linear algebra solution (7) and the frame method are equivalent. By definition the vector of samples is given by Uc=f|x fix. The sampled energy of fεV(φ) is 78$\begin{array}{cc}\sum _{j\in Z}{\&LeftBracketingBar;f\left(x_j\right)\&RightBracketingBar;}^2={⟨f\&LeftBracketingBar;x,f\&RightBracketingBar;x⟩}_{l^2}={⟨\mathrm{Uc},\mathrm{Uc}⟩}_{l^2}={⟨U^{*}\mathrm{Uc},c⟩}_{l^2}& \left(12\right)\end{array}$

[0403] Thus X is a set of sampling if and only if U*U is invertible on l2(Z)

[0404] D. Iterative Frame Methods. In nonuniform sampling problems it is usually difficult to calculate the entire dual frame, therefore one often resorts to iterative methods. Since the Richardson-Landweber iteration in the original paper of Duffin-chaeffer [40] is slow and requires good estimates of the frame bounds, we recommend the conjugate gradient acceleration of the frame algorithm for all problems without additional structure [56]. It converges optimally and does not require the estimate of auxiliary parameters.

[0405] 3. Exploit the Structure of the Problem

[0406] So far we have discussed general purpose methods for the reconstruction of the function. These could be applied in any situation involving frames and do riot take into consideration the particular structure of the sampling problem in shift-invariant spaces.

[0407] 3.1. A Localization Property. We now exploit the special structure of shift-invariant spaces. The following lemma is simple, but crucial. It is a consequence of the assumption that the generator of V(φ) has compact support.

[0408] Lemma 3. If supp φ[−S, S], then T=U*U is a band matrix of (upper and lower) band-width 2S.

[0409] Proof. By definition the entries of U*U are 79${\left(U^{*}U\right)}_{kl}=\sum _{j\in Z}\stackrel{\_}{U_{jk}}{U}_{jl}=\sum _{j\in Z}\stackrel{\_}{\phi \left(x_j-k\right)}\phi \left(x_j-l\right).$

[0410] Since φ has compact support, the sum is always locally finite and its convergence does not pose any problem. Since φ(xj−k)=0 if |xj−k|≧S, we find that (U*U)kl can be non-zero only if both |xj−k|<S and |xj−l|<S. In other words, (U*U)kl≠0 implies that

|k−l|≦|k−xj|+|xj−l|<2S.

[0411] This means that only 4S−1 diagonals of U*U contain non-zero entries.

[0412] REMARKS. A. Banded matrices and the resulting numerical advantages occur in a number of related problems. For instance, in the interpolation of scattered data by radial functions with compact support, the interpolation matrix is banded, see [121] and the references there. Likewise, the calculation of the optimal smoothing spline on a finite set of arbitrary nodes requires the inversion of a banded matrix [125].

[0413] B. Lemma 3 combined with a result of Demko, Moss, and Smith [123] or of Jaffard [126] implies that the inverse matrix possesses exponential decay off the diagonal, i.e., there exist C, A>0 such that 80$\&LeftBracketingBar;{\left(U^{*}U\right)}_{kl}^{-1}\&RightBracketingBar;\le {\mathrm{Ce}}^{-A\&LeftBracketingBar;k-l\&RightBracketingBar;}\forall k,l\in Z.$

[0414] To make our treatment more realistic, we take into account that in any real problem only a finite (albeit large) number of samples is given. It turns out that the model of shift-invariant spaces with compactly supported generator possesses excellent localization properties. These are quantified in the next lemma.

[0415] Lemma 4. The restriction of fεV(φ) to the interval [M0,M1] is determined completely by the coefficients ck for kε(M0−S, M1+S)∩Z.

[0416] Proof. Since φ(x−k)=0 for |x−k|≧S and xε[M0,M1], we obtain that

M 0 −S−x−S<k<x+S≦M 1 +S.

[0417] Consequently, as Sε, we have 81$\begin{array}{c}f\left(x\right)=\sum _{k\in Z}{c}_k\phi \left(x-k\right)\\ =\sum _{\&LeftBracketingBar;x-k\&RightBracketingBar;<S}{c}_k\phi \left(x-k\right)\\ =\sum _{k=M_0-S+1}^{M_1+S-1}{c}_k\phi \left(x-k\right).\end{array}$

[0418] In other words, the exact reconstruction of fεV(φ) on [M0,M1] requires only the M1−M0+2S−1 unknown coefficients ck with kε(M0−S,M1+S)∩Z. By counting dimensions, we find that we need at least M1−M0+2S−1 samples in [M0,M1] for the coefficients to be determined uniquely. Usually the length M1−M0 is large compared to S, therefore the additional 2S−1 coefficients amount to a negligible oversampling.

[0419] Lemma 4 demonstrates an important theoretical and practical advantage of shift-invariant spaces with compactly supported generators. A function fεV(φ) can be reconstructed exactly on an arbitrary interval solely from samples in that interval. In contrast, the restriction of a band-limited function to an interval is not uniquely determined by any finite number of samples in that interval, but can only be approximated by these samples. The localization property expressed in Lemma 4 is one of the main reasons to work with shift-invariant spaces with compactly supported generators as a sampling model!

[0420] Finally we remark that uniform sampling at critical density is not local and may even be unstable in this model. If fεV(φ) is sampled at ξ+k, kεZ, for some ξε[0,1), then there exists an interpolating function ψξ of exponential decay, such that 82$f\left(x\right)=\sum _{k\in Z}f\left(\xi +k\right){\psi }_{\xi }\left(x-k\right)\left[64\right].$

[0421] In this case the restriction of f to [M0,M1] is not determined exclusively by the M1−M0 values f(ξ+k) for ξ+kε[M0, M1]. Moreover, if φ is continuous, then there always exists a ξε[0,1) such that the reconstruction {f(ξ+k)}→X f is unstable. Jaussen's results in [64] indicate that a small amount of oversampling is an essential hypothesis to guarantee the locality and the stability of its reconstruction.

[0422] 3.2. A Local Reconstruction Algorithm. In practice we perform the calculations with a truncated version of the matrices U and T. We now combine Lemmas 3 and 4 to provide an efficient numerical reconstruction algorithm according to one embodiment of the present invention.

[0423] Algorithm

[0424] Input. We assume that finitely many sampling points x1, . . . , xJε[M0,M1] are given with associated sampling vector y=(y1, . . . , yJ)εJ. Assume that J≧M1−M0+2S−1 and that the truncated matrix z defined below is invertible.

[0425] Step 0. First we define and compute the truncated matrices U=UM0,M1 and T=TM0,M1=U*U, given by their entries 83$\begin{array}{cc}{U}_{jk}=\phi \left(x_j-k\right)T_{kl}=\sum _{j=1}^J\stackrel{\_}{\phi \left(x_j-k\right)}\phi \left(x_j-1\right)& \left(13\right)\end{array}$

[0426] for j=1, . . . , J and k,l=M0−S+1, . . . , M+S−1.

[0427] Step 1. Compute b=U*y, i.e., 84$\begin{array}{cc}{b}_k=\sum _{j=1}^J\stackrel{\_}{\phi \left(x_j-k\right)}{y}_j\mathrm{for}k=M_0-S+1,\dots ,M_1+S-1& \left(14\right)\end{array}$

[0428] Step 2. Solve the system of equations

c=T −1 b.   (15)

[0429] Step 3. Compute the restriction of f to [M0,M1] by 85$\begin{array}{cc}f\left(x\right)=\sum _{k=M_0-S+1}^{M_1+S-1}{c}_k\phi \left(x-k\right)\mathrm{for}x\in \left[M_0,M_1\right].& \left(16\right)\end{array}$

[0430] Then f is the (unique) least square approximation of the given data vector y in the sense that 86$\begin{array}{cc}\sum _{j=1}^J{\&LeftBracketingBar;y_j-f\left(x_j\right)\&RightBracketingBar;}^2=\underset{h\in V\left(\phi \right)}{\min }\sum _{j=1}^J{\&LeftBracketingBar;y_j-h\left(x_j\right)\&RightBracketingBar;}^2.& \left(17\right)\end{array}$

[0431] If y arises as the sampled vector of an f c V(φ), i.e., yj=f (xj), then this algorithm provides the exact reconstruction of f

[0432] Proof. The least square property (17) is clear, since this is exactly the property of the solution of the system of normal equations U*U c=U*y. See [124] for details.

[0433] In the case of B-splines a sufficient condition on the sampling density can be extracted from the proofs of Thin. 2.1 and 2.2 of [120]. Assume that xj+1−xj≦δ, and that 87$\begin{array}{cc}\delta \le \frac{M_1-M_0}{M_1-M_{0}2S-1}<1.& \left(18\right)\end{array}$

[0434] Then f is invertible. Condition (18) guarantees that there are at least M1−M0+2S−1 samples in[M0,M1]. Then the Schoenberg-Whitney theorem [127, p. 167] implies that T is invertible. See [120] for the detailed arguments.

[0435] 3.3. Data Segmentation. In another embodiment of the present invention, a further optimization of the reconstruction procedure is possible by data segmentation. Instead of solving the large system of equations

T m 0 ,m 1 =(UM0,M1)*y,

[0436] with a band matrix of dimension M1−M0+2S−1, we will solve t systems of smaller size. For this purpose we partition the large interval [M0, M1] into t smaller intervals [mr,mr+1], r=0, . . . , t−1 with M0=m0 and M1=mt.

[0437] Now we apply Algorithm 3.2 set forth above to each interval separately. More precisely, given the data (xj, yj) where xjε[mr, mr+1], we set tip the matrices Umr,mr+1 and Tmr,mr+1 and solve t equations

T m r ,m r+1 c (r) =(Umr,mr+1)*y(r)  (19)

[0438] where the vector y(r) consists of those data yj for which xjε[mr,mr+1] and the coefficient vector c(r)=(cmr−S+1, . . . ,Cmr+1+S−1)

[0439] Among other things, the segmentation technique has a number of practical advantages as follows:

[0440] A. The dimension of vectors and matrices can be reduced drastically. Using data segmentation, we solve t small systems of size (M−M0)/t+2S−1 instead of the large system of size M1−M0+2S−1.

[0441] B. Parallel processing can be applied because non-adjacent intervals can be handled simultaneously.

[0442] C. The function can be reconstructed on specified subintervals at smaller cost. See discussions associated with FIG. 6, infra.

[0443] On the other hand, data segmentation may also come with some caveats:

[0444] A. The coefficients ck with indices kε[mr−S+1, mr+1+S−1] are computed at least twice because of overlap. Heuristically it has proved best to take averages of the multiply computed coefficients.

[0445] B. For a successful execution of the segmentation method it is necessary that each of the small matrices Tmr,mr+1 is invertible. Again by dimension counts we find that the number of data in the interval [mrmr+1] should exceed the number of variables,

#(X∩[mrmr+1])≧mr+1−mr+2S−1

[0446] Obviously this condition imposes an upper hound for the possible number of segmentations.

[0447] 3.4. Implementation Issues. A. In Algorithm 3.2 as set forth above, the most expensive step is the calculation of the matrix U because it requires the point evaluations of φ. However, if the sampling points xj are given, then U and T can be computed in advance and stored. Thus Step 0 can be taken care of before solving the reconstruction problem.

[0448] We handle the pointwise evaluation of φ by “quantizing” the generator. This means that for δ>0 sufficiently small we create of vector ψ consisting of entries 88$\phi \left(\frac{1}{N}\right)\mathrm{for}l=-\mathrm{NS},\dots ,\mathrm{NS}$

[0449] such that 89$\begin{array}{ccc}\&LeftBracketingBar;\phi \left(x\right)-\phi \left(\frac{l}{N}\right)\&RightBracketingBar;<\delta & & \mathrm{for}\&LeftBracketingBar;x-\frac{1}{N}\&RightBracketingBar;<\frac{1}{2N}\end{array}$

[0450] Thus to build the matrix Ujk=φ(xj−k) amount to selecting the appropriate entries of ψ. This approximation of U works remarkably well and fast in the numerical simulations.

[0451] B. For the solution of the banded system (15) a number of fast algorithms are available. Golub-van Loan [124] offers several efficient algorithms for this task; other options for the inversion of a banded matrix are mentioned in [125]. Since T is assumed to be positive definite, the band Cholesky algorithm seems to be a good choice that minimizes the operation count for Step 2. In one embodiment, MATLAB© is utilized to practice the present invention and it provides the commands SPARSE and CHOL to deal with this task.

[0452] C. Usually f is reconstructed on a grid 90$G=\left\{\frac{1}{N}:l=M_{0}N,\dots ,M_{1}N\right\}.$

[0453] Then (16) amounts to a discrete convolution, and thus Step 3 can be performed quickly. Again, since φ has compact support, we can use the banded structure of the associated matrix to perform this step.

[0454] 3.5. Operation Count. We estimate the number of multiplications for Algorithm 3.2. Recall that J is the number of samples, and D=M1−M0+2S−1 is the dimension of the problem.

[0455] (a) According to (14) each of the D entries of the vector b requires #{j:|xj−k|<S} multiplications. Consequently Step 1 requires 91$\begin{array}{c}\sum _{k=M_0-S+1}^{M_1+S-1}\#\left\{j:\&LeftBracketingBar;x_j-k\&RightBracketingBar;<S\right\}=\sum _{k=M_0-S+1}^{M_1+S-1}\sum _{j=1}^{J}X\left(k-S,k+S\right)\left(X_j\right)\\ =\sum _{j=1}^J\sum _{k=M_0-S+1}^{M_1+S-1}X\left(k-S,k+S\right)\left(X_j\right)\\ =\sum _{j=1}^{J}2S=2SJ\end{array}$

[0456] operations, because a point x is in at most 2S translates of the open interval (−S S).

[0457] (b) Likewise to calculate an entry of T requires #({j:|xj−k|<S}∩{j:|xj−l|<S}) multiplications, see (13). As in (a) we estimate the number of operations to set up the matrix T by 92$\begin{array}{c}\sum _{k=M_0-S+1}^{M_1+S-1}\sum _{l=M_0-S+1}^{M_1+S-1}\#\left(\left\{j:\&LeftBracketingBar;x_j-k\&RightBracketingBar;<S\right\}\bigcap \left\{j:\&LeftBracketingBar;x_j-l\&RightBracketingBar;<S\right\}\right)=\\ \sum _{k=M_0-S+1}^{M_1+S-1}\sum _{l=M_0-S+1}^{M_1+S-1}\sum _{j=1}^{J}X\left(k-S,k+S\right)\left(x_j\right)X\left(l-S,l+S\right)\left(x_j\right)=\\ \sum _{j=1}^J\left(\sum _{k=M_0-S+1}^{M_1+S-1}X\left(k-S,k+S\right)\left(x_j\right)\right)\left(\sum _{l=M_0-S+1}^{M+S-1}X\left(l-S,l+S\right)\left(x_j\right)\right)\end{array}$

[0458] <J·(2S)2.

[0459] (c) For the solution of the banded system T c=b by means of the band Cholesky algorithm we need at most

D ((2S)2+16S+1)=(M1−M0+2S−1)((2S)2+16S+1)≦J(4S2+16S+1)

[0460] operations (and no square roots), see [124, Ch. 4.3.6].

[0461] (d) To compute the reconstruction f on a grid 93$\left\{\frac{l}{N}:l=M_{0}N,\dots ,M_{1}N\right\}$

[0462] we need to calculate (M1−M0) N point evaluations of f via (16). Since, each point evaluation requires at most 2S multiplications. Thus for the reconstruction on the grid we need at most

(M1−M0)N·2S≦J·2SN

[0463] multiplications.

[0464] Combining these estimates, we find that Algorithm 3.2 requires

O(J(S2+SN))  (20)

[0465] operations. In other words, the cost of the algorithm according to this embodiment of the present invention is linear in the number of data and quadratic in the size of the generator.

[0466] 4. Numerical Simulations

[0467] In one embodiment, several simulations performed according to the methods and algorithm(s) of the present invention. In our simulation we have used MATLAB. We used the shift-invariant spline spaces with the B-spline of order 3 94$\phi =\underset{\underset{4\mathrm{times}}{⏟}}{{X\left[-1/2,1/2\right]}^{*}\dots {\hspace{0.17em}}^{*}X\left[-1/2,1/2\right]}$

[0468] as the generator of V(φ). Thus supp φ[−2,2] and S=2. Referring now to FIG. 10, the shift-invariant spline spaces with the B-spline of order 3 with φ[−2,2] and S=2 that can be utilized to practice the present invention is shown. Curve 1101 in FIG. 11 is shown as a plot of the operation count as a function of the number of sampling points. For this simulation, we have reconstructed examples of size 114, 226, 444, 667, 887, 1085, identified by the circles 1103 in FIG. 11, respectively, and used the MATLAB function FLOPS to count the number of operations.

[0469] Referring now to FIG. 12, a reconstruction of a function f with additive noise according to one embodiment of the present invention with a signal on the interval [0,128] is shown, where M0=0, and M1=128. Since S=2, one needs at least M1−M0+2S−1=131 samples. The actual sampling set of FIG. 12 has 200 sample points as calculated and satisfies the maximum gap condition maxj(xj+1−xj)≈0.67<1. More specifically, curve 1201 as shown in FIG. 12(a) represents original signal.

[0470] To make the example more realistic, we have added white noise to the sampled values of a given function fεV(φ). Instead of missing the correct values f(xj) in the reconstruction algorithm, we use the noisy values ferr(xj)=f(xj)+ej so that

ferr|x=f|+e.

[0471] The relative error between the original signal and the noisy signal is measured by 95${\mathrm{err}}_{\mathrm{samp}}=\frac{{\&LeftDoubleBracketingBar;\mathrm{ferr}\&LeftBracketingBar;x-f\&RightBracketingBar;x\&RightDoubleBracketingBar;}_2}{{\&LeftDoubleBracketingBar;f❘x\&RightDoubleBracketingBar;}_2}={\left(\frac{\sum _{j=1}^J{\&LeftBracketingBar;e_j\&RightBracketingBar;}^2}{\sum _{j=1}^J{\&LeftBracketingBar;f\left(x_j\right)\&RightBracketingBar;}^2}\right)}^{1/2}.$

[0472] In our example errsamp=63.8%. Curve 1211 as shown in FIG. 12(b) represents original signal with noise.

[0473] Moreover, as shown in FIG. 12(c), curve 1221 represents noisy signal sampled on a nonuniform grid. Additionally, as shown in FIG. 12(d), curve 1231 represents reconstructed original signal after 10 iterations (continuous line) in comparison with original signal dotted line from curve 1201. The relative error errrec of the reconstruction measured at the sampling points with respect to the correct samples is now 96${\mathrm{err}}_{\mathrm{rec}}=\frac{{\&LeftDoubleBracketingBar;f_{\mathrm{rec}}\&LeftBracketingBar;x-f\&RightBracketingBar;x\&RightDoubleBracketingBar;}_2}{{\&LeftDoubleBracketingBar;f❘x\&RightDoubleBracketingBar;}_2}=18.5\%$

[0474] The noise reduction is thus

err samp =63.8%→errrec=18.5%

[0475] Referring now to FIGS. 13A and 13B, the dependence of the reconstruction on the generator φ is investigated and shown, where the generator φ is a B-spline of order N, i.e. φN=X[−1/2,1/2]* . . . *X[−1/2,1/2] (N+1-fold convolution). The data set (xj,yj)j<1, . . . ,J is generated by sampling a function fεV(φ5) and then noise is added. Curve 1301 and 1303 in FIG. 13A(a) represent the original signal and the noisy sampled data, respectively. Then the optimal approximation of these data in the spline space V(φN), N=0, . . . , 6, starting with an approximation by a step function frecεV(φ0) via an approximation by a piecewise linear function frecεV(φ1) and ending with a smooth approximation frecεV(φ6), are shown by curves 1311, 1321, 1331, 1341, 1351, 1361 and 1371 in FIG. 13A(b)-(d) and FIG. 13B(e)-(h), respectively.

[0476] In each case we have also plotted the original function f (dotted line curve 1301) for comparison. In addition, the relative error errsamp is indicated. The dependency of this error of N is a typical L-curve as it occurs in regularization procedures. In all our examples the best approximation is obtained in the correct space in which f was originally generated. This observation is consistent with the extended literature on smoothing splines, e.g., [122, 125, 128]. One of the main differences between those methods and the algorithm of Section 3.2 of the present invention is in the underlying function space. The reconstruction algorithm 3.2 finds the best local reconstruction in the shift-invariant spline spaceV(φ), whereas the smoothing spline of [122] is based on the nodes xj and does not belong to a fixed function space.

[0477] Referring now to FIG. 14, the associated banded matrix T of the linear system according to one embodiment of the present invention is shown as curve (substantially a linear line/lines) 1401. In FIG. 14, white squares correspond to zero entries, dark squares signify large entries of T, the shading being proportional to the size. The banded structure is clearly visible.

[0478] Referring now to FIG. 15, the power of the method of data segmentation according to one embodiment of the present invention is shown. Instead of reconstructing the entire signal f, we have reconstructed only the restriction to two disjoint intervals, represented by curves 1501, 1503 respectively, in comparison with original data curve 1301. In the absence of noise the reconstruction is exact. Since supp(φ)⊂[−S, S], the calculation for the two intervals can be done locally and simultaneously. This property can be used for parallel processing.

5. An Exemplary Code

[0479] By way of example, the following is a code for fast local data reconstruction in shift invariant spaces according to one embodiment of the present invention as set forth above:

[0480] function f_rec=reconstruction(xp, xs, x_rec, t);

function f_rec = reconstruction(xp, xs, x_rec, t);
% xp	... sampling positions
% xs	... sampling values
% x_rec	... positions, where the function should be reconstructed
% gen	... generator with support supp(gen)=[−S,S]
% t	... number of segmentation
step = round((max(x_rec)−min(x_rec))/t);
for k = min(x_rec):step:max(x_.rec)
	xp_rel = {xp: k−S < xp < k + step + S}	% relevant sampling positions
		% or the interval [k, k + step]
	xp_mm = min(xp.rel);
	xp_max = max(xprel);
	J = length(xp_rel);	% number of sampling points
% calculation of the left side b:
	for j=1:J
	mi=ceil (xp(j))−S;
	ma=floor (xp(j) )+S;
	for l = mi : ma
	b(l−mi+1) = b(l−mi+1) + xs(j)*gen(xp(j)−l);
	end
	end
% calculation of the matrix T:
	T=zeros(xp_max−xp_min+1+2*S,xp_max−xp_min+1+2*S);
	for j=1:J
	mi=ceil (xp(j))−S;
	ma=floor (xp(j) )+S;
	for k =mi : ma
	for l=ma : ma
	T(k−mi+1,l−mi+1) = T(k−mi+1,l−mi+1) + gen(xp(j)−
	l)*gen(xp(j)−l);
	end
	end
	end
%	calculation of the coefficients
	c_part = chol(T,b);	% solving the system T*c_part=b with a
		% banded Cholesky algorithm
	c(xp_min−S:xp_max+S) = c(xp_min−S:xp_max+S) + c_part;
	n(xp_min−S :xp_max+S) = n(xp_min−S:xp_max+S) + ones(xp_max−
xp_min+1+2*S);
	%n . . . normalization of coefficients because of overlapping
end
c = c ./ n;	%normalization of coefficients because of overlapping
% calculation of the reconstruction
for i = 1	length(x_rec)
	for j = floor(x_rec(i)−S)	:	ceil(x_rec(i)+S)	% |x_rec−j| <= S
	f_rec(x_rec(i)) = f_rec(x_rec(i)) + gen(x_rec(i)−j) * c(j);
	end
end

[0481] While this exemplary code is written in C language, it can be written in other languages such as Fortran, Pascal, C++, or the like so long as it can be executed by at least one computer such as Computer 8.

[0482] While there has been shown various embodiments of the present invention, it is to be understood that certain changes can be made in the form and arrangement of the elements of the system and steps of the methods to practice the present invention as would be known to one skilled in the art without departing from the underlying scope of the invention as is particularly set forth in the claims. Furthermore, the embodiments described above are only intended to illustrate the principles of the present invention and are not intended to limit the claims to the disclosed elements.

LIST OF REFERENCES

[0483] [1] A. ALDROUBI, Non-Uniform Weighted Average Sampling and Exact Reconstruction in Shift-Invariant Spaces, preprint, 2001.

[0484] [2] A. ALDROUBI, Portraits of frames, Proc. Amer. Math. Soc., 123 (1995), pp. 1661-1668.

[0485] [3] A. ALDROUBI, Oblique projections in atomic spaces, Proc. Amer. Math. Soc., 124 (1996), pp. 2051-2060.

[0486] [4] A. ALDROUBI AND H. FEIGHTINGER. Exact iterative reconstruction algorithm for multivariate irregularly sampled functions in spline-like spaces: The L˜theory, Proc. Amer. Math. Soc., 126 (1998), pp. 2677-2686.

[0487] [5] A. ALDROUBI AND K. GRÖCHENIG, Nonuniform sampling and reconstruction in shift variant space, SIAM Rev. 43 (2001), No. 4, pp. 585-620.

[0488] [6] A. ALDROUBI. Q. SUN AND W. S. TANG, p-frames and shift invariant spaces of LP, J. Fourier Anal. Appl., 7 (2001), pp. 1-21.

[0489] [7] A. ALDROUBI AND M. UNSER. Families of wavelet transforms in connection with Shannon's sampling theory and the Gabor transform, in Wavelets: A Tutorial in Theory and Applications, C. K. Chul. ed., Academic Press. San Diego, Calif., 1992, pp. 509-528.

[0490] [8] A. ALDROUBI AND M. UNSER, Families of multiresolution and wavelet spaces with optimal properties, Numer. Funct. Anal. Optim., 14 (1993), pp. 417-446.

[0491] [9] A. ALDROUBI AND M. UNSER, Sampling procedure in function spaces and asymptotic equivalence with Shannon's sampling theory, Numer. Funct. Anal. Optim., 15 (1994), pp. 1-21.

[0492] [10] A. ALDROUBI. M. UNSER AND M. EDEN, Asymptotic properties of least square spline filters and application to multi-scale decomposition of signals, in Proceedings of the International Conference on Information Theory and Its Applications, Waikiki, Hi., 1990, pp. 271-274.

[0493] [11] A. ALDROUBI. M. UNSER AND M. EDEN, Cardinal spline filters: Stability and convergence to the ideal sinc interpolator, Signal Processing, 28 (1992), pp. 127-138.

[0494] [12] R. BALAN, Equivalence relations and distances between Hilbert frames, Proc. Amer. Math. Soc., 127 (1999), pp. 2353-2366.

[0495] [13] E. BELLER AND G. DE HAAN, New algorithms for motion estimation on interlaced video, in Proc. SPIE-Visual Communication and Image Processing, 3309 (1998), pp. 111-121.

[0496] [14] J. J. BENEDETTO, Irregular sampling and frames, in Wavelets: A Tutorial in Theory and Applications. C. K. Chui, ed., Academic Press. San Diego, Calif., 1992, pp. 445-507.

[0497] [15] J. J. BENEDETTO, Frame decompositions, sampling, and uncertainty principle inequalities, in Wavelets-Mathematics and Applications, J. J. Benedetto and M. W. Frazier, eds., CRC, Boca Raton, Fla., 1993, pp. 247-304.

[0498] [16] J. J. BENEDETTO AND P. J. S. G. FERREIRA, Modern Sampling Theory, Birkhäuser, Boston, 2000.

[0499] [17] J. J. BENEDETTO AND M. W. FRAZIER. ED., Wavelets-Mathematics and Applications, CRC, Boca Raton, Fla., 1993.

[0500] [18] J. J. BENEDETTO, C. HEIL. AND D. F. WALNUT, Gabor systems and the Ealian-Low theorem, in Gabor Analysis and Algorithms, H. G. Feichtinger and T. Strohmer, eds., Birkhäuser, Boston, 1998, pp. 85-122.

[0501] [19] J. J. BENEDETTO AND S. Li, The theory of multiresolution frames and applications to filter banks, Appl. Comput. Harmon. Anal., 5 (1998), pp. 389-427.

[0502] [20] J. J. BENEDETTO AND D. WALNUT, Gabor frames of L2 and related spaces, in Wavelets-Mathematics and Applications, J. J. Benedetto and M W. Frazier, eds., CRC, Boca Raton, Fla., 1993, pp. 247-304.

[0503] [21] J. J. BENEDETTO AND H.-C. Wu, Non-uniform sampling and spiral MRI reconstruction, in Proc. SPIE-Wavelet Applications in Signal and Image Processing VIII, 4119 (2000), pp. 130-141.

[0504] [22] C. A. BERNSTEIN AND E. V. PATRICK, Exact deconvolution for multiple convolution operators-an overview, plus performance characterizations for imaging sensors, in Proceedings in Multidimensional Signal Processing, IEEE, Piscataway, N.J., 1990, pp. 723-734.

[0505] [23] A. BEURLING, The Collected Works of Arne Beurling. Vol. 1, Complex Analysis, L. Carleson, P. Malliavin, J. Neuberger, and J. Wermer, eds., Birkhäuser, Boston, 1989.

[0506] [24] A. BEURLING, The Collected Works of Arne Beurling. Vol. 2, Harmonic Analysis, L. Carleson, P. Malliavin, J. Neuberger, and J. Wermer, eds., Birkhäuser, Boston, 1989.

[0507] [25] A. BEURLING AND P. MALLIAVIN, On the closure of characters and the zeros of entire functions, Acta Math., 118 (1967), pp. 79-93.

[0508] [26] T. BLU AND M. UNSER, Quantitative Fourier analysis of approximation techniques: Part I—interpolators and projectors, IEEE Trans. Signal Process., 47 (1999), pp. 2783-2795.

[0509] [27] P. L. BUTZER, A survey of the Whittaker-Shannon sampling theorem and some of its extensions, J. Math. Res. Exposition, 3 (1983), pp. 185-212.

[0510] [28] P. G. CASAZZA AND O. CHRISTENSEN, Frames containing a Riesx basis and preservation of this property under perturbations, SIAM J. Math. Anal., 29 (1998), pp. 266-278.

[0511] [29] P. CASAZZA, D. RAN AND D. LARSON, Frames for Banach spaces, in The Functional and Harmonic Analysis of Wavelets and Frames, Contemp. Math. 247, AMS, Providence, R.I., 1999, pp. 149-182.

[0512] [30] W. CHEN, S. ITOH AND J. SHIKI, Irregular sampling theorems for wavelet subspaces, IEEE Trans. Inform. Theory, 44 (1998), pp. 1131-1142.

[0513] [31] O. CHRISTENSEN, Moment problems for frames and applications to irregular sampling and Gabor frames, Appl. Comput. Harmon. Anal., 3 (1996), pp. 82-86.

[0514] [32] C. K. CHUI, ED., Wavelets: A Tutorial in Theory and Applications, Academic Press, San Diego, Calif., 1992.

[0515] [33] I. DAUBEGHIES, Ten Lectures on Wavelets, SIAM, Philadelphia, 1992.

[0516] [34] C. DE BOOR, R. DEVORE AND A. RON, The structure of finitely generated shift-invariant spaces in L2(Rd), J. Funct. Anal., 119 (1994), pp. 37-78.

[0517] [35] C. DE BOOR AND R. DEVORE, Partition of unity and approximation, Proc. Amer. Math. Soc., 93 (1985), pp. 705-709.

[0518] [36] R. A. DEVORE, B. JAWERTH AND B. J. LUCIER, Image compression through wavelet transform coding, IEEE Trans. Inform. Theory, 38 (1992), pp. 719-746.

[0519] [37] I. DJOKOVIC AND P. P. VAIDYANATHAN, Generalized sampling theorems in multiresolution subspaces, IEEE Trans. Signal Process., 45 (1997), pp. 583-599.

[0520] [38] Y. DOMAR, Harmonic analysis based on certain commutative Banach algebras, Acta Math., 96 (1956), pp. 1-66.

[0521] [39] C. C. DONOVAN, J. S. GERONIMO AND D. P. HARDIN, Intertwining multiresolution analyses and the construction of piecewise-polynomial wavelets, SIAM J. Math. Anal., 27 (1996), pp. 1791-1815.

[0522] [40] R. J. DUFFIN AND A. C. SCHAEFFER, A class of nonharmonic Fourier series, Trans. Amer. Math. Soc., 72 (1952), pp. 341-366.

[0523] [41] H. G. FEICHTINGER, Gewichtsfunktionen auf lokalkompakten Gruppen, Sitzber. d. Österr. Akad. Wiss., 188 (1979), pp. 451-471.

[0524] [42] H. G. FEICHTINGER, Banach convolution algebras of Wiener type, in Functions, Series, Operators, Vols. I, H (Budapest, 1980), North-Holland, Amsterdam, 1983, pp. 509-524.

[0525] [43] H. G. FEICHTINGER, Generalized amalgams, with applications to Fourier transform, Canad. J. Math., 42 (1990), pp. 395-409.

[0526] [44] H. G. FEICHTINGER, New results on regular and irregular sampling based on Wiener amalgams, in Proc. Conf. Function Spaces, K. Jarosz, ed., Lecture Notes in Math. 136, Springer-Verlag, New York, 1991, pp. 107-121.

[0527] [45] H. G. FEICHTINGER, Wiener amalgams over Euclidean spaces and some of their applications, in Proc. Conf. Function Spaces. K. Jarosz, ed., Lecture Notes in Math. 136, Springer-Verlag, New York, 1991, pp. 123-137.

[0528] [46] H. G. FEICHTINGER AND K. GRÖCHENIG, Banach spaces related to integrable group representations and their atomic decompositions, I, J. Funct. Anal., 86 (1989), pp. 307-340.

[0529] [47] H. G. FEICHTINGER AND K. GRÖCHENIG, Iterative reconstruction of multivariate band-limited functions from irregular sampling values, SIAM J. Math. Anal., 23 (1992), pp. 244-261.

[0530] [48] H. G. FEICHTINGER AND K. GRÖCHENIG, Theory and practice of irregular sampling, in Wavelets-Mathematics and Applications. J. J. Benedetto and W. Frazier, eds., CRC, Boca Raton, Fla., 1993, pp. 305-363.

[0531] [49] H. G. FEICHTINGER. K. GRÖCHENIG, AND T. STROHMER, Efficient numerical methods in non-uniform sampling theory, Numer. Math., 69 (1995), pp. 423-440.

[0532] [50] H. G. FEICHTINGER AND T. STROHMER. EDS., Gabor Analysis and Algorithms, Birkhäuser, Boston, 1998.

[0533] [51] S. S. GOH AND I. G. H. ONG, Reconstruction of bandlimited signals from irregular samples, Signal Process., 46 (1995), pp. 315-329.

[0534] [52] T. N. T. GOODMAN, S. L. LEE AND W. S. TANG, Wavelet bases for a set of commuting unitary operators, Adv. Comput. Math. 1 (1993), pp. 109-126.

[0535] [53] T. N. T. GOODMAN, S. L. LEE AND W. S. TANG, Wavelet wandering subspaces, Trans. Amer. Math. Soc., 338 (1993), pp. 639-654.

[0536] [54] K. GRÖCHENIG, Describing functions: Atomic decompositions versus frames, Monatsh. Math., 112 (1991), pp. 1-42.

[0537] [55] K. GRÖCHENIG, Reconstruction algorithms in irregular sampling, Math. Comp., 59 (1992), pp. 181-194.

[0538] [56] K. GRÖCHENIG, Acceleration of the frame algorithm, IEEE Trans. Signal Process., Special Issue on Wavelets and Signal Processing, 41 (1993), pp. 3331-3340.

[0539] [57] K. GRÖCHENIG, Invertibility of the Frame Operator and Banach Frames, preprint, 2001.

[0540] [58] K. GRÖCHENIG AND H. RAZAFINJATOVO, On Landau's necessary density conditions for sampling and interpolation of band-limited functions, J. London Math. Soc., 54 (1996), pp. 557-565.

[0541] [59] K. GRÖCHENIG AND T. STROHMER, Numerical and theoretical aspects of non-uniform sampling of band-limited images, in Theory and Practice of Nonuniform Sampling, F. Marvasti, ed., Kluwer/Plenum, New York. 2001, in press.

[0542] [60] C. E. HEIL AND D. F. WALNUT, Continuous and discrete wavelet transforms, SIAM Rev., 31 (1989), pp. 628-666.

[0543] [61] J. R. HIGGINS, Five short stories about the cardinal series, Bull. Amer. Math. Soc., 121 (1985), pp. 45-89.

[0544] [62] J. R. HIGGINS, Sampling theory for Paley-Wiener spaces in the Riesz basis settings, Proc. Roy. Irish Acad. Sect., 2 (1994), pp. 219-235.

[0545] [63] H. S. Hou AND H. C. ANDREWS, Cubic splines for image interpolation and digital filtering, IEEE Trans. Acoust. Speech Signal Process. 26 (1978), pp. 508-517.

[0546] [64] A. J. E. M. JANSSEN, The Zak transform and sampling theorems for wavelet subspaces, IEEE Trans. Signal Process., 41 (1993), pp. 3360-3364.

[0547] [65] A. JERRI, The Shannon sampling theorem-its various extensions and applications: A tutorial review, Proc. IEEE, 65 (1977), pp. 1565-1596.

[0548] [66] K. JETTER AND J. STÖCKLER, Topics in scattered data interpolation and non-uniform sampling, in Surface Fitting and Multiresolution Methods, A. Le Mehautee, C. Rabut, and L. L. Schumaker, eds., Vanderbilt University Press, Nashville, Tenn., 1997, pp. 191-207.

[0549] [67] R.-Q. JIA, Shift-invariant spaces and linear operator equations, Israel Math. J., 103 (1998), pp. 259-288.

[0550] [68] R.-Q. JIA, Stability of the shifts of a finite number of functions, J. Approx. Theory, 95 (1998), pp. 194-202.

[0551] [69] R.-Q. JIA AND C. A. MICCHELLI, On linear independence of integer translates of a finite number of functions, Proc. Edinburgh Math. Soc., 36 (1992), pp. 69-85.

[0552] [70] P. JORGENSEN, A geometric approach to the cascade approximation operator for wavelets, Integral Equations Operator Theory, 35 (1999), pp. 125-171.

[0553] [71] M I. KADEG, The exact value of the Paley-Wiener constant, Soviet Math. DokI., 5 (1964), pp. 559-561.

[0554] [72] H. P. KRAMER, A generalized sampling theorem, J. Math. Phys., 38 (1959), pp. 68-72.

[0555] [73] H. LANDAU, A sparse regular sequence of exponentials closed on large sets, Bull. Amer. Math. Soc., 70 (1964), pp. 566-569.

[0556] [74] H. LANDAU, Necessary density conditions for sampling and interpolation of certain entire functions, Acta Math., 117 (1967), pp. 37-52.

[0557] [75] S. Li, Iterative Irregular Sampling and Useful Irregular Sampling Conditions, preprint.

[0558] [76] Y. Liu, Irregular sampling for spline wavelet subs paces, IEEE Trans. Inform. Theory, 42 (1996), pp. 623-627.

[0559] [77] Y. Liu AND G. G. WALTER, Irregular sampling in wavelet subs paces, J. Fourier Anal. Appl., 2 (1996), pp. 181-189.

[0560] [78] Y. LYUBARSKI{haeck over (I)} AND W. R. MADYCH, The recovery of irregularly sampled band limited junctions via tempered splines, J. Funct. Anal., 125 (1994), pp. 201-222.

[0561] [79] Y. LYUBARSKI{haeck over (I)} AND K. SEIP, Convergence and summability of Gabor expansions at the Nyqyuist density, J. Fourier Anal. AppI., 5 (1999), pp. 127-157.

[0562] [80] S. MALLAT, Multiresolution approximations and wavelet orthonormal bases of L2(R), Trans. Amer. Math. Soc., 315 (1989), pp. 69-97.

[0563] [81] S. MALLAT, A theory for multiresolutiori signal decomposition: The wavelet representation, IEEE Trans. PAMI, 11 (1989), pp. 674-693.

[0564] [82] S. MALLAT, A Wavelet Tour, Academic Press, New York, 1996.

[0565] [83] Y. MEYER, Ondelettes et Opérateurs, Hermann, Paris, France, 1990.

[0566] [84] M. Z. NASHED AND G. G. WALTER, General sampling theorems for functions in reproducing kernel Hilbert spaces, Math. Control Signals Systems, 4 (1991), pp. 363-390.

[0567] [85] J. P. OAKLEY, M. J. CUNNIGHAM AND G. LITTLE, A Fourier-domain formula for the least squares projection of a function onto a repetitive basis in N-dimensional space, IEEE Trans. Acoust. Speech Signal Process., 38 (1990), pp. 114-120.

[0568] [86] J. L. OSTUNI, A. K. S. SANTHAAND, V. S. MATTAY, D. R. WEINBERGER, R. L. LEVIN AND J. A. FRANK, Analysis of interpolation effects in the reslicing of functional MR-Images, J. Computer Assisted Tomography, 21 (1997), pp. 803-810.

[0569] [87] R. E. A. C. PALEY AND N. WIENER, Fourier transform in the complex domain, in Amer. Math. Soc. Colloq. PubI., AMS, Providence, R.I., 1934.

[0570] [88] A. PAPOULIS, Generalized sampling expansions, Circuits Systems, 24 (1977), pp. 652-654.

[0571] [89] B. S. PAVLOV, The basis property of a system of exponential.s and the condition of Muckenhoupt, Dokl. Akad. Nauk SSSR, 247 (1979), pp. 37-40.

[0572] [90] D. POTTS AND C. STEIDL, New Fourier reconstruction algorithms for computerized tomography, in Proc. Wavelet Applications in Signal and Image Processing VIII, 4199 (2000), pp. 13-23.

[0573] [91] H. REITER, Classical Harmonic Analysis and Locally Compact Groups, Oxford Univ. Press, Oxford, UK, 1968.

[0574] [92] M. RAUTH AND T. STROHMER, Smooth approximation of potential fields from noisy scattered data, Geophysics, 63 (1998), pp. 85-94.

[0575] [93] K. SKIP, An irregular sampling theorem for functions bandlemeted en a generalezed sense, SIAM J. Appl. Math., 47 (1987), pp. 1112-1116.

[0576] [94] K. SKIP, On the connection between exponential bases and certain related sequences in L2(−π, π), J. Funct. Anal., 130 (1995), pp. 131-160.

[0577] [95] I. W. SELESNICK, Interpolating multiwavelet and the sampling theorem, IEEE Trans. Signal Process., 46 (1999), pp. 2898-2908.

[0578] [96] C. E. SHANNON, Communications in the presence of noise, Proc. IRE, 37 (1949), pp. 10-21.

[0579] [97] J. STÖCKLER, Multivariate Affine Frames, Shaker Verlag, Aachen, Germany, 1998.

[0580] [98] G. STRANG, The finite element method and approximation theory, in Proc. Sympos. Numerical Solution of Partial Differential Equations, Academic Press, New York, 11(1971), pp. 547-583.

[0581] [99] G. STRANG, Wavelets and dilation equations: A brief introduction, SIAM Rev., 31 (1989), pp. 614-627.

[0582] [100] G. STRANG AND T. NGUYEN, Wavelets and Filter Banks, Wellesley-Cambridge Press, Wellesley, Mass., 1996.

[0583] [101] T. STROHMER, Numerical analysis of the non-uniform sampling problem, J. Comput. Appi. Math., 122 (2000), pp. 297-316.

[0584] [102] W. SUN AND X. ZHOU, On the stability of multivariate trigonometric systems, J. Math. Anal. Appl., 235 (1999), pp. 159-167.

[0585] [103] P. THÉVENAZ, T. BLU, AND M. UNSER, Image interpolation and resampling, in Handbook of Medical Image Processing, Processing and Analysis, I. N. Bankman, ed., Academic Press, San Diego, Calif. 2000, pp. 393-420.

[0586] [104] C. THOMAS, A comparison of motion-compensated interlace-to-progressive conversion methods, Image Communication, 12 (1998), pp. 209-229.

[0587] [105] M. UNSER AND A. ALDROUBI, Polynomial splines and wavelets-a signal processing perspective, in Wavelets: A Tutorial in Theory and Applications, C. K. Chui, ed., Academic Press, San Diego, Calif., 1992, pp. 543-601.

[0588] [106] M. UNSER AND A. ALDROUBI, Ceneralized sampling with application to the wavelet transform, in Proceedings of the Conference on Information Sciences and Systems, The Johns Hopkins University Press, Baltimore, Md., 1993.

[0589] [107] M. UNSER AND A. ALDROUBI, A general sampling theory for non-ideal acquisition devices, IEEE Trans. Signal Process., 42 (1994), pp. 2915-2925.

[0590] [108] M. UNSER, A. ALDROUBI, AND M. EDEN, A sampling theory for polynomial splines, in Proceedings of the International Conference on Information Theory and Its Applications, Waikiki, Hi., 1990, pp. 279-282.

[0591] [109] M. UNSER, A. ALDROUBI, AND M. EDEN, Polynomial spline signal approximations: Filter design and asymptotic equivalence with Shannon's sampling theorem, IEEE Trans. Image Process., 38 (1992), pp. 95-103.

[0592] [110] M. UNSER AND J. ZERUBIA, A generalized sampling theory without bandlimiting constraints, Trans. Circuits and Systems-II: Analog and Digital Signal Processing, 45 (1998), pp. 959-969.

[0593] [111] R. VIO, T. STROHMER. AND W. WAMSTEKER. On the reconstruction of irregularly sampled time series. PubI. Astronom. Soc. Pac., 112 (2000), pp. 74-90.

[0594] [112] D. WALUNT, Nonperiodic sampling of bandlimited functions on union of rectangular, J. Fourier Anal. Appl., 2 (1996), pp. 436-451.

[0595] [113] G. G. WALTER, A sampling theorem for wavelet subspaces, IEEE Trans. Inform. Theory, 38 (1992), pp. 881-884.

[0596] [114] J. M. WHITTAKER, Interpolatory Function Theory, Cambridge University Press, London, 1935.

[0597] [115] X. G. XIA, J. S. CERONIMO, D. P. HARDIN AND B. W. SUTTER, Design of prefilters for discrete multiwavelet transforms, IEEE Trans. Signal Process., 44 (1996), pp. 25-35.

[0598] [116] X.G. XIA AND Z.Z. ZHANG, On sampling theorem, wavelets, and wavelet transforms, IEEE Trans. Signal Process., 41(1993), pp. 3524-3535.

[0599] [117] K. YAO, Applications of reproducing kernel Hilbert spaces-bandlimited signal models, Inform. and Control, 11 (1967), pp. 429-444.

[0600] [118] R. M. YOUNG, An Introduction to Nonharmonic Fourier Series, Academic Press, New York, 1980.

[0601] [119] X. ZHOU AND W. SUN, On the sampling theorem for wavelet subspaces, J. Fourier Anal. Appl., 5 (1999), pp. 347-354.

[0602] [120] A. ALDROUBI AND K. GROCHENIG, Beurling-Landau-type theorems for non-uniform sampling in shift invariant spline spaces, J. Fourier Anal. Appl. 6 (2000), pp. 93-103.

[0603] [121] M. D. BUHMANN, Radial basis functions, Acta numerica, (2000) Cambridge Univ. Press, Cambridge, pp. 1-38.

[0604] [122] P. CRAVEN AND GRACE WAHBA, Smoothing noisy data with spline functions. Estimating the correct degree of smoothing by the method of generalized cross-validation, Numer. Math. 31, (1978/79), pp. 377-403

[0605] [123] S. DEMKO, W. F. Moss AND P. W. SMITH, Decay rates for inverses of band matrices, Math. Comp. 43 (1984), pp. 491-499.

[0606] [124] G. H. GOLUB AND C. F. VAN LOAN, Matrix computations, third ed., Johns Hopkins Universtiy Press, Baltimore, Md. (1996)

[0607] [125] M. F. HUTCHINSON AND F. R. DEHOOG, Smoothing noisy data with spline functions, Numer. Math., 47 (1985), pp. 99-106

[0608] [126] S. JAFFARD, Propriétés des matrices “bien localisées” près de leur diagonale et quelques applications, Ann. Inst. H. Poincaré Anal. Non Liéaire, 7, (1990), pp. 461-476

[0609] [127] L. L. SCHUMAKER, Spline functions: basic theory, John Wiley & Sons Inc., New York, 1981, Pure and Applied Mathematics, A Wiley-Interscience Publication.

[0610] [128] G. WAHBA, Spline models for observational data, Society for Industrial and Applied Mathematics SIAM, Philadelphia, Pa. (1990)

Claims

1. A method for sampling data for an object of interest, which is characterized by a functions, from a data set {f(xj)}obtained from data sample points xjε[M0, M1I], j being an indexing number, [M0, M1] representing an interval where the function f is sampled, comprising the steps of: a. constructing a first matrix =M0,M1 with matrix elements in the form of jk=φ(xj−k), wherein φ is a continuous function with compact support of size S, S being an integer, j being an integer from 1 to J, J being an integer ≧M1−M0+2S−1, and k being an integer from M0−S+1 to M1+S−1; b. obtaining a second matrix =M0,M1==, with matrix elements in the form of 97𝒥k⁢ ⁢l=∑j=1J⁢φ⁡(xj-k)_⁢φ⁡(xj-l),wherein each of k, l is an integer from M0−S+1 to M1+S−1; c. obtaining a first vector b=*y with vector elements in the form of 98bk=∑j=1J⁢φ⁡(xj-k)_⁢yjfor k=M0−S+1, . . . , M1+S−1, wherein y=(y1, . . . , yJ) is a sampling vector; d. solving a plurality of equations in the form of c=1b to obtain a second vector c=(ck), wherein k is an integer from M0−S+1 to M1+S−1; e. computing the restriction of the function f to [M0, M1] to obtain a new data set, f(x), in the form of 99f⁡(x)=∑k=M0-S+1M1+S-1⁢ck⁢φ⁡(x-k)for xε[M0,M1];and f. reconstructing the object of interest f from the new data set.

2. The method of claim 1, wherein the step of reconstructing the object of interest f from the new data set comprises the step of obtaining the object of interest f as the least square approximation of the sampling vector y from

3. The method of claim 2, wherein when y arises as the sampling vector, y={f(xj)}, an exact reconstruction of f is obtained from the new data set.

4. The method of claim 1, wherein the data set {f(xj)}comprises a plurality of digital data related to sound.

5. The method of claim 1, wherein the data set {f(xj)}comprises a plurality of digital data related to at least one image.

6. The method of claim 5, wherein the at least one image is one of a plurality of magnetic resonance image, a plurality of computerized tomography images, a plurality of optical image, a plurality of ultra sound images, a plurality of electronic images that are transmittable over a network, a plurality of satellite images, a plurality of three-dimensional images, a plurality of spectral images, a plurality of n-dimensional images with n being an integer greater than 1, or any combination of them.

7. The method of claim 1, wherein the data set {f(xj)}comprises a plurality of biological data.

8. The method of claim 7, wherein the plurality of biological data contain a plurality of data from the group of a plurality of biological data related to genes, a plurality of biological data related to proteins, a plurality of biological data related to cells, a plurality of biological data related to bacteria, or a plurality of biological data related to tissues.

9. The method of claim 1, wherein the data set {f(xj)} comprises a plurality of data from at least one measurement of the object of the interest.

10. The method of claim 9, wherein the at least one measurement is one of at least one astronomical measurement, at least one seismic measurement, at least one marine measurement, at least one geophysical measurement, at least one atmospheric measurement, at least one engineering measurement, at least one physical measurement, or at least one chemical measurement.

11. The method of claim 1, wherein the data set {f(xj)} comprises a plurality of data related to time dependent signals.

12. A system for sampling data for an object of interest, which is characterized by a function f, from a data set {f(xj)}obtained from data sample points xjε[M0, M1], j being an indexing number, [M0, M1] representing an interval where the function f is sampled, comprising a processor, wherein the processor is capable of performing the steps of: a. constructing a first matrix =M0,M1 with matrix elements in the form of jk=φ(xj−k), wherein φ is a continuous function with compact support of size S, S being an integer, j being an integer from 1 to J, J being an integer ≧M1−M0+2S−1, and k being an integer from M0−S+1 to M1+S−1; b. obtaining a second matrix =M0,M1=*, with matrix elements in the form of 101𝒥k⁢ ⁢l=∑j=1J⁢φ⁡(xj-k)_⁢φ⁡(xj-l),wherein each of k, l is an integer from M0−S+1 to M1+S−1; c. obtaining a first vector b=*y with vector elements in the form of 102bk=∑j=1J⁢φ⁡(xj-k)_⁢yjfor k=M0−S+1, . . . , M1+S−1, wherein y=(y1, . . . , yJ) is a sampling vector; d. solving a plurality of equations in the form of c=−1b to obtain a second vector c=(ck), wherein k is an integer from M0−S+1 to M1+S−1; i. computing the restriction of the function f to [M0,M1] to obtain a new data set, f(x), in the form of 103f⁡(x)=∑k=M0-S+1M1+S-1⁢ck⁢φ⁡(x-k)⁢ ⁢for⁢ ⁢x∈[M0,M1];andii. reconstructing the object of interest f from the new data set.

13. The system of claim 12, wherein the step of reconstructing the object of interest f from the new data set comprises obtaining the object of interest f as the least square approximation of the sampling vector y from

14. The system of claim 13, wherein when y arises as the sampling vector, y={f(xj)}, an exact reconstruction of f is obtained from the new data set.

15. The system of claim 12, wherein the data set {f(xj)}comprises a plurality of digital data related to sound.

16. The system of claim 12, wherein the data set {f(xj)}comprises a plurality of digital data related to at least one image.

17. The system of claim 16, wherein the at least one image is one of a plurality of magnetic resonance image, a plurality of computerized tomography images, a plurality of optical image, a plurality of ultra sound images, a plurality of electronic images that are transmittable over a network, a plurality of satellite images, a plurality of three-dimensional images, a plurality of spectral images, a plurality of n-dimensional images with n being an integer greater than 1, or any combination of them.

18. The system of claim 12, wherein the data set {f(xj)}comprises a plurality of biological data.

19. The system of claim 18, wherein the plurality of biological data contain a plurality of data from the group of a plurality of biological data related to genes, a plurality of biological data related to proteins, a plurality of biological data related to cells, a plurality of biological data related to bacteria, or a plurality of biological data related to tissues.

20. The system of claim 12, wherein the data set {f(xj)} comprises a plurality of data from at least one measurement of the object of the interest.

21. The system of claim 20, wherein the at least one measurement is one of at least one astronomical measurement, at least one seismic measurement, at least one marine measurement, at least one geophysical measurement, at least one atmospheric measurement, at least one engineering measurement, at least one physical measurement, or at least one chemical measurement.

22. The system of claim 12, wherein the data set {f(xj)} comprises a plurality of data related to time dependent signals.

23. The system of claim 12, wherein the processor comprises a microprocessor.

24. The system of claim 12, further comprising an input device coupled to the processor.

25. The system of claim 24, wherein the input device comprises at least one device selected from the group of a processor interface, a GUI, a scanner, a CD-ROM, a diskette, a computer coupled to a network, and a networking device.

26. The system of claim 12, further comprising an output device coupled to the processor.

27. The system of claim 26, wherein the output device comprises at least one device selected from the group of a GUI, a printer, a CD-ROM, a diskette, a memory device, a computer coupled to a network, and a networking device.

28. A method for converting data for an object of interest, which is characterized by a functions between a digital form and an analog form, comprising the steps of: a. locally selecting a plurality of data sample points in the order of N, N being an integer, wherein the data sample points are in a first form of data type; b. performing a transformation in a shift-invariant space to the locally selected data sample points to obtain a new data set that is in a second form of data type different from the first form of data type; and c. reconstructing the object of interest f from the new data set.

29. The method of claim 28, wherein the first form of data type is one of the digital form and the analog form.

30. The method of claim 28, wherein the second form of data type is one of the digital form and the analog form.

31. The method of claim 28, wherein the step of locally selecting a plurality of data sample points comprises the step of obtaining a data set {f(xj)} from data sample points xjε[M0, M1], j being an indexing number up to the order of N, [M0, M1] representing an interval where the function fis sampled.

32. The method of claim 31, wherein the step of performing a transformation in a shift-invariant space comprises the steps of: a. constructing a first matrix =M0,M1 with matrix elements in the form of jk=φ(xj−k), wherein φ is a continuous function with compact support of size S, S being an integer, j being an integer from 1 to J, J being an integer ≧M1−M0+2S−1, and k being an integer from M0−S+1 to M1+S−1; b. obtaining a second matrix =M0,M1=*, with matrix elements in the form of 105𝒥k⁢ ⁢l=∑j=1J⁢φ⁡(xj-k)_⁢φ⁡(xj-l),wherein each of k, l is an integer from M0−S+1 to M1+S−1; c. obtaining a first vector b=*y with vector elements in the form of 106bk=∑j=1J⁢φ⁡(xj-k)_⁢yjfor k=M0−S+1, . . . , M1+S−1, wherein y=(y1, . . . , yJ) is a sampling vector; d. solving a plurality of equations in the form of c=−1b to obtain a second vector c=(ck), wherein k is an integer from M0−S+1 to M1+S−1; and e. computing the restriction of the function f to [M0, M1] to obtain the new data set, f(x), in the form of 107f⁡(x)=∑k=M0-S+1M1+S-1⁢ck⁢φ⁡(x-k)⁢ ⁢for⁢ ⁢x∈[M0,M1].

33. The method of claim 32, wherein the step of reconstructing the object of interest f from the new data set comprises the step of obtaining the object of interest f as the least square approximation of the sampling vector y from

34. The method of claim 28, wherein the data set {f(xj)}comprises a plurality of data related to sound.

35. The method of claim 28, wherein the data set {f(xj)}comprises a plurality of data related to at least one image.

36. The method of claim 35, wherein the at least one image is one of a plurality of magnetic resonance image, a plurality of computerized tomography images, a plurality of optical image, a plurality of ultra sound images, a plurality of electronic images that are transmittable over a network, a plurality of satellite images, a plurality of three-dimensional images, a plurality of spectral images, a plurality of n-dimensional images with n being an integer greater than 1, or any combination of them.

37. The method of claim 28, wherein the data set {f(xj)} comprises a plurality of biological data.

38. The method of claim 37, wherein the plurality of biological data contain a plurality of data from the group of a plurality of biological data related to genes, a plurality of biological data related to proteins, a plurality of biological data related to cells, a plurality of biological data related to bacteria, or a plurality of biological data related to tissues.

39. The method of claim 28, wherein the data set {f(xj)} comprises a plurality of data from at least one measurement of the object of the interest.

40. The method of claim 39, wherein the at least one measurement is one of at least one astronomical measurement, at least one seismic measurement, at least one marine measurement, at least one geophysical measurement, at least one atmospheric measurement, at least one engineering measurement, at least one physical measurement, or at least one chemical measurement.

41. The method of claim 28, wherein the data set {f(xj)} comprises a plurality of data related to time dependent signals.

42. A system for converting data for an object of interest, which is characterized by a functions, between a digital form and an analog form, comprising a processor, wherein the processor is capable of performing the steps of: a. locally selecting a plurality of data sample points in the order of N, N being an integer, wherein the data sample points are in a first form of data type; b. performing a transformation in a shift-invariant space to the locally selected data sample points to obtain a new data set that is in a second form of data type different from the first form of data type; and C. reconstructing the object of interest f from the new data set.

43. The system of claim 42, wherein the first form of data type is one of the digital form and the analog form.

44. The system of claim 42, wherein the second form of data type is one of the digital form and the analog form.

45. The system of claim 42, wherein the step of locally selecting a plurality of data sample points comprises the step of obtaining a data set {f(xj)} from data sample points xjε[M0, M1], j being an indexing number up to the order of N, [M0, M1] representing an interval where the function f is sampled.

46. The system of claim 45, wherein the step of performing a transformation in a shift-invariant space comprises the steps of: a. constructing a first matrix =M0,M1 with matrix elements in the form of jk=φ(xj−k), wherein φ is a continuous function with compact support of size S, S being an integer, j being an integer from 1 to J, J being an integer ≧M1−M0+2S−1, and k being an integer from M0−S+1 to M1+S−1; b. obtaining a second matrix =M0,M1=*, with matrix elements in the form of 109𝒥k⁢ ⁢l=∑j=1J⁢ ⁢φ⁡(xj-k)_⁢φ⁡(xj-l),wherein each of k, l is an integer from M0−S+1 to M1+S−1; c. obtaining a first vector b=*y with vector elements in the form of 110bk=∑j=1J⁢ ⁢φ⁡(xj-k)_⁢yjfor k=M0−S+1, . . . , M1+S−1, wherein y=(y1, . . . , yj) is a sampling vector; d. solving a plurality of equations in the form of c=−1b to obtain a second vector c=(ck), wherein k is an integer from M0−S+1 to M1+S−1; and e. computing the restriction of the function f to [M0, M1] to obtain the new data set, f(x), in the form of 111f⁡(x)=∑k=M0-S+1M1+S-1⁢ck⁢ ⁢φ⁡(x-k)⁢ ⁢for⁢ ⁢x∈[M0,M1].

47. The system of claim 46, wherein the step of reconstructing the object of interest f from the new data set comprises the step of obtaining the object of interest f as the least square approximation of the sampling vector y from

48. The system of claim 45, wherein the data set {f(xj)}comprises a plurality of data related to sound.

49. The system of claim 45, wherein the data set {f(xj)} comprises a plurality of data related to at least one image.

50. The system of claim 49, wherein the at least one image is one of a plurality of magnetic resonance image, a plurality of computerized tomography images, a plurality of optical image, a plurality of ultra sound images, a plurality of electronic images that are transmittable over a network, a plurality of satellite images, a plurality of three-dimensional images, a plurality of spectral images, a plurality of n-dimensional images with n being an integer greater than 1, or any combination of them.

51. The system of claim 45, wherein the data set {f(xj)}comprises a plurality of biological data.

52. The system of claim 51, wherein the plurality of biological data contain a plurality of data from the group of a plurality of biological data related to genes, a plurality of biological data related to proteins, a plurality of biological data related to cells, a plurality of biological data related to bacteria, or a plurality of biological data related to tissues.

53. The system of claim 45, wherein the data set {f(xj)} comprises a plurality of data from at least one measurement of the object of the interest.

54. The system of claim 53, wherein the at least one measurement is one of at least one astronomical measurement, at least one seismic measurement, at least one marine measurement, at least one geophysical measurement, at least one atmospheric measurement, at least one engineering measurement, at least one physical measurement, or at least one chemical measurement.

55. The system of claim 45, wherein the data set {f(xj)} comprises a plurality of data related to time dependent signals.

56. The system of claim 42, wherein the processor comprises a microprocessor.

57. The system of claim 42, further comprising an input device coupled to the processor.

58. The system of claim 57, wherein the input device comprises at least one device selected from the group of a processor interface, a GUI, a scanner, a CD-ROM, a diskette, a computer coupled to a network, and a networking device.

59. The system of claim 42, further comprising an output device coupled to the processor.

60. The system of claim 59, wherein the output device comprises at least one device selected from the group of a GUI, a printer, a CD-ROM, a diskette, a memory device, a computer coupled to a network, and a networking device.

61. A method for recovering information about an object of interest f from a data set {f(xj)} obtained from nonuniform data sample points xjεX, j being an indexing number, X representing the set of points where the object of interests is sampled, comprising the steps of: a. selecting data f(xj) from the data set {f(xj)}; b. constructing an approximation operator QX; c. applying the approximation operator QX to the data f(xj) to obtain an approximation QXf(xj); d. constructing a projection operator P; e. applying the projection operator P to the approximation QXf(xj) to obtain a first approximations f1=P QXf(xj); f. obtaining an error e=f(xj)−f1; g. applying the projection operator P to the error e to obtain a first approximation of error e1=P e; h. obtaining a second approximations f2=f1+e1; and i. returning to step (e) until a sequence fn=f1+e1+e2+e3+ . . . +en−1 is obtained, wherein function fn converges to the object of interest f.

62. The method of claim 61, wherein the approximation operator QX comprises an interpolation operator.

63. The method of claim 62, wherein the interpolation operator provides piecewise linear interpolation.

64. The method of claim 62, wherein the approximation operator QX comprises a quasi-interpolation operator.

65. The method of claim 64, wherein the quasi-interpolation operator provides step-wise approximation.

66. The method of claim 61, wherein the projection operator P comprises a bounded projection.

67. The method of claim 61, wherein the projection operator P comprises a universal projection.

68. The method of claim 61, wherein the data set {f(xj)}comprises a plurality of digital data related to sound.

69. The method of claim 61, wherein the data set {f(xj)}comprises a plurality of digital data related to at least one image.

70. The method of claim 69, wherein the at least one image is one of a plurality of magnetic resonance image, a plurality of computerized tomography images, a plurality of optical image, a plurality of ultra sound images, a plurality of electronic images that are transmittable over a network, a plurality of satellite images, a plurality of three-dimensional images, a plurality of spectral images, a plurality of n-dimensional images with n being an integer greater than 1, or any combination of them.

71. The method of claim 61, wherein the data set {f(xj)}comprises a plurality of biological data.

72. The method of claim 71, wherein the plurality of biological data contain a plurality of data from the group of a plurality of biological data related to genes, a plurality of biological data related to proteins, a plurality of biological data related to cells, a plurality of biological data related to bacteria, or a plurality of biological data related to tissues.

73. The method of claim 61, wherein the data set {f(xj)} comprises a plurality of data from at least one measurement of the object of the interest.

74. The method of claim 73, wherein the at least one measurement is one of at least one astronomical measurement, at least one seismic measurement, at least one marine measurement, at least one geophysical measurement, at least one atmospheric measurement, at least one engineering measurement, at least one physical measurement, or at least one chemical measurement.

75. The method of claim 61, wherein the data set {f(xj)} comprises a plurality of data related to time dependent signals.

76. A system for recovering information about an object of interest f from a data set {f(xj)} obtained from nonuniform data sample points xjεX, j being an indexing number, X representing the set of points where the object of interest f is sampled, comprising: (i). an input device for receiving the information in the data set {f(xj)}; and (ii). a processor communicating to the input device and performing the steps of: a. selecting data f(xj) from the data set {f(xj)}; b. constructing an approximation operator QX; c. applying the approximation operator QX to the data f(xj) to obtain an approximation QX f(xj); d. constructing a projection operator P; e. applying the projection operator P to the approximation QXf(xj) to obtain a first approximation f1=PQXf(xj); f. obtaining an error e=f(xj)−f1; g. applying the projection operator P to the error e to obtain a first approximation of error e1=P e; h. obtaining a second approximation f2=f1+e1; and i. returning to step (e) until a sequence fn=f1+e1+e2+e3+ . . . +en−1 is obtained, wherein function fn converges to the object of interest f.

77. The system of claim 76, wherein the approximation operator QX comprises an interpolation operator.

78. The system of claim 77, wherein the interpolation operator provides piecewise linear interpolation.

79. The system of claim 77, wherein the approximation operator QX comprises a quasi-interpolation operator.

80. The system of claim 79, wherein the quasi-interpolation operator provides step-wise approximation.

81. The system of claim 76, wherein the projection operator P comprises a bounded projection.

82. The system of claim 76, wherein the projection operator P comprises a universal projection.

83. The system of claim 76, wherein the processor comprises a microprocessor.

84. The system of claim 76, wherein the input device comprises at least one device selected from the group of a processor interface, a GUI, a scanner, a CD-ROM, a diskette, a computer coupled to a network, and a networking device.

85. The system of claim 76, further comprising an output device coupled to the processor.

86. The system of claim 85, wherein the output device comprises at least one device selected from the group of a GUI, a printer, a CD-ROM, a diskette, a memory device, a computer coupled to a network, and a networking device.

87. The system of claim 76, wherein the data set {f(xj)}comprises a plurality of digital data related to sound.

88. The system of claim 76, wherein the data set {f(xj)}comprises a plurality of digital data related to at least one image.

89. The system of claim 88, wherein the at least one image is one of a plurality of magnetic resonance image, a plurality of computerized tomography images, a plurality of optical image, a plurality of ultra sound images, a plurality of electronic images that are transmittable over a network, a plurality of satellite images, a plurality of three-dimensional images, a plurality of spectral images, a plurality of n-dimensional images with n being an integer greater than 1, or any combination of them.

90. The system of claim 76, wherein the data set {f(xj)}comprises a plurality of biological data.

91. The system of claim 90, wherein the plurality of biological data contain a plurality of data from the group of a plurality of biological data related to genes, a plurality of biological data related to proteins, a plurality of biological data related to cells, a plurality of biological data related to bacteria, or a plurality of biological data related to tissues.

92. The system of claim 76, wherein the data set {f(xj)} comprises a plurality of data from at least one measurement of the object of the interest.

93. The system of claim 92, wherein the at least one measurement is one of at least one astronomical measurement, at least one seismic measurement, at least one marine measurement, at least one geophysical measurement, at least one atmospheric measurement, at least one engineering measurement, at least one physical measurement, or at least one chemical measurement.

94. The system of claim 76, wherein the data set {f(xj)} comprises a plurality of data related to time dependent signals.

95. A method for recovering information about an object of interests from a data set {f(xj)} obtained from nonuniform data sample points xjεX, j being an indexing number, X representing the set of points where the object of interest f is sampled in presence of noise, comprising: a. selecting data fj′ from the data set {fj′=f(xj)+ηj}, wherein ηj represents a corresponding noise component; b. constructing an initialization function f as the summation of function fj′: f=ΣjεJfj′βj, wherein βj represents the jth component of a partition of unity; c. constructing an approximation operator QX; d. applying the approximation operator QX to the function f′ to obtain an approximation QXf′; e. constructing a projection operator P; f. applying the projection operator P to the approximation QXf′ to obtain a first approximation f1=P QXf′; g. obtaining an error e=f−f1; h. applying the projection operator P to the error e to obtain a first approximation of error e1=P e; i. obtaining a second approximation f2=f1+e1; and j. returning to step (f) until a sequence fn=f1+e1+e2+e3+ . . . +en−1 is obtained, wherein function fn converges to a function f that adequately describes the object of interest f.

96. The method of claim 95, wherein the approximation operator QX comprises an interpolation operator.

97. The method of claim 96, wherein the interpolation operator provides piecewise linear interpolation.

98. The method of claim 95, wherein the approximation operator QX comprises a quasi-interpolation operator.

99. The method of claim 98, wherein the quasi-interpolation operator provides step-wise approximation.

100. The method of claim 95, wherein the projection operator P comprises a bounded projection.

101. The method of claim 95, wherein the projection operator P comprises a universal projection.

102. The method of claim 95, wherein the function f satisfies the relationship of P QXf∞=P QX{fj′}.

103. The method of claim 95, wherein the data set {f(xj)}comprises a plurality of digital data related to sound.

104. The method of claim 95, wherein the data set {f(xj)}comprises a plurality of digital data related to at least one image.

105. The method of claim 104, wherein the at least one image is one of a plurality of magnetic resonance image, a plurality of computerized tomography images, a plurality of optical image, a plurality of ultra sound images, a plurality of electronic images that are transmittable over a network, a plurality of satellite images, a plurality of three-dimensional images, a plurality of spectral images, a plurality of n-dimensional images with n being an integer greater than 1, or any combination of them.

106. The method of claim 95, wherein the data set {f(xj)} comprises a plurality of biological data.

107. The method of claim 106, wherein the plurality of biological data contain a plurality of data from the group of a plurality of biological data related to genes, a plurality of biological data related to proteins, a plurality of biological data related to cells, a plurality of biological data related to bacteria, or a plurality of biological data related to tissues.

108. The method of claim 95, wherein the data set {f(xj)} comprises a plurality of data from at least one measurement of the object of the interest.

109. The method of claim 108, wherein the at least one measurement is one of at least one astronomical measurement, at least one seismic measurement, at least one marine measurement, at least one geophysical measurement, at least one atmospheric measurement, at least one engineering measurement, at least one physical measurement, or at least one chemical measurement.

110. The method of claim 95, wherein the data set {f(xj)} comprises a plurality of data related to time dependent signals.

111. A system for recovering information about an object of interest f from a data set {f(xj)} obtained from nonuniform data sample points xjεX, j being an indexing number, X representing the set of points where the object of interest f is sampled in presence of noise, comprising: (i). an input device for receiving the information in the data set {f(xj)}; and (ii). a processor communicating to the input device and performing the steps of: a. selecting data fj′ from the data set {fj′=f(xj)+ηj}, wherein ηj represents a corresponding noise component; b. constructing an initialization function f′ as the summation of function fj′:f′=93jεJfj′βj, wherein βj represents the jth component of a partition of unity; c. constructing an approximation operator QX; d. applying the approximation operator QX to the function f′ to obtain an approximation QXf; e. constructing a projection operator P; f. applying the projection operator P to the approximation QXf′ to obtain a first approximation f1=P QXf′; g. obtaining an error e=f−f1; h. applying the projection operator P to the error e to obtain a first approximation of error e1=P e; i. obtaining a second approximation f2=f1+e1; and j. returning to step (f) until a sequence fn=f1+e1+e2+e3+ . . . +en−1 is obtained, wherein function fn converges to a function f∞ that adequately describes the object of interest f.

112. The system of claim 111, wherein the approximation operator QX comprises an interpolation operator.

113. The system of claim 112, wherein the interpolation operator provides piecewise linear interpolation.

114. The system of claim 111, wherein the approximation operator QX comprises a quasi-interpolation operator.

115. The system of claim 114, wherein the quasi-interpolation operator provides step-wise approximation.

116. The system of claim 111, wherein the projection operator P comprises a bounded projection.

117. The system of claim 111, wherein the projection operator P comprises a universal projection.

118. The system of claim 111, wherein the function f∞ satisfies the relationship of PQXf∞=PQX{fj′}.

119. The system of claim 111, wherein the processor comprises a microprocessor.

120. The system of claim 111, wherein the input device comprises at least one device selected from the group of a processor interface, a GUI, a scanner, a CD-ROM, a diskette, a computer coupled to a network, and a networking device.

121. The system of claim 111, further comprising an output device coupled to the processor.

122. The system of claim 121, wherein the output device comprises at least one device selected from the group of a GUI, a printer, a CD-ROM, a diskette, a memory device, a computer coupled to a network, and a networking device.

123. The system of claim 111, wherein the data set {fj′=f(xj)+ηj} comprises a plurality of digital data related to sound.

124. The system of claim 111, wherein the data set {fj′=f(xj)+ηj} comprises a plurality of digital data related to at least one image.

125. The system of claim 124, wherein the at least one image is one of a plurality of magnetic resonance image, a plurality of computerized tomography images, a plurality of optical image, a plurality of ultra sound images, a plurality of electronic images that are transmittable over a network, a plurality of satellite images, a plurality of three-dimensional images, a plurality of spectral images, a plurality of n-dimensional images with n being an integer greater than 1, or any combination of them.

126. The system of claim 111, wherein the data set {fj′=f(xj)+ηj} comprises a plurality of biological data.

127. The system of claim 126, wherein the plurality of biological data contain a plurality of data from the group of a plurality of biological data related to genes, a plurality of biological data related to proteins, a plurality of biological data related to cells, a plurality of biological data related to bacteria, or a plurality of biological data related to tissues.

128. The system of claim 111, wherein the data set {fj′=f(xj)+ηj} comprises a plurality of data from at least one measurement of the object of the interest.

129. The system of claim 128, wherein the at least one measurement is one of at least one astronomical measurement, at least one seismic measurement, at least one marine measurement, at least one geophysical measurement, at least one atmospheric measurement, at least one engineering measurement, at least one physical measurement, or at least one chemical measurement.

130. The system of claim 111, wherein the data set {fj′=f(xj)+ηj} comprises a plurality of data related to time dependent signals.

131. A method for recovering information about an object of interest f from a data set {f(xj)}obtained from nonuniform data sample points xjεX, j being an indexing number, X representing the set of points where the object of interest f is sampled in presence of noise, comprising the steps of: a. sampling a set of data samples from the data set {f(xj)}; b. iteratively framing the sampled set of data samples until the sampled set of data samples converges to a new set of data samples that has a size significantly smaller than the data set {f(xj)}; and c. reconstructing the object of interest f from the new set of data samples.

132. The method of claim 131, wherein the data set {f(xj)}comprises a plurality of digital data related to sound.

133. The method of claim 131, wherein the data set {f(xj)}comprises a plurality of digital data related to at least one image.

134. The method of claim 133, wherein the at least one image is one of a plurality of magnetic resonance image, a plurality of computerized tomography images, a plurality of optical image, a plurality of ultra sound images, a plurality of electronic images that are transmittable over a network, a plurality of satellite images, a plurality of three-dimensional images, a plurality of spectral images, a plurality of n-dimensional images with n being an integer greater than 1, or any combination of them.

135. The method of claim 131, wherein the data set {f(xj)}comprises a plurality of biological data.

136. The method of claim 135, wherein the plurality of biological data contain a plurality of data from the group of a plurality of biological data related to genes, a plurality of biological data related to proteins, a plurality of biological data related to cells, a plurality of biological data related to bacteria, or a plurality of biological data related to tissues.

137. The method of claim 131, wherein the data set {f(xj)} comprises a plurality of data from at least one measurement of the object of the interest.

138. The method of claim 137, wherein the at least one measurement is one of at least one astronomical measurement, at least one seismic measurement, at least one marine measurement, at least one geophysical measurement, at least one atmospheric measurement, at least one engineering measurement, at least one physical measurement, or at least one chemical measurement.

139. The method of claim 131, wherein the data set {f(xj)} comprises a plurality of data related to time dependent signals.

140. A system for recovering information about an object of interests from a data set {f(xj)}obtained from nonuniform data sample points xjεX, j being an indexing number, X representing the set of points where the object of interest f is sampled in presence of noise, comprising: a. means for sampling a set of data samples from the data set {f(xj)}; b. means for iteratively framing the sampled set of data samples until the sampled set of data samples converges to a new set of data samples that has a size significantly smaller than the data set {f(xj)}; c. means for reconstructing the object of interest f from the new set of data samples.

141. The system of claim 140, wherein the data set {f(xj)}comprises a plurality of digital data related to sound.

142. The system of claim 140, wherein the data set {f(xj)}comprises a plurality of digital data related to at least one image.

143. The system of claim 142, wherein the at least one image is one of a plurality of magnetic resonance image, a plurality of computerized tomography images, a plurality of optical image, a plurality of ultra sound images, a plurality of electronic images that are transmittable over a network, a plurality of satellite images, a plurality of three-dimensional images, a plurality of spectral images, a plurality of n-dimensional images with n being an integer greater than 1, or any combination of them.

144. The system of claim 140, wherein the data set {f(xj)}comprises a plurality of biological data.

145. The system of claim 144, wherein the plurality of biological data contain a plurality of data from the group of a plurality of biological data related to genes, a plurality of biological data related to proteins, a plurality of biological data related to cells, a plurality of biological data related to bacteria, or a plurality of biological data related to tissues.

146. The system of claim 144, wherein the data set {f(xj)} comprises a plurality of data from at least one measurement of the object of the interest.

147. The system of claim 146, wherein the at least one measurement is one of at least one astronomical measurement, at least one seismic measurement, at least one marine measurement, at least one geophysical measurement, at least one atmospheric measurement, at least one engineering measurement, at least one physical measurement, or at least one chemical measurement.

148. The system of claim 140, wherein the data set {f(xj)} comprises a plurality of data related to time dependent signals.

149. The system of claim 140, wherein the sampling means comprises an input device.

150. The system of claim 140, wherein the framing means comprises a processor.

151. The system of claim 140, wherein the reconstructing means comprises an output device.

152. A method for recovering information about an object of interest f from a data set of averages {f, ψxj} obtained from nonuniform data sample points xjεX, j being an indexing number, X representing the set of points where the object of interest f is sampled, comprising the steps of: a. selecting data f, ψxj from the data set {f, ψxj}; b. constructing a quasi-reconstruction operator AX for obtaining a quasi-reconstruction AXf wherein AXf=ΣjεJf, ψxjβj, wherein βj represents the jth component of a partition of unity, and {ψxk:xjεX} is a set of functionals that act on f; c. constructing a projection operator P; d. applying the projection operator P to the quasi-reconstruction AXf to obtain a first approximation f1=PAXf; e. obtaining an error e=f−f1; f. applying the projection operator P and the quasi-reconstruction operator AX to the error e to obtain a first approximation of error e1=P AX e; g. obtaining a second approximation f2=f1+e1; and h. returning to step (e) until a sequence fn=f1+e1+e2+e3+ . . . +en'1=P AX(f−fn−1)+fn−1 is obtained, wherein function fn converges to the object of interest f.

153. The method of claim 152, wherein the data set of averages {f, ψxj} comprises a plurality of digital data related to sound.

154. The method of claim 152, wherein the data set of averages {f, ψxj} comprises a plurality of digital data related to at least one image.

155. The method of claim 154, wherein the at least one image is one of a plurality of magnetic resonance image, a plurality of computerized tomography images, a plurality of optical image, a plurality of ultra sound images, a plurality of electronic images that are transmittable over a network, a plurality of satellite images, a plurality of three-dimensional images, a plurality of spectral images, a plurality of n-dimensional images with n being an integer greater than 1, or any combination of them.

156. The method of claim 152, wherein the data set of averages {f, ψxj} comprises a plurality of biological data.

157. The method of claim 156, wherein the plurality of biological data contain a plurality of data from the group of a plurality of biological data related to genes, a plurality of biological data related to proteins, a plurality of biological data related to cells, a plurality of biological data related to bacteria, or a plurality of biological data related to tissues.

158. The method of claim 152, wherein the data set of averages {f, ψxj} comprises a plurality of data from at least one measurement of the object of the interest.

159. The method of claim 158, wherein the at least one measurement is one of at least one astronomical measurement, at least one seismic measurement, at least one marine measurement, at least one geophysical measurement, at least one atmospheric measurement, at least one engineering measurement, at least one physical measurement, or at least one chemical measurement.

160. The method of claim 152, wherein the data set of averages {f, ψxj} comprises a plurality of data related to time dependent signals.

161. A system for recovering information about an object of interest f from a data set of averages {f, ψxj} obtained from nonuniform data sample points xjεX, j being an indexing number, X representing the set of points where the object of interest f is sampled, comprising: (i). an input device for receiving the information in the data set of averages {f, ψxj}; and (ii). a processor communicating to the input device and performing the steps of: a. selecting data f, ψxj from the data set {f, ψxj}; b. constructing a quasi-reconstruction operator AX for obtaining a quasi-reconstruction AXf wherein AXf=ΣjεJf, ψxjβj, wherein βj; represents the jth component of a partition of unity, and {ψxj:xjεX}is a set of functionals that act on f; c. constructing a projection operator P; d. applying the projection operator P to the quasi-reconstruction AXf to obtain a first approximation f1=PAXf; e. obtaining an error e=f−f1; f. applying the projection operator P and the quasi-reconstruction operator Ax to the error e to obtain a first approximation of error e1=P AX e; g. obtaining a second approximation f2=f1+e1; and h. returning to step (e) until a sequence fn=f1+e1+e2+e3+ . . . +en−1=P AX(f−fn−1)+fn−1 is obtained, wherein function fn converges to the object of interest f.

162. The system of claim 161, wherein the processor comprises a microprocessor.

163. The system of claim 161, wherein the input device comprises at least one device selected from the group of a processor interface, a GUI, a scanner, a CD-ROM, a diskette, a computer coupled to a network, and a networking device.

164. The system of claim 161, further comprising an output device coupled to the processor.

165. The system of claim 164, wherein the output device comprises at least one device selected from the group of a GUI, a printer, a CD-ROM, a diskette, a memory device, a computer coupled to a network, and a networking device.

166. The system of claim 161, wherein the data set of averages {f, ψxj} comprises a plurality of digital data related to sound.

167. The system of claim 161, wherein the data set of averages {f, ψxj} comprises a plurality of digital data related to at least one image.

168. The system of claim 167, wherein the at least one image is one of a plurality of magnetic resonance image, a plurality of computerized tomography images, a plurality of optical image, a plurality of ultra sound images, a plurality of electronic images that are transmittable over a network, a plurality of satellite images, a plurality of three-dimensional images, a plurality of spectral images, a plurality of n-dimensional images with n being an integer greater than 1, or any combination of them.

169. The system of claim 161, wherein the data set of averages {f, ψxj} comprises a plurality of biological data.

170. The system of claim 169, wherein the plurality of biological data contain a plurality of data from the group of a plurality of biological data related to genes, a plurality of biological data related to proteins, a plurality of biological data related to cells, a plurality of biological data related to bacteria, or a plurality of biological data related to tissues.

171. The system of claim 161, wherein the data set of averages {f, ψxj} comprises a plurality of data from at least one measurement of the object of the interest.

172. The system of claim 171, wherein the at least one measurement is one of at least one astronomical measurement, at least one seismic measurement, at least one marine measurement, at least one geophysical measurement, at least one atmospheric measurement, at least one engineering measurement, at least one physical measurement, or at least one chemical measurement.

173. The system of claim 161, wherein the data set of averages {f, ψxj} comprises a plurality of data related to time dependent signals.

174. A method for recovering information about an object of interest f from a data set of averages {f, ψxj} obtained from nonuniform data sample points xjεX, j being an indexing number, X representing the set of points where the object of interest f is sampled in presence of noise, comprising: a. selecting data fj′ from the data set {fj′=z,903 f,ψxj+ηj}, wherein ηj represents a corresponding noise component; b. constructing an initialization function f as the summation of function fj′: f′=ΣjεJfj′βj, wherein ok represents the jth component of a partition of unity; c. constructing a quasi-interpolant operator QX; d. applying the quasi-interpolant operator QX to the function f to obtain an approximation QXf′; e. constructing a projection operator P; f. applying the projection operator P to the approximation QXf′ to obtain a first approximation f1=P QXf; g. obtaining an error e=f−f1; h. constructing a quasi-reconstruction operator AX for obtaining a quasi-reconstruction AXf wherein AXf=ΣjεJf,ψxjβj, and {ψxj:xjεX} is a set of functionals that act on f; i. applying the projection operator P and the quasi-reconstruction operator AX to the error e to obtain a first approximation of error e1=PAXe; j. obtaining a second approximation f2=f1+e1; and k. returning to step (g) until a sequence fn=f1+e1+e2+e3+ . . . +en−1=f1+(I−PAX)fn−1 is obtained, wherein I is an unit operator and function fn converges to a function f that adequately describes the object of interest f.

175. The method of claim 174, wherein the data set {fj′=f,ψxj+ηj} comprises a plurality of digital data related to sound.

176. The method of claim 174, wherein the data set {fj′=f,ψxj+ηj} comprises a plurality of digital data related to at least one image.

177. The method of claim 176, wherein the at least one image is one of a plurality of magnetic resonance image, a plurality of computerized tomography images, a plurality of optical image, a plurality of ultra sound images, a plurality of electronic images that are transmittable over a network, a plurality of satellite images, a plurality of three-dimensional images, a plurality of spectral images, a plurality of n-dimensional images with n being an integer greater than 1, or any combination of them.

178. The method of claim 174, wherein the data set {fj′=f,ψxj+ηj} comprises a plurality of biological data.

179. The method of claim 178, wherein the plurality of biological data contain a plurality of data from the group of a plurality of biological data related to genes, a plurality of biological data related to proteins, a plurality of biological data related to cells, a plurality of biological data related to bacteria, or a plurality of biological data related to tissues.

180. The method of claim 174, wherein the data set {fj′=f,ψxj+ηj} comprises a plurality of data from at least one measurement of the object of the interest.

181. The method of claim 180, wherein the at least one measurement is one of at least one astronomical measurement, at least one seismic measurement, at least one marine measurement, at least one geophysical measurement, at least one atmospheric measurement, at least one engineering measurement, at least one physical measurement, or at least one chemical measurement.

182. The method of claim 174, wherein the data set {fj′=f,ψxj+ηj} comprises a plurality of data related to time dependent signals.

183. A system for recovering information about an object of interest f from a data set of averages {f,ψxj} obtained from nonuniform data sample points xjεX, j being an indexing number, X representing the set of points where the object of interest f is sampled in presence of noise, comprising: (i). an input device for receiving the information in the data set {f(xj)}; and (ii). a processor communicating to the input device and performing the steps of: a. selecting data fj′ from the data set {fj′=f,ψxj+ηj}, wherein Tl represents a corresponding noise component; b. constructing an initialization function f′ as the summation of function fj′:f′=ΣjεJfj′βj, wherein βj represents the jth component of a partition of unity; c. constructing a quasi-interpolant operator QX; d. applying the quasi-interpolant operator QX to the functions to obtain an approximation QXf′; e. constructing a projection operator P; f. applying the projection operator P to the approximation QXf′ to obtain a first approximation f1=P QXf′; g. obtaining an error e=f−f1; h. constructing a quasi-reconstruction operator AX for obtaining a quasi-reconstruction AXf wherein AXf=τjεJf,ψxjβj, and {ψxj:xjεX}is a set of functionals that act on f(xj); i. applying the projection operator P and the quasi-reconstruction operator AX to the error e to obtain a first approximation of error e1=PAX e; j. obtaining a second approximation f2=f1+e1; and k. returning to step (h) until a sequence fn=f1+e1+e2+e3+ . . . +en−1=f1+(I−P AX)fn−1 is obtained, wherein I is an unit operator and function fn converges to a function f∞ that adequately describes the object of interest f.

184. The system of claim 183, wherein the processor comprises a microprocessor.

185. The system of claim 183, wherein the input device comprises at least one device selected from the group of a processor interface, a GUI, a scanner, a CD-ROM, a diskette, a computer coupled to a network, and a networking device.

186. The system of claim 183, further comprising an output device coupled to the processor.

187. The system of claim 186, wherein the output device comprises at least one device selected from the group of a GUI, a printer, a CD-ROM, a diskette, a memory device, a computer coupled to a network, and a networking device.

188. The system of claim 183, wherein the data set {fj′=f,ψxj+ηj} comprises a plurality of digital data related to sound.

189. The system of claim 183, wherein the data set {fj′=f,ψxj+ηj} comprises a plurality of digital data related to at least one image.

190. The system of claim 189, wherein the at least one image is one of a plurality of magnetic resonance image, a plurality of computerized tomography images, a plurality of optical image, a plurality of ultra sound images, a plurality of electronic images that are transmittable over a network, a plurality of satellite images, a plurality of three-dimensional images, a plurality of spectral images, a plurality of n-dimensional images with n being an integer greater than 1, or any combination of them.

191. The system of claim 183, wherein the data set {fj′=f,ψxj+ηj} comprises a plurality of biological data.

192. The system of claim 191, wherein the plurality of biological data contain a plurality of data from the group of a plurality of biological data related to genes, a plurality of biological data related to proteins, a plurality of biological data related to cells, a plurality of biological data related to bacteria, or a plurality of biological data related to tissues.

193. The system of claim 183, wherein the data set {fj′=f,ψxj+ηj} comprises a plurality of data from at least one measurement of the object of the interest.

194. The system of claim 193, wherein the at least one measurement is one of at least one astronomical measurement, at least one seismic measurement, at least one marine measurement, at least one geophysical measurement, at least one atmospheric measurement, at least one engineering measurement, at least one physical measurement, or at least one chemical measurement.

195. The system of claim 183, wherein the data set {fj′=f,ψxj+ηj} comprises a plurality of data related to time dependent signals.