Frequency Transforms based on Nonuniform
Sampling — Basic Stochastic Properties
Frida Eng
,
Fredrik Gustafsson
Control & Communication
Department of Electrical Engineering
Link¨
opings universitet
, SE-581 83 Link¨
oping, Sweden
WWW:
http://www.control.isy.liu.se
E-mail:
frida@isy.liu.se
,
fredrik@isy.liu.se
1st October 2004
AUTOMATIC CONTROL
COM
MUNICATION SYSTEMS
LINKÖPING
Report no.:
LiTH-ISY-R-2631
Submitted to ICASSP 2005
Abstract
We investigate two natural linear estimates of a signal’s continuous time Fourier transform (FT) based on non-uniform samples, where sam-ple value and samsam-ple time come in pairs where at least the sampling time is a stochastic process. Such stochastic sampling occurs in several applications, such as queue theory and event based sampling. Analyti-cal expressions of mean and covariance are derived for the FT estimate. Intended use of these are to give uncertainty bounds on transforms and spectral estimates, but they are also useful in signal estimation and system identification.
Keywords: nonuniform sampling, frequency analysis, stochastic window
Frequency Transforms based on Nonuniform
Sampling — Basic Stochastic Properties
Frida Eng
∗and Fredrik Gustafsson
Departement of Electrical Engineering
Link¨
opings universitet, SE-58183 Link¨
oping, Sweden
Email: frida.eng@isy.liu.se
September 30, 2004
Abstract
We investigate two natural linear estimates of a signal’s continuous time Fourier transform (FT) based on non-uniform samples, where sam-ple value and samsam-ple time come in pairs where at least the sampling time is a stochastic process. Such stochastic sampling occurs in several applications, such as queue theory and event based sampling. Analyti-cal expressions of mean and covariance are derived for the FT estimate. Intended use of these are to give uncertainty bounds on transforms and spectral estimates, but they are also useful in signal estimation and system identification.
1
Introduction
The signal processing theory is to a large extent focused on uniform sampling, where the measurements
yk = s(kT ) + ek
are observations of a signal s(t) at each sample interval T , corrupted with mea-surement noise ek. Non-uniform sampling concerns the more general case where yk = s(tk) + ek. The time instants tk may for instance be generated by an ex-ternal process, in which case we have event based sampling. This is the case for example when measuring the angle of rotating axles. The measurement may also be commanded by a non-real time system, where a certain sampling jitter is added. Further, arrivals and departures in queue processes have a stochastic time (but deterministic signal change of plus or minus one unit).
Nonuniform sampling, both deterministic and stochastic, is described in Bilinskis and Mikelsons (1992), Marvasti (2001). The literature describes among other things how the Fourier Transform (FT) can be approximated when the
sampling times are deterministic. For instance, the Lomb-Scargle periodogram Lomb (1976), Scargle (1982) is introduced as a least square fitting at each fre-quency. The value of the periodogram at frequency f is
Y0(f ) = arg min a N X k=1 |yk− ae−2πf tk|2= 1 N N X k=1 yke−2πf tk. (1) This coincides, apart from a scaling of 1/N , with the Dirichlet transform, which is suggested as the accurate way of performing frequency transform of non-uniform time samples in Bilinskis and Mikelsons (1992). Another alternative is a Riemann approximation of the Fourier integral
Y1(f ) = N X k=1 yk(tk− tk−1)e−2πf tk≈ Z y(t)e2πf tkdt. (2) This work is concerned with additive random sampling, i.e., sampling times are formed as
tk = tk−1+ τk, (3)
with τkbeing realizations from a stationary stochastic process with known prop-erties, e.g., E[τk] = µT. The problem under consideration is how this kind of stochastic sampling time process tk affects the two FT approximations in (1) and (2).
The outline is as follows: Section 2 introduces notation. Section 3 derives the main results. In short, the mean and covariance will turn out to obey the relations
E(Yi(f )) = S ∗ E(Wi)(f )
Cov(Yi(f1), Yi(f2)) = S ∗ ∗ Cov(Wi)(f1, f2),
(4) where S is the FT of the signal s(t), and ∗ and ∗∗ we denote single and double convolution defined as f ∗ g(x) = Z f (ξ)g(x − ξ)dξ, f ∗ ∗h(x, y) = Z Z f (ξ)h(x − ξ, y − ψ)f∗(ψ)dξdψ.
The stochastic data window Withus describes both the FT approximation and the sampling process, and W0coincides with the usual leakage window for uni-form sampling. Section 4 extends the results to the case of additive measurement noise, and Section 5 illustrates with an example which also motivates a Gaussian distribution in Eq. (4).
The problem that is treated in this work was introduced in Gunnarsson (2003), where a lot of different transform approximations were investigated. In this work two transforms will be further studied, the Dirichlet transform and the
Riemann transform. They both coincide with the DFT for uniformly sampled data.
In Gunnarsson et al. (2004), stochastic sampling was introduced into the problem and initial investigations of the stochastic window, W , were made.
The frequency transform of nonuniformly sampled data can also be used for parameter estimation, especially in the case of jitter sampling with limited but unknown jitter, yk = s(kT + τk). This problem is investigated in Eng and Gustafsson (2004).
2
Frequency Transform of Nonuniform Samples
A signal with the Fourier transform representation s(t) =
Z ∞ −∞
S(φ)ei2πφtdφ, (5)
is sampled at nonuniform times giving signal samples, yk = s(tk), k = 1, . . . , N . The frequency transform is given as a linear transform of the measurements
Y (f ) = 1 N N X k=1 ykwk(f ) (6) = 1 N N X k=1 Z ∞ −∞ S(φ)ei2πφtkdφw k(f ) (7) = Z ∞ −∞ S(φ) 1 N N X k=1 wk(f )ei2πφtk | {z } W(f,φ) dφ, (8)
which gives the result as a convolution of the true transform with a stochastic window, W (·). The weights, wk, can be varied to produce different kinds of transform approximations. In this work we consider two types with,
wk(f ) = τk µT n e−i2πf tk (9) = e −i2πf tk, n = 0, Dirichlet τk µTe −i2πf tk, n = 1, Riemann , (10)
where n ∈ {0, 1} is used to chose between the transform types. The Dirichlet transformis produced when n = 0, and n = 1 gives the Riemann transform, both seen before in e.g., Gunnarsson et al. (2004), Bilinskis and Mikelsons (1992), Persson (2002). We will study the properties of the window, W , in order to gain knowledge of Y . The above definitions of the weights wk imply that W (f, φ) = W (f − φ), which results in a nice convolution formulation.
3
Expected Value and Covariance
The expected value and covariance of the transform Y can be found as E[Y (f )] = Z S(φ) E[W (f − φ)]dφ, (11) RY(f, g) = Cov(Y (f ), Y (g)) = Z Z S(φ) Cov(W (f − φ), W (g − ϕ))S(ϕ)∗dφdϕ, (12) which shows that the first moments of the stochastic window, W , is important. The expected value and covariance of the stochastic window, W , are both prod-ucts of a sampling part and a transform type part. These two parts correspond to two types of leakage.
Lemma 1 (Expected value of frequency window) Let the setting in Eqs. (3), (6) and (9) hold. Then, the stochastic window has a mean
E[W (f )] = Γtp(f )ΓN(f ),
where thesampling leakage ΓN(f ) = N1 P γ(f)k−1and theinterpolation leakage Γtp(f ) = γn(f ).
Proof: γn is defined as:
γn(f ) = ET[ τk µT n e−i2πf τk], (13) γ(f ) = γ0(f ) = ET[e−i2πf τk]. (14) The function γ(f ) is the moment generating function (MGF) for the distribution of τk, and the exact expression can be found in textbooks on stochastic variables, for various distributions. The calculations are straightforward and given in
Appendix A.
Lemma 2 (Covariance of frequency window) Given the
background in Lemma 1, the missing part for the covariance is given as E[W (f )W∗(g)] = ΛT
tp(f, g)ΛN(f, g). The sampling covariance is given as
ΛN(f, g) = ΣN(f, f − g) 1 NΓN(f − g) ΣN(−g, g − f )
and, the interpolation covariance as
Λtp(f, g) = Γtp(f )Γtp(f − g) γ2n(f − g) Γtp(−g)Γtp(−g, g − f ) .
Proof: The calculations are straightforward and involve the following block: ΣN(f, g) = 1 N2 N X k=1 k−1 X l=1 γ(f )k−l−1γ(g)l−1 = 1 N(γ(f )−γ(g))[ΓN(f ) − ΓN(g)] , f 6= g (N −1)γ(f )N−N γ(f )N −1+1 N2(1−γ(f ))2 , f = g 6= 0 (N −1) 2N , f = g = 0 (15)
The explicit calculations can be found in Appendix A.
3.1
Asymptotic analysis
¿From the simple fact that |γ(f )| ≤ 1 and γ(f ) = 1 only if f = 0, properties of the transforms can be derived for when N → ∞.
Lemma 3 The transform is asymptotically unbiased, i.e., E[Y (f )] → S(f )Γtp(0) = S(f ), N → ∞ and its covariance asymptotically vanishes, i.e.,
RY(f, g) → 0, N → ∞.
Proof: The proof exploits the limit features of ΓN, and is sketched in
Ap-pendix A.
4
Measurement Noise
Usually a measured signal complies of a deterministic wanted part corrupted with measurement noise, i.e.,
yk= sk+ vk, (16)
with E[vk] = 0 and E[vkvl] = σ2δk−l. This gives Y (f ) = 1 N X skwk(f ) | {z } S∗W (f ) + 1 N X vkwk(f ) | {z } V(f ) .
Lemma 4 (Additive measurement noise) If the measurements fulfill Eq. (16) the noise affects the transform distribution as
E[Y (f )] = S ∗ E[W (f )]
Cov(Y (f ), Y (g)) = S ∗ ∗ Cov(W )(f, g) +σ 2
Nγ2n(f − g)ΓN(f − g). Thus, the measurement noise only affects the covariance of the transform, with a term vanishing as the number of measurements increase.
Proof: The mean value is easily seen E[Y ] = E[S ∗ W ] + 1 N X E[vkwk(f )] = E[S ∗ W ] + 1 N X EE[vk|τk] wk(f ) = S ∗ E[W ].
The covariance consists of four terms
Cov(Y ) = E[(Y (f ) − E[Y (f )])(Y (g) − E[Y (g)])∗] = Cov(S ∗ W (f ), S ∗ W (g))+
+ E[(S ∗ W (f ) − E[S ∗ W (f )])] E[V (g)∗]+ + E[V (f )] E[(S ∗ W (g) − E[S ∗ W (g)])∗]+ + E[V (f )V (g)∗].
Since E[V ] = 0, the only unknown term is E[V (f )V (g)∗] = 1 N2 X X E[vkvlwk(f )wl(g)∗] = σ 2 N2 X k E[wk(f )wk(g)∗], = σ 2 N2 X k γ2n(f − g)γ(f − g)k−1,
which, together with the definition of ΓN in Lemma 1, concludes the proof.
5
Numerical illustrations
To be able to give confidence bounds on the transform approximation, we need its distribution. It basically follows from the central limit theorem for dependent stochastic variables that this is the case for (6). This will be further validated in the example. Furthermore, for spectral analysis, the distribution of the trans-form magnitude is needed. Once this is known, confidence bounds on the pe-riodogram |Y (f )|2 and smoothed spectral estimates based on the periodogram are readily computed.
We start by noting the following fact: If z comes from a complex Gaussian distribution with mean µz, and variance σz2, then |z| is Ricean distributed with known mean µ|z| and variance R|z|. The stochastic variable r = |z − µz| is Rayleigh distributed with known mean µr and variance Rr. These values are related as
|µz| ≤ µ|z|, max(Rr− 2µr|µz|, 0) ≤ R|z| ≤ Rz.
0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 0 1 2 3 4 5 6 7 8 9 R |z| lower bound upper bound R r Ri 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1 2 3 4 5 6 7 8 9 10 |µz| µ|z| µrRi 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 0 1 2 3 4 5 6 7 8 9 10 R|z| lower bound upper bound RrRi 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1 2 3 4 5 6 7 8 9 |µz| µ|z| µrRi
Figure 1: The results from Example 1 are shown together with the limits in Eq. (17). The plots from left to right, variance and mean for the Dirichlet transform and variance and mean for the Riemann transform. Comparison of mean value, µ|z| and variance, R|z|, of the transform amplitude, |Y |, with theoretical values, µRi
r and RRir , from the Rice distribution. The upper and lower limits were given in Eq. (17), but it can be noted that for the mean value, µ|z|, the lower limit is |µz| = | E[Y (f )]| from Lemma 1.
The achieved value for Y will vary depending on the realization of the sam-pling times, and the distribution will be studied using Monte Carlo runs for a specific signal.
Example 1 For simplicity the transform of a single sinusoid is studied, i.e., s(t) = sin(2πf0t). For each value of f and f0, Y (f ) = S ∗ W (f ) will be a complex stochastic variable.
4000 Monte Carlo evaluations of the transform have been done, with N = 20 andτk ∈ U (τL, τH). For a specific value of f0 the following was performed:
1. YM(f ) is M instantiations of the transform Y (f ).
2. NumericallyRz(f ) = Var(YM(f )) and µz(f ) = mean(YM(f )) is found. 3. The theoretical mean values µRa
r (f ) = p
Rz(f )π/2 and µRir (f ) for the Rayleigh and Rice distribution corresponding to the parameters above are computed.
4. The theoretical variances RRa
r (f ) = Rz(f )(1 − π/4) and RRir (f ) for the Rayleigh and Rice distribution are calculated.
5. Numerically R|z|(f ) = Var(|YM(f )|) and µ|z|(f ) = mean(|YM(f )|) are produced.
The calculated values can be used for numerical comparison.
The results from Example 1 are shown in Figure 1, together with the limits given in Eq. (17), for the two transform types. The plots illustrate the following: • It can be concluded that the magnitude of the transforms have their first two moments equal to a those of Ricean distribution, which validates the Gaussian assumption and that the conditions in the central limit theorem are fulfilled in this example.
• It is also interesting to note that the upper and lower limits are quite conservative, which demonstrates the importance of knowing the true dis-tribution for Y . Using the mean value and variance from Section 3 for Y , the same quantities for |Y | can be found from the above discussion. • The magnitude of expected value, |µz|, and the expected value of the
magnitude, µ|z|, are identical close to the true frequency, f0. This indicates that |Y | is Rayleigh distributed for f = f0.
6
Conclusions
This paper presented analysis of the stochastic properties of two frequency trans-forms of nonuniformly sampled signals. Analytic expressions for the mean value and covariance were given, based upon the distribution of the inter sampling times. These are important to know in practice, to be able to give confidence
bounds on the transform approximation, and to plot confidence limits on the periodogram and smoothed versions thereof. Future work will utilize these re-sults for signal estimation and system identification in the case of non-uniform stochastic sampling.
References
Ivar Bilinskis and Arnolds Mikelsons. Randomized Signal Processing. Prentice Hall, London, 1992.
Frida Eng and Fredrik Gustafsson. System identification using measurements subject to stochastic time jitter. Technical report, Department of Electrical Engineering, Link¨opings Universitet, Sep 2004.
Frida Gunnarsson. On modeling and control of network queue dynamics. Licen-ciate Thesis, No. 1048, Dep. of Electrical Engineering, Link¨opings Universitet, 2003.
Frida Gunnarsson, Fredrik Gustafson, and Fredrik Gunnarsson. Frequency anal-ysis using nonuniform sampling with applications to adaptive queue man-agement. In Acoustics, Speech, and Signal Processing., IEEE International Conference on, May 2004.
N R Lomb. Least-squares frequency analysis of unequally spaced data. Astro-physics and Space Science, February 1976.
Farokh Marvasti, editor. Nonuniform Sampling: Theory and Practice. Kluwer Academic Publishers, 2001.
Niclas Persson. Event Based Sampling with Application to Spectral Estima-tion. Licentite Thesis No. 981, Dep. of Electrical Engineering, Link¨opings Universitet, 2002.
Jeffrey D Scargle. Studies in astronomical time series analysis II: Statistical as-pects of spectral analysis of unevenly spaced data. The Astrophysical Journal, December 1982.
A
Proofs of the Stated Lemmas
This section states the proofs of the previous Lemmas. The notation and the signal model is the same as before.
A.1
Proof of Lemma 1
The inter sampling times τk are independent and identically distributed, which means that the expected value of the stochastic window can be calculated as follows: E[W (f )] = 1 N N X k=1 E[ τk µT n e−i2πf tk] = 1 N N X k=1 E[ τk µT n k Y l=1 e−i2πf τl] = 1 N N X k=1 γn(f )γ(f )k−1 = γn(f )ΓN(f ),
which gives the result of the Lemma by identification of terms.
A.2
Proof of Lemma 2
The covariance for W is complicated to read but simple to produce. E[W (f )W∗(g)] = 1 N2 N X k=1 N X l=1 E[ τkτl µ2 T n e−i2πf tkei2πgtl] = N X k=1 k−1 X l=1 · · · + N X k=1 k X l=k · · · + N X k=1 N X l=k+1 · · · = Φ1(f, g) + Φ2(f, g) + Φ3(f, g).
Each term is treated separately for convenience,
Φ1(f, g) = 1 N2 N X k=1 k−1 X l=1 E[ τkτl µ2 T n k Y l+1 e−i2πf τj l Y 1 e−i2π(f −g)τj].
The expected value is taken given that the inter sampling times are independent,
Φ1(f, g) = 1 N2 N X k=1 k−1 X l=1 γn(f )γ(f )k−1−lγn(f − g)γ(f − g)l−1, = γn(f )γn(f − g)ΣN(f, f − g)
The second term is the sum when k and l are equal, Φ2(f, g) = 1 N2 N X k=1 E[τk2n k Y j=1 e−i2π(f −g)τj], = 1 N2 N X k=1 γ2n(f − g)γ(f − g)k−1 = 1 Nγ2n(f − g)ΓN(f − g) and the last term is given by symmetry,
Φ3(f, g) = 1 N2 N X k=1 N X l=k+1 E[ τkτl µ2 T n l Y k+1 ei2πgτj k Y 1 ei2π(f −g)τj] = 1 N2 N X l=1 l−1 X k=1 E[ τkτl µ2 T n l Y k+1 ei2πgτj k Y 1 ei2π(f −g)τj] = Φ1(−g, −f )
Identification of terms proves the Lemma.
A.3
Proof of Lemma 3
Since |γ(f )| < 1, f 6= 0 and γ(0) = 1, and ΓN(f ) = 1 N N X k=1 γ(f )k−1 = ( 1 N 1−γ(f )N 1−γ(f ) , f 6= 0, 1, f = 0,
and integral such as
Z H(f )ΓN(f )df → H(0), when N → ∞. Similarly, Z Z H(f, g)ΣN(f, g)df dg → 1 2H(0, 0),
based on the definition of ΣN(f, g) in Eq. (15). Including this into Eqs. (11) and (12) gives the result of the Lemma:
E[ ˆY (f )] = Z Y (φ)γn(f − φ)ΓN(f − φ)dφ → Y (f )γn(0) = Y (f ). RY(f, g) = Z Z Y (φ)Γ1 tp(f − φ, g − ϕ)ΣN(f − φ, f − φ − g + ϕ)Y (ϕ)∗dφdϕ+ + 1 N Z Z Y (φ)Γ2tp(f − φ, g − ϕ)ΓN(f − φ − g + ϕ)Y (ϕ)∗dφdϕ+ + Z Z Y (φ)Γ3tp(f − φ, g − ϕ)ΣN(ϕ − g, g − ϕ − f + φ)Y (ϕ)∗dφdϕ+ − E[ ˆY (f )]E[ ˆY (g)]∗ → 1 2Y (f )Γ 1 tp(0, 0)Y (g)∗+ 0 + 1 2Y (f )Γ 3 tp(0, 0)Y (g)∗− Y (f )Y (g)∗ ≡ 0, since Γ1tp(0, 0) = Γ3tp(0, 0) = E[ τk µT n
] = 1. The last equality, of course, only holds when n ∈ {0, 1}, which is exactly the set we restricted the transform component, n, to, in the definition.
