Observability of Spectral Components beyond Nyquist Limit in Nonuniformly Sampled Signals

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 [2–6].

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𝑥𝑡_𝑛𝑒^{−𝑗𝜔𝑡𝑛}|||||²,(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𝑛=0}sin2𝜔𝑡_𝑛∑^{𝑁−1𝑛=0}cos2𝜔𝑡_𝑛,(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] 𝑃_LS(𝜔)=12⎧⎪⎨⎪⎩∑^{𝑁−1𝑛=0}𝑥𝑡_𝑛cos𝜔𝑡_𝑛−𝜏²+∑∑^{𝑁−1𝑛=0}cos²𝜔𝑡_𝑛−𝜏_{𝑁−1𝑛=0}𝑥𝑡_𝑛sin𝜔𝑡_𝑛−𝜏²∑^{𝑁−1𝑛=0}sin²𝜔𝑡_𝑛−𝜏⎫⎪⎬⎪⎭.(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 𝜔₀, sampled at the time instants 𝑡_𝑛: 𝑥𝑡_𝑛=𝐶𝑒^{𝑗𝜔0𝑡𝑛}.(5) The sampling instants are supposed to be the random numbers obeying the following model: 𝑡_𝑛=𝑛𝑇₀+𝜉_𝑛,(6) where 𝑇₀ 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: 𝑁𝐷_𝑁𝜔−𝜔₀𝑇₀,(9) where the kernel is 𝐷_𝑁(𝑥)=1𝑁sisnin(𝑁(𝑥𝑥/2/2)).(10) The mean FT (8) can be rewritten in terms of the characteristic function shifted by the frequency 𝜔₀: 𝑋⁻(𝜔)=𝐶𝑁𝐷_𝑁𝜔−𝜔₀𝑇₀𝑇𝜔−𝜔₀,(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 𝜔₀.

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 𝐶²𝑁𝐸^𝑁−1_{𝑛=0𝑁−1𝑘=0}𝑒^{𝑗𝜔0(𝑛−𝑘)𝑇0}𝑒^{−𝑗𝜔(𝑛−𝑘)𝑇0}𝑒^{−𝑗(𝜔−𝜔0)𝜉𝑛}𝑒^{𝑗(𝜔−𝜔0)𝜉𝑘}.(14) For 𝑛≠𝑘, we have due to the presumed independence: 𝐸𝑒^{𝑗𝜔𝜉𝑛}𝑒^{−𝑗𝜔𝜉𝑘}=𝑇(−𝜔)𝑇(𝜔)=𝑇²(𝜔),(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=𝑇²𝜔−𝜔₀+1−𝑇²𝜔−𝜔₀.(16) After substitution 𝑚=𝑛−𝑘 in (14) and taking into account (15), (16), the summation (14) simplifies to 𝑃⁻(𝜔)=𝐶+²𝐶1₂−𝑇²𝜔−𝜔₀𝑇²𝜔−𝜔₀^𝑁−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 𝑆_𝑁𝜔−𝜔₀=𝑁𝐷^2𝑁𝜔−𝜔₀𝑇₀.(18) The final formula then becomes 𝑃⁻(𝜔)=𝑇²𝜔−𝜔₀𝐶²𝑆_𝑁𝜔−𝜔₀+𝐶²1−𝑇²𝜔−𝜔₀.(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.

Figure 1 

Components of the nonuniformly sampled PSD—the 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+𝑚(𝑡)𝑇₀𝑑𝑡,(20) where 𝑘 is an integer, 𝑡₀=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/𝑇₀.

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: 𝑋_𝑠(𝜔)=12𝜋𝑋(𝜔)∗SPC(𝜔).(23) The spectrum of counts has been extensively analyzed in several works [12–14]. 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 SPC(𝜔)=1𝑇₀∗⎧⎪⎨⎪⎩(𝛿2𝜋(𝜔𝛿)(𝜔+)1+𝑀(𝜔))𝜋^∞_𝑛0=1𝐹𝑇cos2𝜋𝑛⁰𝑇₀𝑡+^𝑡−∞𝑚(𝜏)𝑑𝜏⎫⎪⎬⎪⎭.(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: 𝑚(𝑡)=𝑀₁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 𝐽₀, the argument of which is given by the gradually increased index of the frequency modulation: 𝛽=𝑛⁰𝜔_𝑠𝑀₁𝜔_𝑚,(26) where 𝜔_𝑠 is the mean sampling frequency corresponding to the first carrier (𝑛₀=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 𝜔=+𝜔_𝑖=𝜔₀, while the component aliased into 0−𝑓_𝑠/2 is originated from the negative input frequency −𝜔_𝑖, shifted by value 𝜔_𝑠, that is, −𝜔_𝑖=𝜔₀ and 𝜔=−𝜔_𝑖+𝜔_𝑠. The ARI can be then written in terms of the characteristic function as ARI_RND=𝑇(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+12!𝑇^(0)𝜔²=1−12𝜎²𝜔².(29) Using 𝜔_𝑠=2𝜋/𝑇₀ in (29), 𝑇(0)=1 in (28), we obtain an approximate formula for the ARI in term of the relative timing standard deviation 𝜎/𝑇₀: ARI_RND≈2𝜋²𝑇𝜎₀²,(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), 𝑛₀=1. The resulting formula for the ARI estimation, analogical to (28), becomes ARI_IPFM=𝐽⁰(0)−𝐽₀(𝛽)𝐽₀(0).(31) Since 𝐽₀(0)=1, and using expansion of 𝐽₀: 𝐽₀(𝑥)=^∞_𝑘=0(−1)^𝑘(𝑘!)²41𝑥²^𝑘≈1−14𝑥²,(32) ARI can be approximated by ARI_IPFM≈14𝛽².(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𝑡_𝑘−𝑘𝑇₀=−𝑀¹𝜔_𝑚sin𝜔_𝑚𝑡_𝑘.(34) The deviations of the sampling times from a regular scheme is thus sinusoidal, with amplitude 𝑀₁/𝜔_𝑚, and the standard deviation estimated as the RMS value of a continuous-time sinusoid is 𝜎=1√2𝑀𝜔¹_𝑚.(35) Combining (33) with (35) and (26) results in the formula: ARI_IPFM≈2𝜋²𝑇𝜎₀².(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𝑇₀). 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𝑇₀. 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𝑇₀ 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.

(a)

(b)

(c)

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.

(a)

(b)

(c)

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 𝑀₁ 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).

(a)

(b)

(c)

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.

(a)

(b)

(c)

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 [16–18]. 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.

(a)

(b)

(c)

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.

(a)

(b)

(c)

(d)

(e)

(f)

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.