Table of Contents
ISRN Signal Processing
Volume 2012 (2012), Article ID 643563, 10 pages
http://dx.doi.org/10.5402/2012/643563
Research Article

Observability of Spectral Components beyond Nyquist Limit in Nonuniformly Sampled Signals

Institute of Electronics and Photonics, Slovak University of Technology in Bratislava, Ilkovičova 3, 81219 Bratislava, Slovakia

Received 27 March 2012; Accepted 3 May 2012

Academic Editors: A. M. Peinado and G. A. Tsihrintzis

Copyright © 2012 Jozef Púčik et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract

Identification of a signal component with the frequency exceeding the Nyquist limit is a challenging problem in signal theory as well as in some specific applications areas like astronomy and biosciences. A consequence of the well-known sampling theorem for a uniformly sampled signal is that the spectral component above the Nyquist limit is aliased into lower frequency range, making a distinction between true and aliased components impossible. The nonuniform sampling, however, offers a possibility to reduce aliased components and uncover information above the Nyquist limit. In this paper, we provide a theoretical analysis of the aliased components reduction in the nonparametric periodogram for two sampling schemes: the random sampling pattern and the sampling pattern generated by the integral pulse frequency modulation (IPFM), the latter widely accepted as the heart rate timing model. A general formula that relates the variance of timing deviations from a regular scheme with the aliased component suppression is proposed. The derived relation is illustrated by Lomb-Scargle periodograms applied on simulated data. Presented experimental data consisting of the respiration signal derived from the electrocardiogram and the heart rate signal also support possibility to detect frequencies above the Nyquist limit in the condition known as the cardiac aliasing.

1. Introduction

Aliasing is a well-known phenomenon in signal theory that occurs if a continuous-time signal is sampled at a frequency 𝑓𝑠 below twice the maximum frequency 𝑓𝑚 of the sampled signal, that is, 𝑓𝑠<2𝑓𝑚. The sampling frequency halved, 𝑓𝑁=𝑓𝑠/2, is referred to as the Nyquist frequency and can be interpreted as the presumed maximum frequency of the sampled signal provided that a sampling was performed at a sufficiently high rate. If the frequency of the sampled signal exceeds the Nyquist limit, this frequency appears at a wrong position, in the range 0𝑓𝑠/2, unless an antialiasing filter was applied. These rules are, however, not fully applicable in nonuniformly sampled data.

Nonuniform sampling can be exploited intentionally in a system design, or it can reflect specific restrictive conditions in measurement and exploration of data in some scientific fields, such as astronomy, seismology, genetics, biophysical measurements, and cardiac physiology. Data cannot be acquired at prescribed time instants, they are missing or corrupted. In other situations, essentially irregular sampling times are dictated by nature, such as it is in analysis of the heart rate variability. Nonuniform sampling makes analysis of a signal more complicated, because great majority of data analysis and digital signal processing procedures have been developed for uniformly sampled signals. On the other hand, nonuniformly sampled signals are more informative and components beyond the Nyquist limit can be detected.

Notion of the Nyquist frequency needs some clarification when considering nonuniformly sampled signals. There is no well-accepted definition if the sampling frequency is not constant [1]. Some authors define the Nyquist frequency in the nonuniform case as 1/(2𝑡min), where 𝑡min is the smallest time interval between data, while others infer the Nyquist frequency from average sampling times 𝑡mean and defines it as 1/(2𝑡mean). Amore general way of defining the Nyquist frequency for nonuniform samples uses a concept of the spectral window, that is, magnitude squared Fourier transform (FT) of the sampling pattern. The spectral window reaches its maximum at the zero frequency and a secondary maximum close to the main maximum at the frequency that is stated to be the Nyquist frequency doubled. In this paper, we will consider sampling times that deviates from regular instants in a homogenous way, that is, without marked clusters or long gaps. In such a situation, both average sampling time and spectral window concepts give practically identical results.

The spectral analysis of nonuniformly sampled data has relatively old history due to a theoretical and practical significance of this problem. Various approaches have been proposed to deals with nonuniform samples. A comprehensive review of the state of the art can be found in [1]. We will restrict our consideration to the methods commonly used in the research of the heart rate data which are used in this paper as demonstrative data to illustrate analyzed phenomena. The simplest method to deal with nonuniformly sampled data is to ignore sampling irregularity, relocate the samples into a regular grid, and treat them as a uniformly sampled signal. Such approach results in a reduction of spectral amplitudes and sampling irregularity is translated into a spectral noise. Components above the Nyquist frequency cannot be distinguished from aliased frequencies. Interpolation-based methods attenuate higher frequencies making them unsuitable for analysis of components above the Nyquist frequency. A widely used method that does not assume a signal resampling is known as the Lomb-Scargle periodogram, which is basically least-squares fitting of measured data by sinusoids [26].

The Lomb-Scargle periodogram (LS) [2] is a method especially designed for nonuniformly acquired samples. This method was originally proposed for the analysis of astronomic time series and also found successful application in the analysis of cardiovascular data, such as those in the analyses of the heart rate variability [7, 8], circadian rhythms [9], and gene sequence analysis. Results of LP are in many situations almost identical to the classical (Schuster) periodogram generalized to nonuniform samples [10] and thus they can be inferred from the FT of nonuniformly sampled data. Such an approach is used in [6] where several properties of the nonuniformly sampled FT are presented, although they are referred imprecisely as properties of the Lomb power spectral density.

Unlike our paper, the analyses [6] have been performed without a special attention to quantify an effect of the aliasing, that is, ability to uncover components above the Nyquist limit. This paper studies the utility of nonuniform sampling for both reducing aliased components and revealing components above the Nyquist rate. We will provide theoretical analysis of a model situation, when a sinusoidal signal is sampled according to a prespecified nonuniform sampling pattern. A general formula that relates the variance of timing deviations from a regular scheme with the aliased component suppression is proposed. A theoretically proven observability of components beyond the Nyquist limit is illustrated by simulation examples and cardio-respiratory data obtained from an experiment especially arranged for this purpose.

Outline of the paper is as follows. Section 2 summarizes definitions of periodograms used in the spectral analysis of nonuniform data. Subsequent sections present theoretical spectral analysis of the sinusoid sampled with two sampling schemes. In Section 3 of the paper, we will derive the formulas for power spectral densities (PSDs) of nonuniformly sampled sinusoidal signals provided that the sampling times obey a random model. Section 4 provides a treatment of sinusoid spectrum when the sampling pattern was generated by the integral pulse frequency modulation (IPFM). A quantitative formula describing an aliased component reduction is derived in Section 5. Demonstrative analyses of simulated data are presented in Section 6. Spectral analyses of experimental data—electrocardiogram- (ECG-) derived respiration and the heart rate variability (HRV) signals—illustrate performance of the LS method in the case of so-called cardiac aliasing, when components above the Nyquist frequency were elicited by the respiratory sinus arrhythmia.

2. Periodograms for Nonuniformly Sampled Signals

The conventional periodogram based on the magnitude squared Fourier transform can be generalized also for nonuniformly sampled signals in the form [4]: 1𝑃(𝜔)=𝑁|||||𝑁1𝑛=0𝑥𝑡𝑛𝑒𝑗𝜔𝑡𝑛|||||2,(1) where 𝑥(𝑡𝑛) are data samples taken at the time instants 𝑡𝑛. Notice that roots of a Fourier-type analysis with complex exponentials reside in a least-squares fitting of data to the complex sinusoid 𝑒𝑗𝜔𝑡. Lomb [2] formulated the periodogram by means of the least-squares fitting of observed real-valued data to the model of the form: 𝜔𝑡𝐴cos𝑖𝜔𝑡𝜏+𝐵sin𝑖𝜏.(2) The redundant time delay variable 𝜏, computed as tan(2𝜔𝜏)=𝑁1𝑛=0sin2𝜔𝑡𝑛𝑁1𝑛=0cos2𝜔𝑡𝑛,(3) was introduced into the model (2) in order to orthogonalize sine and cosine terms, which arise in a formulation of the normal equations of the fitting problem. The LS PSD is then computed as [4] 𝑃LS1(𝜔)=2𝑁1𝑛=0𝑥𝑡𝑛𝜔𝑡cos𝑛𝜏2𝑁1𝑛=0cos2𝜔𝑡𝑛+𝜏𝑁1𝑛=0𝑥𝑡𝑛𝜔𝑡sin𝑛𝜏2𝑁1𝑛=0sin2𝜔𝑡𝑛.𝜏(4) Provided that the summed sine and cosine squared terms in (4) are both equal and approach to 𝑁/2, and 𝑥 is shifted by 𝜏, the LS (4) approaches to the classical periodogram (1). Generally, (1) and (4) are not equivalent. The LS is preferred due to statistical properties useful for significance testing of spectral peaks [4, 11]. On the other hand, (1) is mathematically more tractable.

3. Random Sampling Pattern

The analyzed signal will take a form of the complex sinusoid with the amplitude 𝐶, frequency 𝜔0, sampled at the time instants 𝑡𝑛: 𝑥𝑡𝑛=𝐶𝑒𝑗𝜔0𝑡𝑛.(5) The sampling instants are supposed to be the random numbers obeying the following model: 𝑡𝑛=𝑛𝑇0+𝜉𝑛,(6) where 𝑇0 is the mean sampling period, 𝑛=0,,𝑁1, and 𝜉𝑛 is a sequence of independent identically distributed random variables. Average behavior of (1), approximately valid also for the LS spectrum, can be deduced by applying the expectation on the Fourier transform of (5): 𝑋[](𝜔)=𝐸𝑋(𝜔)=𝐸𝑁1𝑛=0𝑥𝑡𝑛𝑒𝑗𝜔𝑡𝑛.(7) After substitution (5), (6) into (7) and rearrangement, we obtain 𝐸[]𝑋(𝜔)=𝐶𝑁1𝑛=0𝑒𝑗(𝜔𝜔0)𝑛𝑇0𝐸𝑒𝑗(𝜔𝜔0)𝜉𝑛.(8) In (8), we can recognize FT of a uniformly sampled sinusoid expressible by the periodic Dirichlet kernel: 𝑁𝐷𝑁𝜔𝜔0𝑇0,(9) where the kernel is 𝐷𝑁(1𝑥)=𝑁sin(𝑁𝑥/2).sin(𝑥/2)(10) The mean FT (8) can be rewritten in terms of the characteristic function shifted by the frequency 𝜔0: 𝑋(𝜔)=𝐶𝑁𝐷𝑁𝜔𝜔0𝑇0𝑇𝜔𝜔0,(11) where the characteristic function 𝑇 of the timing irregularity 𝜉, characterized by the probability density function 𝑝𝜉, is defined as 𝑒𝑇(𝜔)=𝐸𝑗𝜔𝜉=𝑝𝜉(𝜉)𝑒𝑗𝜔𝜉𝑑𝜉.(12) Therefore, the characteristic function attenuates periodic peaks—aliases or spectral images in (11), while it retains the true peak at the frequency 𝜔0.

The resultant expectation of the FT (11) does not account thoroughly for a random effect of the sampling irregularity. Therefore, the expectation of a magnitude squared FT, that is, PSD by its definition, is evaluated to have a more complete though on the spectrum of an irregularly sampled signal: 𝑃[]=1(𝜔)=𝐸𝑃(𝜔)𝑁𝐸𝑁1𝑛=0𝑥𝑡𝑛𝑒𝑗𝜔𝑡𝑛𝑁1𝑘=0𝑥𝑡𝑘𝑒𝑗𝜔𝑡𝑘=1𝑁𝐸𝑁1𝑛=0𝑁1𝑘=0𝑥𝑡𝑛𝑥𝑡𝑘𝑒𝑗𝜔(𝑡𝑛𝑡𝑘).(13) The summation term in (13) after substitution (5) is 𝐶2𝑁𝐸𝑁1𝑛=0𝑁1𝑘=0𝑒𝑗𝜔0(𝑛𝑘)𝑇0𝑒𝑗𝜔(𝑛𝑘)𝑇0𝑒𝑗(𝜔𝜔0)𝜉𝑛𝑒𝑗(𝜔𝜔0)𝜉𝑘.(14) For 𝑛𝑘, we have due to the presumed independence: 𝐸𝑒𝑗𝜔𝜉𝑛𝑒𝑗𝜔𝜉𝑘=𝑇(𝜔)𝑇(𝜔)=𝑇2(𝜔),(15) where we have assumed a symmetrical probability function in (12) and thus a real-valued characteristic function.

In the specific case, 𝑛=𝑘, expectation (15) is 1 and it can be formally expanded as 𝐸𝑒𝑗0=𝑇2𝜔𝜔0+1𝑇2𝜔𝜔0.(16) After substitution 𝑚=𝑛𝑘 in (14) and taking into account (15), (16), the summation (14) simplifies to 𝑃(𝜔)=𝐶21𝑇2𝜔𝜔0+𝐶2𝑇2𝜔𝜔0𝑁1𝑚=𝑁+1|1𝑚|𝑁𝑒𝑗𝜔0𝑚𝑇0𝑒𝑗𝜔𝑚𝑇0.(17) The summation term in (17) can be recognized as the FT of a uniformly sampled sinusoid multiplied by the triangular window. In order to facilitate the formula readability, this FT will be denoted as 𝑆𝑁𝜔𝜔0=𝑁𝐷2𝑁𝜔𝜔0𝑇0.(18) The final formula then becomes 𝑃(𝜔)=𝑇2𝜔𝜔0𝐶2𝑆𝑁𝜔𝜔0+𝐶21𝑇2𝜔𝜔0.(19) Unlike (11), the PSD expression (19), besides the sinusoid-related peaks, includes a smooth spectrum portion related to the randomness of the timing instants. Since 𝑇(0)=1, this smooth portion exhibits a valley near the true frequency of the complex sinusoid. Composition of PSD from the terms present in (19) is clarified in Figure 1.

643563.fig.001
Figure 1: Components of the nonuniformly sampled PSDthe random sampling pattern.

4. IPFM Model-Generated Sampling Pattern

The integral pulse frequency modulation (IPFM) model is widely accepted as the sinoatrial node timing model in modeling the heart rate variability. It is an integrate-and-fire model, also used to describe a neuronal activity. The IPFM model assumes that a modulating signal is integrated and a beat trigger pulse is generated when the integrated signal reaches a threshold (Figure 2). The beat occurrence times (representing sampling times) 𝑡𝑘 generated by the IPFM model can be defined as [12] 𝑘=𝑡𝑘𝑡01+𝑚(𝑡)𝑇0𝑑𝑡,(20) where 𝑘 is an integer, 𝑡0=0, 𝑇0 is the mean heart period, and 𝑚(𝑡) is a dimensionless modulating signal interpreted as the fractional change of the instantaneous heart rate relative to 1/𝑇0.

643563.fig.002
Figure 2: IPFM model for generating the beat occurrence times.

The sampling pattern can be written as a sequence of the Dirac pulses: spc(𝑡)=𝑘=𝛿𝑡𝑡𝑘.(21) The Fourier transform of (21) is referred to as the spectrum of counts (SPC) [12]: SPC(𝜔)=𝑘=𝑒𝑗𝜔𝑡𝑘.(22) The spectrum of a sampled signal 𝑋𝑠(𝜔) is related to the true spectrum 𝑋(𝜔) by the frequency-domain convolution: 𝑋𝑠1(𝜔)=2𝜋𝑋(𝜔)SPC(𝜔).(23) The spectrum of counts has been extensively analyzed in several works [1214]. We will summarize here facts required to assess the effect of the spectral aliasing. The SPC defined by (22) can be rewritten in terms of the modulating signal as 1SPC(𝜔)=𝑇01(2𝜋𝛿(𝜔)+𝑀(𝜔))𝛿(𝜔)+𝜋𝑛0=1𝐹𝑇cos2𝜋𝑛0𝑇0𝑡+𝑡.𝑚(𝜏)𝑑𝜏(24) Interestingly, the spectrum of the modulating signal 𝑀(𝜔) is directly included in (24). Beside the Dirac pulse located at the frequency origin, infinite numbers of frequency modulated (FM) terms convolved with 𝑀(𝜔) are present. The carriers of these FM modulated terms producing peaks in (24) are located at multiples of the mean sampling frequency, which explain an origin of the aliased components in the spectrum (23). For the sinusoidal modulating signal: 𝑚(𝑡)=𝑀1𝜔cos𝑚𝑡,(25) where 𝜔𝑚 is the modulating frequency, the spectrum (24) can be expanded using the Bessel functions of the first kind [15]. Spectral amplitudes of the carriers follow the 0th-order Bessel function 𝐽0, the argument of which is given by the gradually increased index of the frequency modulation: 𝑛𝛽=0𝜔𝑠𝑀1𝜔𝑚,(26) where 𝜔𝑠 is the mean sampling frequency corresponding to the first carrier (𝑛0=1).

5. Aliasing Reduction Index

For quantification of the aliased components amplitudes, we will restrict our considerations to the situation when the input frequency is in the range 𝑓𝑠/2 to 𝑓𝑠. Unlike the derivation presented in Section 3, a real-valued sinusoid, comprising of positive and negative frequencies, will be considered in this treatment. It is well known that the spectrum of such a uniformly sampled signal is mirrored around the Nyquist frequency. If the frequency of the sampled real-valued signal exceeds the Nyquist limit, this frequency is folded over the Nyquist frequency and appears at a wrong position, in the range 0𝑓𝑠/2. Specifically, if an input signal has frequency 𝑓𝑖,𝑓𝑁<𝑓𝑖<𝑓𝑠, the aliased component appears at the frequency 𝑓𝑎=𝑓𝑠𝑓𝑖. If the spectrum is plotted in a range above 𝑓𝑠/2, both the spectral peaks, 𝑓𝑖 and 𝑓𝑎, have equal amplitudes and there is no possibility to distinguish them. For nonuniformly sampled signal, this statement does not hold, and differences can be observed. A component with higher amplitude is considered to be the true component and other component at mirrored frequency as its alias. We introduce a quantity to describe the observed difference between the true and aliased components and name this quantity as the aliasing reduction index (ARI): 𝐴ARI=𝑖𝐴𝑎𝐴𝑖,(27) where 𝐴𝑖 is the spectral amplitude at the true (input) frequency 𝑓𝑖, 𝐴𝑎 is the spectral amplitude at the corresponding aliased frequency 𝑓𝑎. Next, we will be dealing with an evaluation of the proposed ARI for both considered sampling models—the random sampling pattern and the IPFM-generated sampling pattern.

The reduction of an aliased component in the situation with the random sampling according to (6) is deduced directly from (11). For a true component, the argument in (11) is 𝜔=+𝜔𝑖=𝜔0, while the component aliased into 0𝑓𝑠/2 is originated from the negative input frequency 𝜔𝑖, shifted by value 𝜔𝑠, that is, 𝜔𝑖=𝜔0 and 𝜔=𝜔𝑖+𝜔𝑠. The ARI can be then written in terms of the characteristic function as ARIRND=𝑇𝜔(0)𝑇𝑠.𝑇(0)(28)

The result (28) depends on a particular probability distribution function which describes the random sampling jitter. For small-to-moderate values of ARI, a more general version of (28) can be derived using power series expansion of the characteristic function. Taking into account the presumed symmetry of the probability density function and neglecting the high-order terms in the expansion, the approximated characteristic function is expressed by the variance of the timing deviations: 1𝑇(𝜔)1+𝑇2!(0)𝜔21=12𝜎2𝜔2.(29) Using 𝜔𝑠=2𝜋/𝑇0 in (29), 𝑇(0)=1 in (28), we obtain an approximate formula for the ARI in term of the relative timing standard deviation 𝜎/𝑇0: ARIRND2𝜋2𝜎𝑇02,(30) which is independent of a particular probability distribution type.

The proposed index can be evaluated also for the IPFM-generated sampling pattern. The spectrum (23) of a real-valued sinusoid, nonuniformly sampled in IPFM-generated time instants, is composed from SPC shifted in the frequency direction by +𝜔𝑖 and by 𝜔𝑖. A dominant aliased term is associated with the 1st FM carrier, that appears, after 𝜔𝑖 shifting by 𝜔𝑠, at the frequency 𝜔𝑎=𝜔𝑠𝜔𝑖. Its amplitude is determined by the 0th-order Bessel function of the first kind for 𝛽 given by (26), 𝑛0=1. The resulting formula for the ARI estimation, analogical to (28), becomes ARIIPFM=𝐽0(0)𝐽0(𝛽)𝐽0.(0)(31) Since 𝐽0(0)=1, and using expansion of 𝐽0: 𝐽0(𝑥)=𝑘=0(1)𝑘(𝑘!)214𝑥2𝑘114𝑥2,(32) ARI can be approximated by ARIIPFM14𝛽2.(33) In order to facilitate a comparison of the obtained formula (33) with the ARI derived for the random sampling model, we will rewrite (33) by using the standard deviation of the timing instants. Integration of the IPFM model formula (20) gives𝑡𝑘𝑘𝑇0𝑀=1𝜔𝑚𝜔sin𝑚𝑡𝑘.(34) The deviations of the sampling times from a regular scheme is thus sinusoidal, with amplitude 𝑀1/𝜔𝑚, and the standard deviation estimated as the RMS value of a continuous-time sinusoid is 1𝜎=2𝑀1𝜔𝑚.(35) Combining (33) with (35) and (26) results in the formula: ARIIPFM2𝜋2𝜎𝑇02.(36) Surprisingly, (36) has exactly the same form as (30), despite the essentially different sampling pattern.

6. Simulation Example

In order to illustrate a mechanism of the aliasing in nonuniformly sampled signal and show an extent to which an aliased component reduction can be achieved, we have performed a spectral analysis of sequences synthesized according to both described sampling models. Since this work is inspired by analysis of the cardiac aliasing effect, we have arranged the signal parameters—duration, frequency—as to fit the experimental data presented in subsequent section. Each of the simulated signals consists of 120 samples of real-valued, unit-amplitude sinusoid. The frequencies are expressed as the fractions of the mean sampling rate 𝑓𝑠. The mean sampling rate is considered as the reciprocal value of the mean sampling interval. An input frequency exceeding the Nyquist limit was set to 0.65𝑓𝑠 which is aliased to the frequency 0.35𝑓𝑠. Conversely, an input frequency 0.35𝑓𝑠 below the Nyquist limit 0.5𝑓𝑠 will be imaged to 0.65𝑓𝑠 and an informationless peak in the spectrum can be observed when a periodogram is plotted above the Nyquist limit. Both spectral peaks appear as identical when the sampling is uniform. The nonuniform sampling, however, introduces differences between the true component and the unwanted—aliased or imaged components. This effect can be observed in the spectra plotted in the range 0𝑓𝑠 instead of a conventionally used range 0𝑓𝑠/2.

Figures 3 and 4 show LS periodograms for both scenarios (below and above Nyquist limit) for the sampling times randomly deviated from regular positions. Spectral amplitudes are plotted as the square roots of the LS periodograms and scaled to directly reflect amplitudes of the generated sinusoids. The sampling irregularity follows the Gaussian distribution with the variance set according to a required irregularity level. Three levels were used in simulation: low, medium, high. The low irregularity is represented by sampling jitter with the standard deviation 5% of the mean sampling time (𝜎=0.05𝑇0). Such a degree of the sampling nonuniformity caused only a minor spectral amplitude difference (about 5%). The medium level is represented by two-fold reduction of alias power (3 dB), that is, ARI = 0.707. The formula derived in Section 5 allowed to estimate 𝜎 required in simulation as 0.122𝑇0. A maximum deviation of a sampling instant is theoretically unlimited; therefore, we incorporate a practical limit as the irregularity level which ensures a monotonic increase of the sampling times. This practical maximum is used as the high level of the irregularity. Since Gaussian variables are not confident to finite interval, determination of the high level deviations must be made in probabilistic sense. The choice 𝜎=0.3𝑇0 ensures that randomly chosen proximate sampling times 𝑡𝑘1,𝑡𝑘, fulfill the monotonicity condition 𝑡𝑘1<𝑡𝑘 with probability 0.99. Presented figures show a gradual destruction of the aliased component (Figure 4) and spectral images (Figure 3) as the sampling irregularity level increases. In the high-level irregularity, the aliasing is effectively destroyed—the aliased peak is buried in a noise caused by sampling irregularity. Notice that a noise-like spectrum increase is evident as the sampling irregularity increases. A valley near the true frequency in noise PSD, which may be expected by examination of (19), is not evident, because it is overlapped by the smooth noise spectrum coming from the negative frequency.

fig3
Figure 3: Spectral amplitudes of the LS spectra—nonuniformly sampled signal with the random sampling pattern. The input frequency 0.35𝑓𝑠 lies below the Nyquist limit. Different levels of sampling irregularity: (a) low, (b) medium, (c) high.
fig4
Figure 4: Spectral amplitudes of the LS spectra—nonuniformly sampled signal with the random sampling pattern. The input frequency 0.65𝑓𝑠 exceeds the Nyquist limit. Different levels of sampling irregularity: (a) low, (b) medium, (c) high.

A similar set of simulations (Figures 5 and 6) was performed for the sinusoid sampled by the pattern generated by the IPFM model driven by a sinusoidal modulating signal. The sampling time instants were obtained by means of numerical solution of the IPFM equation (34) using MATLAB function fzero. The frequency of the sampled sinusoid was set equal to the frequency of the modulating signal, a condition occurred in heart rate modulation by the respiration signal. The low and medium irregularity levels were established in the same way as in the random sampling pattern. The maximum irregularity level is considered when the normalized modulation signal amplitude 𝑀1 reaches the value1. Due to this restriction, the achieved attenuation cannot be as high as in the random sampling pattern scenarios. Notice that numbers of spurious peaks can be found in the spectrum which are explained by a complex structure of the SPC (24).

fig5
Figure 5: Spectral amplitudes of periodograms—nonuniformly sampled signal with IPFM generated sampling pattern. Input frequency 0.35𝑓𝑠 lies below the Nyquist limit. Different levels of sampling irregularity: (a) low, (b) medium, (c) high.
fig6
Figure 6: Spectral amplitudes of periodograms—nonuniformly sampled signal with IPFM generated sampling pattern. Input frequency 0.65𝑓𝑠 exceeds the Nyquist limit. Different levels of sampling irregularity: (a) low, (b) medium, (c) high.

7. Cardiorespiratory Data Analysis

Beat-to-beat variations of consecutive heart periods or the instantaneous heart rate, known as the heart rate variability (HRV), become intensively studied over the last decades. Analysis of the HRV provides a noninvasive insight into the cardiovascular neural control and a spectral analysis of HRV signals is a widely used method for assessment of the autonomous nervous system. Generally, the power spectra of HRV can be divided into the high-frequency (HF, 0.15–0.4 Hz), low-frequency (LF, 0.04–0.15 Hz), and very-low-frequency range (VLF, 0.003–0.04 Hz). The LF power is modulated by both sympathetic and parasympathetic controls, while the HF power mainly reflects the parasympathetic influence linked to the respiration. Therefore, a concurrent measurement of the respiration signal is helpful for HRV interpretation. To obtain the respiration activity signal, an additional specialized instrument is not necessary because the respiration can be estimated from modulation of an electrocardiogram amplitude. Such a procedure is referred to as ECG-derived respiration [1618]. Notice that a raw ECG signal itself includes spectral components related to heart rate oscillations and respiration [19], but they are mixed in the spectrum in a complex way [20], and thus, appropriate signals must be derived from measured ECG. Both the HRV as well as the respiration signal derivation incorporate heart beat detection, and the beat occurrence times define sampling instants, irregularly spaced by nature. Although a heart rate (HR) is essentially a discrete-time signal, computed from beat-to-beat occurrence times, it can be considered as a hypothetically continuous signal representing the autonomous nervous system activity which continually modulates the heart rate. Unlike an artificially designed system, a physiological system does not include antialiasing filter. For example, 0.4 Hz frequency considered as the highest limit of the standardized HF range requires sampling at a mean heart rate 0.8 Hz, and thus, the mean heart rate of at least 48 bpm is required for reliable spectral analysis. Therefore, a bradycardia or an extended HF range caused by the elevated respiration rate can lead to the aliasing. For this kind of “physiological aliasing,” the term cardiac aliasing has been naturalized. High frequency components modulating the heart rate have been actually observed in neonates, transplanted heart patients, and animals [21, 22]. Improper signal processing or ignoring rarely occurred conditions can then result in incorrect interpretation of a spectral analysis.

To test an ability of different methods to detect frequencies beyond the Nyquist rate in a real-world measured data, we have performed an experiment exploiting respiratory sinus arrhythmia. The measured subject was asked to breathe at an elevated rate in synchrony with a moving pattern displayed on the computer screen at the frequency of 0.8 Hz. An electrocardiogram (ECG) was recorded and two signals were extracted: the respiration and the instantaneous heart rate (HR). The respiration signal was obtained as the amplitude in bandpass-filtered ECG curve taken at maxima of QRS complexes [23]. The HR was computed from beat-to-beat time distances of the QRS maxima. Both signals are thus sampled at a variable rate dictated by the heart rate. The mean heart rate of 74 bpm (1.23 Hz) was sufficiently low to observe an aliasing.

In Figure 7(a), the respiration signal is analyzed as a sequence of values equidistantly spaced at multiples of the mean heart period. A meaningful spectral range is thus up to 0.615 Hz provided that the mean heart rate was 1.23 Hz. The respiration frequency 0.8 Hz is aliased near 0.4 Hz and spreads over a wider spectral range due to wandering of the mean heart rate. The respiration sequence was then interpolated to 4 Hz sampling frequency by means of the cubic spline interpolation. A small peak at the true respiratory frequency can be detected but the aliased component near 0.4 Hz largely exceeds the true peak in its amplitude—Figure 7(b). The Lomb-Scargle periodogram manifests a clearly visible peak at 0.8 Hz, Figure 7(c), albeit the components at aliased frequency region do not disappear completely.

fig7
Figure 7: Spectral analysis of the ECG-derived respiration signal: (a)equidistantly repositioned QRS maxima samples with frequency axis scaled according to the mean heart period, (b) spline interpolated unevenly sampled QRS maxima, (c) the Lomb-Scargle periodogram of unevenly sampled QRS maxima.

The results of application of different methods for HRV spectrum computation are summarized in Figure 8. Repositioning of unevenly spaced samples, similarly as it was observed in Figure 7(a), yields a symmetrical spectrum that is unable to detect components above the mean heart rate halved. The cubic spline interpolation, although frequently used in HRV analysis, noticeably attenuates high-frequency components above the mean Nyquist frequency. The simplest interpolation type, nearest neighbor method, gives surprisingly better outcome, comparable to the Berger resampling method [24]. As in the case of the respiration signal analysis, only the last two methods presented, the Lomb-Scargle periodogram and the SPC, are capable to show a higher amplitude at the true position 0.8 Hz than near the aliased frequency (about 0.4 Hz). The apparent increase of the amplitudes at higher frequencies in Figure 8(f) can be explained by low-pass filtering effect of the integrator in the IPFM model that does not affect the spectrum of counts. In the low frequency range, all methods give similar results.

fig8
Figure 8: HRV spectra obtained by means of different methods: (a)equidistantly repositioned HR samples with frequency axis scaled according to the mean heart period; unevenly sampled heart rate sequence resampled to evenly spaced signal- (b) cubic-spline interpolation, (c) nearest neighbour interpolation, (d) Berger method; spectra computed without resampling-, (e) Lomb-Scargle periodogram, (f) spectrum of counts.

8. Conclusion

The fact that a spectrum of nonuniformly sampled signal conveys useful information above the Nyquist limit has been pointed out in several works [25, 26]. In this work, we have presented a theoretical treatment which quantitatively explains this phenomenon. The analyses of simulated and experimental data illustrate a possibility to identify components beyond the Nyquist limit in real-world problems.

Although the spectrum of a nonuinformly sampled signal is not alias free, unlike a uniform sampling situation, the aliased spectral amplitudes are attenuated when compared to the true components amplitudes. Analogically, components below the Nyquist limit can be observed also as images above this limit. The analytic treatment of these unwanted components constitutes a key contribution of this paper. We have studied the aliasing effect in two sampling schemes: random deviations from regular times and sampling times generated by the IPFM model. The former was chosen as a basic model which allowed to explain a mechanism of aliasing suppression. The later is inspired by a model used to describe generation of events in biological systems, such as the sinus node or the neuronal firing.

The two essentially distinctive sampling models required different way of analysis. The random model analysis relies on terms common in probability/statistics. We have shown that the characteristic function of timing jitter plays a key role in composition of the resulting spectrum. The characteristic function modulates amplitudes of aliased or imaged spectral lines. Besides the line portion of the spectrum, a smooth noise portion shaped by a term of characteristic function is present due to a random nature of the sampling process, despite the sampled signal is deterministic. The IPFM sampling scheme uses concepts originated in communication systems theory and presented analysis is purely deterministic. Consequently, the smooth spectrum is not present; on the other hand, spurious peaks frequently occur due to a complex nature of the spectrum of the modulated signal. In the both situations, the relative reduction of aliased components was evaluated by a quantity named in this paper as the aliasing reduction index. Interestingly, the derived formula approximating this index for small-to-moderate sampling irregularity is identical for the both sampling models. The formula which we have found seems to be quite universal and relates a decrease of unwanted aliased components with the standard deviation of timing irregularity measured as the fraction of the mean sampling period. Validity of the formula in more general sampling schemes, such as those with dependent sampling irregularities and additive random sampling, is intended to be studied in a future work.

Mitigation of the aliasing effect by means of nonuniform sampling not only attracts scientific interest but also has practical implications. We have demonstrated the possibility to detect frequencies above the Nyquist limit in the ECG-derived respiration signal and the heart rate signal. The condition of aliasing, so-called cardiac aliasing, was elicited by increasing the respiration rate, and signals derived from recorded electrocardiogram were analyzed by different methods. Conventionally used interpolation/resampling methods were found to be unsuitable in the condition of cardiac aliasing. The Lomb-Scarglre method, the classical periodogram used for nonuniform samples and the spectrum of counts are all capable to uncover components above the Nyquist frequency. Since real-world measurements are usually contaminated by noise, a signal to be detected needs not to be a pure sinusoid, or it can be random and even nonstationary, development of a universal method able to resolve between true and aliased spectral components, with properly assigned significance measure, is another challenging task intended for future work.

References

  1. P. Babu and P. Stoica, “Spectral analysis of nonuniformly sampled data—a review,” Digital Signal Processing, vol. 20, no. 2, pp. 359–378, 2010. View at Publisher · View at Google Scholar · View at Scopus
  2. N. R. Lomb, “Least-squares frequency analysis of unequally spaced data,” Astrophysics and Space Science, vol. 39, no. 2, pp. 447–462, 1976. View at Publisher · View at Google Scholar · View at Scopus
  3. P. Stoica, J. Li, and H. He, “Spectral analysis of nonuniformly sampled data: a new approach versus the periodogram,” IEEE Transactions on Signal Processing, vol. 57, no. 3, pp. 843–858, 2009. View at Publisher · View at Google Scholar · View at Scopus
  4. J. D. Scargle, “Studies in astronomical time series analysis. II—statistical aspects of spectral analysis of unevenly spaced data,” Astrophysical Journal, vol. 263, pp. 835–853, 1982. View at Publisher · View at Google Scholar
  5. P. Vaníček, “Approximate spectral analysis by least-squares fit-Successive spectral analysis,” Astrophysics and Space Science, vol. 4, no. 4, pp. 387–391, 1969. View at Publisher · View at Google Scholar · View at Scopus
  6. P. Laguna, G. B. Moody, and R. G. Mark, “Power spectral density of unevenly sampled data by least-square analysis: performance and application to heart rate signals,” IEEE Transactions on Biomedical Engineering, vol. 45, no. 6, pp. 698–715, 1998. View at Publisher · View at Google Scholar · View at Scopus
  7. “Heart rate variability: standards of measurement, physiological interpretation and clinical use. Task Force of the European Society of Cardiology and the North American Society of Pacing and Electrophysiology,” Circulation, vol. 93, no. 5, pp. 1043–1065, 1996. View at Scopus
  8. G. G. Berntson, J. T. Bigger, D. L. Eckberg et al., “Heart rate variability: origins, methods, and interpretive caveats,” Psychophysiology, vol. 34, no. 6, pp. 623–648, 1997. View at Google Scholar · View at Scopus
  9. M. Teplan, L. Molčan, and M. Zeman, “Spectral analysis of cardiovascular parameters of rats under irregular light-dark regime,” in Proceedings of the 8th International Conference on Measurement, pp. 343–346, Smolenice, Slovak Republic, 2011.
  10. V. V. Vityazev, “Time series analysis of unequally spaced data: intercomparison between estimators of the power spectrum,” in Proceedings of the Astronomical Data Analysis Software and Systems VI ASP Conference Series, vol. 125, pp. 166–169, 1997.
  11. A. Schwarzenberg-Czerny, “The distribution of empirical periodograms: Lomb-Scargle and PDM spectra,” Monthly Notices of the Royal Astronomical Society, vol. 301, no. 3, pp. 831–840, 1998. View at Google Scholar · View at Scopus
  12. J. Mateo and P. Laguna, “Improved heart rate variability signal analysis from the beat occurrence times according to the IPFM model,” IEEE Transactions on Biomedical Engineering, vol. 47, no. 8, pp. 985–996, 2000. View at Publisher · View at Google Scholar · View at Scopus
  13. M. Brennan, M. Palaniswami, and P. Kamen, “Distortion properties of the interval spectrum of IPFM generated heartbeats for heart rate variability analysis,” IEEE Transactions on Biomedical Engineering, vol. 48, no. 11, pp. 1251–1264, 2001. View at Publisher · View at Google Scholar · View at Scopus
  14. J. Púčik, O. Ondráček, E. Cocherová, and A. Sultan, “Spectrum of counts computation for HRV analysis,” in Proceedings of 19th International Conference Radioelektronika (RADIOELEKTRONIKA '09), pp. 255–258, April 2009. View at Publisher · View at Google Scholar · View at Scopus
  15. E. J. Bayly, “Spectral analysis of pulse frequency modulation in the nervous systems,” IEEE Transactions on Biomedical Engineering, vol. 15, no. 4, pp. 257–265, 1968. View at Google Scholar · View at Scopus
  16. G. Moody, R. Mark, A. Zoccola, and S. Mantero, “Derivation of respiratory signals from multi-lead ECGs,” IEEE Computer Society, vol. 12, pp. 113–116, 1985. View at Google Scholar
  17. G. D. Clifford, F. Azuaje, P. E. McSharry et al., Advanced Methods and Tools for ECG Data Analysis, Artech House, 2006.
  18. D. Widjaja, J. Taelman, S. Vandeput et al., “ECG-derived respiration: comparison and new measures for respiratory variability,” in Proceedings of the Computing in Cardiology (CinC '10), pp. 149–152, September 2010. View at Scopus
  19. A. Gersten, O. Gersten, A. Ronen, and Y. Cassuto, “The RR interval spectrum, the ECG signal and aliasing,” . Eprint, http://arxiv.org/abs/physics/9911017v1.
  20. J. Šurda, S. Lovás, J. Púčik, and M. Jus, “Spectral properties of ECG signal,” in Proceedings of the of the 17th International Conference Radioelektronika, pp. 537–541, Brno, Czech Republic, 2007.
  21. E. Toledo, I. Pinhas, D. Aravot, and S. Akselrod, “Very high frequency oscillations in the heart rate and blood pressure of heart transplant patients,” Medical and Biological Engineering and Computing, vol. 41, no. 4, pp. 432–438, 2003. View at Google Scholar · View at Scopus
  22. H. A. Campbell, J. Z. Klepacki, and S. Egginton, “A new method in applying power spectral statistics to examine cardio-respiratory interactions in fish,” Journal of Theoretical Biology, vol. 241, no. 2, pp. 410–419, 2006. View at Publisher · View at Google Scholar · View at Scopus
  23. J. Púčik, M. Uhrík, A. Sultan, and A. Šurda, “Experimental setup for cardio-respiratory interaction study,” in Proceedings of the 8th Czech-Slovak Conference, Trends in Biomedical Engineering, pp. 126–129, Bratislava, Slovakia, 2009.
  24. R. Berger, S. Akselrod, D. Gordon, and R. Cohen, “An efficient algorithm for spectral analysis of heart rate variability,” IEEE Transactions on Bio-Medical Engineering, vol. 33, no. 9, pp. 900–904, 1986. View at Google Scholar · View at Scopus
  25. W. H. Press, B. P. Flannery, S. A. Teukolsky, W. T. Vetterling et al., Numerical Recipes in Fortran 77: The Art of Scientific Computing, vol. 1, Cambridge University Press, 1992.
  26. G. L. Bretthorst, “Nonuniform sampling: bandwidth and aliasing,” Concepts in Magnetic Resonance A, vol. 32, no. 6, pp. 417–435, 2008. View at Publisher · View at Google Scholar · View at Scopus