From periodic sampling to irregular sampling through PNS (Periodic Nonuniform Sampling)
aa r X i v : . [ phy s i c s . d a t a - a n ] M a y From periodic sampling to irregular samplingthrough PNS (Periodic Nonuniform Sampling)
Bernard Lacaze T´eSA, 7 bd de la Gare, 31500 [email protected] 28, 2019
Abstract
Resampling is an operation costly in calculation time and accuracy.It regularizes irregular sampling, replacing N data by N periodic estima-tions. This stage can be suppressed, using formulas built with incomingdata and completed by sequences of elements which influence decreaseswhen the number of data increases. Obviously, some spectral proper-ties (for processes) and some asymptotic properties (for the sampling se-quence) have to be fulfilled. In this paper, we explain that it is possible toelaborate a logical theory, starting from the ordinary periodic sampling,and generalized by the PNS (Periodic Nonuniform Sampling), to treatmore general irregular samplings. The ”baseband spectrum” hypothe-sis linked to the ”Nyquist bound” (or Shannon bound) is generalized tospectra in a finite number of intervals, suited to the ”Landau condition”. keywords : periodic sampling, periodic nonuniform sampling, irregularsampling, Vostok ice core. In signal theory, sampling addresses a time sequence t = { t n , n ∈ A} and values g ( t n ) of some function g ( t ) , belonging to a given class B . We look for approx-imations e g ( t ) close to g ( t ) from the set of g ( t n ) , n ∈ A . Among possibilities,we have linear combinations such as e g ( t ) = X k ∈ Z a k ( t ) g ( t k ) . (1)Coefficients a k ( t ) characterize the class B . e g ( t ) may be equal to g ( t ) or not, andmay minimize some given distance. Practically, formula (1) will be truncated(apodized) when A is not finite, to provide usable results.The sampling sequence t= { nT, n ∈ Z } defines the simplest ”periodic sam-pling”. Bibliography on this basic situation goes back to A. L. Cauchy andE. Borel [1], [2]. The ”sampling formula” of Shannon or Nyquist (or oth-ers) refers to functions with Fourier transform cancelling outside the interval( − / T, / T ). As explained below, the problem is closely linked to usual FSD(Fourier Series Developments). PNS N (Periodic Nonuniform Sampling of or-der N ) are concatenations of N periodic samplings of same period. They were1ntroduced in the 1950s, firstly as variations on the purely periodic case [3],[4]. They allow to bypass aliasing problems, to simplify and to produce cheaperelectronical devices. The particular case of Coset (or Multi-Coset) assumes theexistence of a primary clock. Cosets are produced from subsequences of sameperiod [3], [5], [6]. They used in the domain of communications, but algorithmscan be disputed.Given some function g ( t ) , we define the set g n ( t ) , n = 1 , , .., N, as outputsof bandpass filters on disjoined sets ∆ n , where ∆ = ∪ Nn =1 ∆ n is the ”spectralsupport” of g ( t ) . For a large collection of ∆ n , we explain that a PNS N isable to recover the g n ( t ) . Extending sampling properties from ordinary func-tions to random processes is very easy. In both cases, we have to develop thesame exponential functions, but with convergence criteria which can be dif-ferent. Therefore, we prove that PNS N, when N is large enough, is able todetermine bandpass filterings and power spectra of stationary processes. Par-ticularly, ice cores provide a wide scope of data at irregular times and allow toverify improvements thanks to methods displayed in this paper.This paper is organized as follows. Section 2 explains the periodic samplingas dependent of spectral supports. Section 3 develops the PNS N for particularpartitions of spectra. Subsequent sections extend these results for more generalsamplings, and section 8 specifically addresses data of the Vostok ice score whichis a basis for the study of the last million year climate.
1) We consider real or complex functions g ( t ) such that g ( t ) = Z / T − / T G ( f ) e iπft df (2)with G ( f ) = 0 , f / ∈ ( − / T, / T ) . g ( t ) is the Fourier transform (FT) of G ( f ) , and ( − / T, / T ) is a support of G ( f ) . We assume that the boundsof used sets (open, closed..) have no incidence on computations. The historical”sampling formula” addresses the sampling sequence t = { nT, n ∈ Z } . From (2) we have g ( nT ) = Z / T − / T G ( f ) e iπfnT df. (3)To ”develop” g ( t ) as a function of the g ( nT ) is equivalent to develop e iπft asa function of the e iπfnT when f ∈ ( − / T, / T ) . We know that, whateverthe real te iπft = ∞ X n = −∞ sinc (cid:20) π (cid:18) tT − n (cid:19)(cid:21) e iπfnT , − / T < f < / T (4)where sinc x = (sin x ) /x is the ”sinc function” (see appendix 1 about the Fourierseries). In (2) , we replace e iπft by the rhs of (4) , we change the order ofsummations (assuming that this operation is allowed), to obtain2 ( t ) = ∞ X n = −∞ sinc (cid:20) π (cid:18) tT − n (cid:19)(cid:21) g ( nT ) (5)It is the ”sampling formula” of Shannon, Nyquist and others. Conversely, (5)does not work if the ”support ”of G ( f ) is not included in ( − / T, / T ) . /T, the length of this interval, is the maximum which is compatible with the ”sam-pling period” T. It is the ”Nyquist bound”.2) From (3) g ( nT ) is the Fourier series coefficient of the periodic functionequal to G ( f ) on the period ( − / T, / T ) . The Dirichlet theorem states thatthe set of g ( nT ) determines G ( f ) and then g ( t ) (see appendix). We find inEmile Borel ([1], 1897) the following sentences:”Let’s put g ( z ) = Z π − π φ ( x ) e izx dx, and let’s assume that φ ( x ) fulfills the Dirichlet’s conditions. Therefore, if weknow the values of g ( z ) for z = 0 , ± , ± , ... the function φ ( x ) is determined andconsequently the entire function g ( z ) is known without ambiguousness” ( italics come from Borel).3) Formula (5) is verified in the same time as (2) , assuming that a reversalof summations is possible. Using (4) , we see that (5) is true for the larger setof functions as g ( t ) = Z / T − / T G ( f ) e iπft df + X k a k e iπf k t (6)where f k ∈ ( − / T, / T ) (at least when the number of f k is finite). Froma ”mechanical” point of view, (6) isolates ”masses” distributed with a density G ( f ) among ”masses” located at points f k . For acoustic specialists and others,we have a mix of a ”noise” (with spectral density G ( f )) and of a ”signal” madeof pure frequencies f k . Extracting pure frequencies (or concentrations of power) is a problem of”signal theory”. Therefore, decomposition (6) is natural, because it separatesthings of different nature. Curiously, the community prefers both components(”continuous” and ”discrete”) under the shape (2) . In this form, the rhs (2) is no longer an ordinary integral (for instance aRiemann integral). We consider a new set of mathematical items, the ”Diracfunctions” δ ( f ) which verify for f k ∈ ( − / T, / T ) , and by convention e iπf k t = Z / T − / T δ ( f − f k ) e iπft df. (7)In this way, the last term in (6) enters the integral (it is a shoehorn procedure). G ( f ) is replaced by G ( f ) + X k a k δ ( f − f k ) . (8)In this form, G ( f ) is the continuous part of the FT (it is a density) and theresidual is the discrete part (”the spectral lines”).3 .2 The general case
1) Instead of (2) , we consider the more general shape g ( t ) = Z ∆ G ( f ) e iπft df (9)where ∆ is for instance a union of intervals (it is a G ( f ) support). Using thepartition R = ∪ k ∈ Z (cid:0) kT − T , kT + T (cid:1) , (9) becomes g ( t ) = Z / T − / T e iπft H t ( f ) dfH t ( f ) = H t + T ( f ) = X k ∈ Z G (cid:18) f + kT (cid:19) e iπkt/T . (10)Particularly g ( nT ) = Z / T − / T e iπfnT H ( f ) df.T g ( nT ) is the Fourier coefficient of H ( f ) . Conversely H ( f ) = T X k ∈ Z g ( nT ) e − iπfnT . Knowing g ( nT ) , n ∈ Z , is equivalent to knowing H ( f ) = X k ∈ Z G (cid:18) f + kT (cid:19) (11)and knowing of G ( f ) implies g ( t ) . To obtain g ( t ) from H ( f ) it is necessarythat the rhs of (11) contains at most one non-zero term. It is the ”alias-free”condition. Consequently, the length of ∆ (the G ( f ) support) cannot be largerthan 1 /T, but can be arbitrary close of 1 /T, whatever the bounds of any intervalwhich contain ∆ . Knowing ∆ allows to retrieve g ( t ) if, T being given, the alias-free condition is fulfilled.2) For instance, following formulas are fitted to simple ∆a) ∆ = (cid:18) a, a + 1 T (cid:19) ⇐⇒ g ( t ) = e iπt ( a +1 / T ) X k ∈ Z ( − k e − iπkT sinc π (cid:18) tT − k (cid:19) g ( kT ) . This formula comes from the Fourier development of e iπft on ∆ (see appendix).b) ∆ = (cid:18) − k − T , − kT (cid:19) ∪ (cid:18) kT , k + 12 T (cid:19) , k ∈ N ⇐⇒ g ( t ) = 1sin πt/T [ − B ( t ) sin 2 πkt/T + A ( t ) sin π (2 k + 1) t/T ] A ( t ) = X n ∈ Z sinc π (cid:18) tT − n (cid:19) g ( nT )4 ( t ) = X n ∈ Z ( − n sinc π (cid:18) tT − n (cid:19) g ( nT ) . (12)Both intervals of ∆ are folded on disjoined intervals (a motion of ( k/T, (2 k + 1) / T )towards the left and a motion of opposite amplitude ( − (2 k + 1) / T, k/T ) to-wards ( − / T, . Then, the rhs of (11) is reduced to one term . The hypothesis k ∈ N suppresses foldings.
1) The PNS N (Periodic Nonuniform Sampling of order N ) addresses samplingsequences t = { nT + θ k , n ∈ Z , k = 1 , , .., N } (13)where parameters θ k are such that the elements of t are distinct. t is theconcatenation of periodic subsequences t k = { nT + θ k , n ∈ Z } which are interleaved in t. We already find such cases in [3], [4]. From (4) , wededuce the more general formula e iπft = ∞ X n = −∞ e iπαT ( t − θT − n )sinc (cid:20) π (cid:18) t − θT − n (cid:19)(cid:21) e iπf ( nT + θ ) , (14) α − T < f < α + 12
T , t ∈ R available for all real ( α, θ ) . Therefore, we have iat the same time, for any set( α , α , .., α N , θ , θ , .., θ N ) g j ( t ) = Z αj +1 / Tα j − / T G ( f ) e iπft dfg j ( t ) = ∞ X n = −∞ e iπα j T (cid:16) t − θkT − n (cid:17) sinc (cid:20) π (cid:18) t − θ k T − n (cid:19)(cid:21) g j ( nT + θ k ) . (15)It is the Shannon formula (5) in the ”frequency band” δ j = (cid:0) α j − T , α j + T (cid:1) , for the periodic sequence t k (delayed by θ k ) .
2) If we consider the case g ( t ) = Z ∆ G ( f ) e iπft df, ∆ = ∪ Nj =1 (cid:18) α j − T , α j + 12 T (cid:19) (16)we choose the α j so that (cid:26) α j +1 − α j ≥ /Tα j T ∈ β + Z . (17)With the first condition, the δ k are disjoined. The second one allows to re-place the g j ( nT + θ k ) (unknown quantities) by their sum g ( nT + θ k ) (what is5easured) . Therefore, (15) yields N X j =1 g j ( t ) e − iπα j ( t − θ k ) = ∞ X n = −∞ e − iπnβ sinc (cid:20) π (cid:18) t − θ k T − n (cid:19)(cid:21) g ( nT + θ k )(18)It is a linear system N x N which verifies the matricial equation MG t = H t (19)where the columns matrices G t , H t and the square matrix M (independent of t ) are defined by M = (cid:2) e iπα j θ k (cid:3) , G t = (cid:2) g j ( t ) e − iπα j t (cid:3) H t = " ∞ X n = −∞ e − iπnβ sinc (cid:20) π (cid:18) t − θ k T − n (cid:19)(cid:21) g ( nT + θ k ) (20)Provided that det M = 0, which depends on the choice of θ k , equation (19) hasa unique solution. We obtain the value of G t , which leads to the determinationof g ( t ) = N X j =1 g j ( t ) . (21)3) Let’s summarize. We begin by the a periodic sampling and a ”basebandfunction” (∆ = ( − a, a )) . Simple extentions provide the notion of PNS N . It isa mixing of a finite number N of periodic samplings with a common period T. With this model, it is possible to reconstruct a large choice of functions with”spectral support” ∆ included in the union of N intervals of length 1 /T. Let’sreturn to the last example of section 2.2, relaxing hypotheses on ∆ . If∆ = (cid:18) − T , − T (cid:19) ∪ (cid:18) T , T (cid:19) (22)the periodic sampling { nT, n ∈ Z } is not sufficient, because both intervals are”folded” on ( − / T, / T ) by the (symbolic) operations (cid:18) − T , − T (cid:19) + 1 T = (cid:18) T , T (cid:19) − T Nevertheless, for a / ∈ T Z /2 , we prove that g ( t ) = 1sin 2 πa/T [ − A ( t ) sin 2 π ( t − a ) /T + A a ( t ) sin 2 πt/T ] A x ( t ) = ∞ X n = −∞ sinc (cid:20) π (cid:18) t − x T − n (cid:19)(cid:21) g (2 nT + x ) . (23)Therefore, with the PNS2 , t = { nT, nT + a, n ∈ Z } , we cancel the drawbackof a folding which appears in a periodic sampling [4], [9], [7]. The propertyis generalized in: with a well-chosen PNS N, it is possible to overcome N − Nyquist bound and Landau bound
To linearly reconstruct g ( t ) = Z ∆ G ( f ) e iπft df (24)from the sampling sequence t = { t n , n ∈ Z } , is equivalent to find a formula as e iπft = ∞ X n = −∞ a n ( t ) e iπft n , f ∈ ∆ , t ∈ R (25)because, taking some caution, we obtain from (24) g ( t ) = ∞ X n = −∞ a n ( t ) g ( t n ) . (26)The ”Nyquist condition” is summarized by∆ = (cid:18) − µ , µ (cid:19) = ⇒ lim | n |→∞ t n n = µ ≤ µ . (27)It is a necessary condition for reconstruction: if the sampling sequence verifieslim | n |→∞ t n n = µ for some finite µ, a necessary condition for errorless recon-struction of any g ( t ) fulfilling (24) is µ ≤ µ . It is not sufficient (see section5.1), and the main drawback of this condition comes from the shape of ∆ . TheNyquist condition answers problems of analysis about the completion of com-plex exponentials [8], and sufficient conditions can be found in particular cases(Kadec’s th. for example in case of ”jitter”). The condition (27) becomes verystrong when µ is small, and physical models (particularly in communications)prescribe other shapes for ∆ . The ”Landau condition” [10] is summarized by∆ = ∪ k ( a k , b k ) = ⇒ lim | n |→∞ t n n = µ ≤ | ∆ | . (28)For an errorless reconstruction when ∆ is an union of intervals, the ”samplingrate” µ has to be smaller (or equal) to the length of ∆ . Landau condition isweaker than Nyquist condition, and better suited. For instance, we know thatthe frequency spectrum used in communications is divided in bands correspond-ing to telephony, television... This gives a particular interest to the case∆ = ( − b, − a ) ∪ ( a, b ) . In section 2.1 we consider the simplest situation (a Nyquist one) t n = nT, µ = µ = T. In section 2.2, examples refer to Landau condition, and errorless reconstruction.
The Lagrange interpolation sums up points clouds by polynomials [11], [12]. Theextension to entire functions is a part of the complex analysis [13]. Being given7he increasing sequence t = { t n , n ∈ Z } , and g ( t ) in the form (2) for instance,the Lagrange interpolation brings up the validity problem of the equality (insome sense) g ( t ) = X n ∈ Z H ( t )( t − t n ) H ′ ( t n ) g ( t n )with : H ( t ) = Y n ∈ Z (cid:18) − tt n (cid:19) or H ( t ) = t Y n ∈ Z ∗ (cid:18) − tt n (cid:19) (29)following that 0 / ∈ t or t = 0 ∈ t . In baseband, ∆ = ( − / T, / T ) , theNyquist condition has to be verified [13], but it is not sufficient. For instance(29) is true for t= Z , T = 1. It is the standard formula knowing thatsin πx = πx Y n =0 (cid:16) − xn (cid:17) . (30)Formula (29) is not verified for t= Z - { } , even if G ( f ) cancels at the neigh-bourood of ± / . In this case, it is easy tofind alternative formulas. The problem is linked to the behavior of H ( t ) at theinfinite points. The difference comes fromlim | x |→∞ | sin x | = 1 vs lim | x |→∞ sin xx = 0where t = 0 ( t = Z ) or not ( t= Z - { } ). In real situations, devices compute areduction e g ( t ) of (29) e g ( t ) = [ t ]+ N X n =[ t ] − N e H ( t )( t − t n ) e H ′ ( t n ) g ( t n ) (31)because of limited capacities of devices and the future of t n is not very wellknown ([ t ] is close to t so that it is the neighbourhood of t which is takeninto account). This type of approximation enters the class of operations named” apodization ”. e H ( t ) can take many shapes. The simplest is to retain the t n in some interval around t [14], [15] . Another widely used technique is the ” resampling” . Let’s assume thatformula (5) is available and that we know g ( t k ) in a neighbourhood of t. For N large enough, we will often have g ( t k ) ≈ [ t ]+ N X n =[ t ] − N sinc (cid:20) π (cid:18) t k T − n (cid:19)(cid:21) g ( nT ) , | k | ≤ N. This (approximative) linear system ( N + 1)x( N + 1) provides the g ( nT ) whichallows to determine estimations of g ( t ) whatever t .
1) We consider the (increasing) sampling sequence t = { t n , n ∈ Z } , and a function g ( t ) verifying (16) , (17). We have | ∆ | = N/T, and we take t N = { nT + θ k , n ∈ Z , k = 1 , , .., N } , θ k = t N + k . (32)818) , (19) are true whatever integer N . Let’s assume that t n +1 − t n < T /N, sothat the sequence t N shows a natural order due to the inequality t N + N − T < t N +1 < ... < t N + N < t N +1 + T. For a given t, we choose N = N ( t ) such that (for even N ) t N +1 < ... < t N + N/ ≤ t ≤ t N +1+ N/ < ... < t N + N . In such a situation, when N is large enough, we can replace (19) , (20) by theapproximations M e G t = e H t M = (cid:2) e iπα j θ k (cid:3) , e G t = (cid:2)e g j ( t ) e − iπα j t (cid:3) , e H t = (cid:20) sinc (cid:20) π (cid:18) t − θ k T (cid:19)(cid:21) g ( θ k ) (cid:21) . (33)This means that we suppress in H t the terms n = 0 , which have abscissas thefarthest away from t. These terms are not measured. The estimation of g ( t )is obtained from N measures g ( t k ) , at points t j which are the closest to t. Itis equivalent to an apodization of an errorless formula, not taking into accountunknown quantities (the g ( nT + θ k ) which are not measured for n = 0) .
2) For each (
N, N ) , the sequence t N defined in (32) verifies the Landaucondition, even if lim n −→±∞ t n /n does not exist. The quality of estimations e G t of G t depends on the ”distance” between e H t and H t . Obviously, the technique iswell suited when t is chosen in the interval ( t , t N ) , and far from the nT + t k , n =0 , where we can expect that the omitted terms are negligible.When G ( f ) is regular enough and N large enough, so that the mean theoremof integral can be involved, we have, from (15) g k ( t ) ≈ | ∆ | N e iπα k t G ( α k ) . (34)Approximately, the ”spectral content ” of g ( t ) at α k is close to | ∆ | N G ( α k ) . Equivalently, the filtering of g ( t ) in the band δ k is quasi-monochromatic, itsamplitude and its phase being defined by | ∆ | N G ( α k ) . Z= { Z ( t ) , t ∈ R } is a (wide sense) stationary random process when K ( τ ) = E [ Z ( t ) Z ∗ ( t − τ )] = Z ∆ e iπfτ s ( f ) df (35)where E[ .. ] is the mathematical expectation (ensemble mean), a ∗ is the complexconjugate of a [16], [17]. The integral has the common sense (a Riemann inte-gral) when s ( f ) is regular enough. Then, s ( f ) is real and positive ( ≥
0) and issymmetric for real Z. Elsewhere, s ( f ) is explained by (6) , (7) . It is a sum oftwo terms, the first one is an ordinary function (but positive) the other one is asum of complex exponentials with positive coefficients. s ( f ) defines the ”powerspectrum” of Z. ∆ is a support of the spectrum ( s ( f ) = 0 outside ∆) . t n , Z ( t n )) we have to find a linear reconstruction of Z ( t ) , i.e a formula as Z ( t ) = X n ∈ Z a n ( t ) Z ( t n ) (36)where a n ( t ) depends only of ∆ (and not of values of s ( f ) on ∆) . In the mean-square sense, equality (35) verifies E (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) Z ( t ) − X n ∈ Z a n ( t ) Z ( t n ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) = 0 ⇔ Z ∆ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) e iπfτ − X n ∈ Z a n ( t ) e iπft n (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) s ( f ) df = 0This means that we have to develop e iπft on ∆ as function of e iπft n , n ∈ Z .With mathematical nuances about differences between kinds of convergences,integrals..., the sampling problem remains the same that for functions like (2)which have Fourier transforms on some set ∆ . We have to look for a developmentof e iπft on ∆ as e iπft = X n ∈ Z a n ( t ) e iπft n , f ∈ ∆ , t ∈ R (37)and to replace e iπft by Z ( t ) , for t, t j , j ∈ Z , in (35) . We obtain, from section3 M e G t = e H t M = (cid:2) e iπα j t k (cid:3) , e G t = h e Z j ( t ) e − iπα j t ie H t = (cid:20) sinc (cid:20) π (cid:18) t − t j T (cid:19)(cid:21) Z ( t j ) (cid:21) (38)with the approximation e Z ( t ) = N X k =1 e Z k ( t ) (39) Z k ( t ) reults from a linear filtering of Z on δ k . The Z k are uncorrelated (the δ k are disjoined) and K k ( τ ) = E [ Z k ( t ) Z ∗ k ( t − τ )] = Z δ k e iπfτ s ( f ) df (40)In (34) , the power spectrum s ( f ) (for Z ) has the sense attributed to G ( f ) (for g ( t )) in section 2.1. s ( f ) is the sum of a continuous part and a discrete part.This means that we observe accumulations of power in some places, which arepreserved when N increases. 10 Power spectra measurements
From (34), the power spectrum s ( f ) of Z, is the inverse Fourier transform ofthe correlation function K ( t ) , classically estimated from e K ( nT ) = 1 N − n N − n X k =1 Z ( kT ) Z (( k + n ) T ) , n = 0 , .., N − . The estimation worsens when n increases.Generally, people consider that the accuracy in the research of high frequencies(v.s low) is very dependent on the accuracy in the estimation for low values (v.shigh) of n. If true, we have to dispose of samples in a large interval of time ifwe want low frequencies, and samples very close for high frequencies. Moreover,in the case of irregular sampling, (40) will be achevied when the Z ( kT ) will bethemselves estimated from the Z ( t k ).When integer N is large enough, s ( f ) can be estimated directly from PNS.To simplify notations, we take the particular case∆ = ( − b, b ) , δ k = (cid:18) − b + 2 bN ( k − , − b + 2 bN k (cid:19) , k = 1 , , .., N. We have | δ k | = 2 b/N . Let’s admit that (37) provides good estimations ofelements of G t Z k ( t ) e − iπα k t , α k = − b + bN (2 k − G t belongs to the frequential set ( − b/N, b/N ) . This means that Z j ( t ) e − iπα j t is a demodulation of Z j ( t ) in baseband. When b/N is smallenough, Z k ( t ) appears as the monochromatic wave a ( t ) e iπα k t with an ampli-tude | a ( t ) | which varies slowly. The total power P k of Z k is (it is a definition) P k = E h | Z k ( t ) | i = E h | G k ( t ) | i which is estimated measuring the mean of | G k ( t ) | in (0 , t N ) . Vostok ice cores are columns of ice of near four kilometres length coming from thetop to the bottom of a lake of southern antarctica (near the magnetic south pole).The top (v.s the bottom) corresponds to recent (v.s around 400000 years old)snow. The study of samples collected along the ice core gives many indications,the date, temperature, composition and variations of elements as CO , CH , D (deuterium), O (isotope of oxygen)... Figure 1 shows the evolution of CO concentration on a period of 400000 years (the past at the right). Data are in* and interpolation in - - -. Figure 2 shows the corresponding power spectrumon (cid:0) , . − (cid:1) in yr − , with N = 282 and | δ k | = 7 . − yr − . The originalestimation is in inset [19]. Improvements on the accuracy of estimations areclear, by the sharpeness of spectral rays and by the appearance of low frequencycomponents. Figure 3 shows components 3, 9, 10,... ( Z ( t ) , Z ( t ) .. ) . In figure 411e see an addition of the first three rays in low frequencies. Figure 5 representsthe power spectrum of CH . It is very similar to the former one (figure 2), withthe very strong ray at 10 − yr − (related to the 100000 yr problem [21]).The Vostok ice core is a very telling example of irregular sampling. Dataare not in arithmetic series, due to varying conditions of precipitation and me-chanical loads [22], [23]. Estimations of power spectra have to be compared toresults in [20], where are detailed the best methods by the best specialists. In the Vostok ice core, we have taken ∆ = (cid:0) − − , − (cid:1) and the δ k are adja-cent. Nyquist and Landau conditions are equivalent. It is no longer true when∆ is piecewise. Figure 6 addresses a Gaussian process in baseband modulatedin ( − / , − ∪ (2 , /
2) i.e multiplied by cos (9 πt/ . Sampling times are in theform t n = n (1 − ε ) + A n , < ε < / A n are uniformly distributed in ( − / , − / ∪ (1 / , / . The esti-mation at time t is done using 10 samples before t and 10 after. On figure 6, wefind samples (*) and a perfect reconstruction (at the thickness of the line). Wehave lim n →∞ t n n = 1 − ε < n →∞ t n /n = 5) . This shows that the Nyquist condition is often pessimistic.
This paper addresses irregular sampling of ordinary functions g ( t ) and randomprocesses Z with bounded support spectrum. We explain why PNS (PeriodicNonuniform Samplings) are the right framework when we have at our disposala large enough number of samples (but a finite number). To construct errorlessformulas, we begin by the periodic sampling, suited to an arbitrary interval.This leads to fit PNS N to piecewise power spectra, verifying the Landau con-dition ( N is the number of samples). The reconstruction of g ( t ) or Z is donethrough a N x N matricial system which solutions are the N components corre-sponding to pieces of power spectra. The mathematical context does not exceedthe elementary framework of Fourier analysis (one-dimensional series and inte-grals).Considerations about Vostok ice cores illustrate theoretical results. Theyprovide data which give strong informations about climate in the last 4.10 years. Physical-chemical studies have allowed to measure the concentration inC0 , CH , isotopes of oxygene... of atmosphere trapped in ice at irregular times.Sampling methods explained here allow to reconstitute these quantities andtheir evolution on a long period of the past. The proposed method, thanks toband decomposition in small intervals, shows improvements in the number andthe thickness of spectral rays. Particularly, low frequency components appear,which are hidden in usual computations.12
1) Let’s h ( f ) such that h ( f ) = h ( f + f )for ” h ( f ) is periodic with period f ” . Fourier cefficients c n are defined as c n = 1 T Z f / − f / h ( f ) e − iπnf/f df. We have the following property (provided some conditions about the h ( f ) reg-ularity and the convergence definition) . If e h ( f ) = ∞ X n = −∞ c n e iπnf/f , then h ( f ) = e h ( f ) . (42) e h ( f ) is the Fourier Series Development (FSD) of h ( f ) .
2) This definition of FSD is dangerous. Firstly, a function h ( f ) being given,we may look for a development of h ( f ) verifying (42) on some set ∆ for non-periodic h ( f ).For instance, if h ( f ) = f , f ∈ R , ∆ = ( a, a + b ) , we have c n = 1 b Z a + ba f e − iπnf/b df = − b iπn e − iπna/b , n = 0 e h ( f ) = a + b − bπ ∞ X n =1 n sin (cid:20) πnb ( f − a ) (cid:21) . For n ∈ Z and except at the bounds of intervals : e h ( f ) = (cid:26) h ( f ) , f ∈ ( a, a + b ) h ( f ) − nb, f ∈ ( a + nb, a + ( n + 1) b )Figure 7 explains the result. We have found a trigonometric development of h ( f ) on ∆ (it is e h ( f )) . e h ( f ) is periodic, as a sinus, and then is different from h ( f ) outside ∆ . Changing the value of a and/or b yields a FSD different onanother set ∆ .
3) Then, a function defined on the whole R , provides a different FSD fol-lowing the set ∆ where the development takes place. It is true when h ( f )is periodic with period f . Obviously, h ( f ) can be developped on intervals oflength different of f . In the beginning of this paper, we prove the standardsampling formula (Shannon, Nyquist,...) using FSD when∆ = ( − / T, / T ) and h ( f ) = e iπft .h ( f ) is a periodic function of f ∈ R with period f = 1 /t. Formula ( .. ) isthe DSF of h ( f ) on ∆ . It is natural that coefficients depend on parameter t, but t it is not the used variable. The obtained DSF provides h ( f ) on ∆ butnot elsewhere (except the particular case t = T ) . Figure 8 illustrates this veryfundamental situation. 13 eferences [1] E. Borel,
Sur l’interpolation,
CRAS
124 ( ) La formule d’´echantillonnage et A. L. Cauchy , Traitement dusignal (2) (2004) 129-140.[3] J. L. Yen, On Nonuniform Sampling of Bandwidth-Limited Signals,
IRETrans. on Circ. Th.,
CT-3 (12) (1956) 251-257.[4] A. Kohlenberg,
Exact interpolation of band-limited functions,
J. AppliedPhysics, (12) (1953) 1432-1436.[5] Y-P. Lin, P. P. Vaidyanathan, Periodically Nonuniform Sampling of Band-pass Signals,
IEEE Trans. on Signal and Systems-II (3) (1998) 340-351.[6] R. Venkataramani, Perfect Reconstruction Formulas and Bounds on Alias-ing Error in Sub-Nyquist Nonuniform Sampling of Multibands Signals,
IEEE Trans. on Inf. Th. (6) (2000) 2173-2183.[7] B. Lacaze, Equivalent Circuits for the PNS2 Sampling Scheme,
IEEE Cir-cuits and Systems 1 (11) (2010) 2904-2914.[8] N. Levinson, Gap and density theorems,
AMS (1940).[9] A. J. Jerri, The Shannon sampling theorem. Its various extensions andapplications. A tutorial review,
Proc. IEEE (11) (1977) 1565-1596.[10] H. J. Landau, Sampling, data transmission, and the Nyquist rate,
Proc . IEEE (1967) 1701-1706.[11] J. L. Lagrange, , 1795.[12] E. Waring, Problems concerning interpolations,
Phil. Trans. Royal Soc. ofLondon, (1779) 59-67.[13] B. Ja. Levin, Distribution of Zeros of Entire Functions,
American Mathe-matical Society (1966).[14] B. Lacaze,
Reconstruction formula for irregular sampling,
Samp. Th. Sign.Image Proc. (1) (2005) 33-43.[15] B. Lacaze, The Ghost Sampling Sequence Method,
Samp. Th. Sign. ImageProc. (1) (2009) 13-21.[16] A. Papoulis, Signal Analysis,
Mc-Graw Hill (1977).[17] H. Cram´er, M. R. Leadbetter,
Stationary and related stochastic processes,
Wiley, New-York (1966).[18] B. Lacaze,
A Theoretical Exposition of Stationary Processes Sampling,
Samp. Th. Sign. Image Proc. (3) (2005) 201-230.[19] J. R. Petit et all, Climate and atmospheric history of the past 420,000 yearsfrom the Vostok ice core, Antartica,
Nature (6735) (1999) 429-436.1420] P. Babu, P. Stoica,
Spectral analysis of nonuniformly sampled data - areview,
Digital Signal Processing (2010) 359-378.[21] J. A. Rial, J. Oh, E. Reischmann, Synchronization of the climate systemto eccentricity forcing and the 100000-year problem,
Nature Geoscience (2013) 289-293.[22] B. Lacaze, Reconstruction from Periodic Nonstationary Sampling (PNS)without Resampling . Th. Sign. Image Proc. (2017) 1-21.[23] D. Bonacci, B. Lacaze, New CO concentration predictions and spectralestimation applied to the Vostok ice core , IEEE Trans. on geoscience andremote sensing, (1) (2018) 145-151.15(1) (2018) 145-151.15