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, . The sampling frequency halved, , 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 , 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 , where is the smallest time interval between data, while others infer the Nyquist frequency from average sampling times and defines it as . 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]: 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: The redundant time delay variable , computed as 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] Provided that the summed sine and cosine squared terms in (4) are both equal and approach to , 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 : The sampling instants are supposed to be the random numbers obeying the following model: where is the mean sampling period, , 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): After substitution (5), (6) into (7) and rearrangement, we obtain In (8), we can recognize FT of a uniformly sampled sinusoid expressible by the periodic Dirichlet kernel: where the kernel is The mean FT (8) can be rewritten in terms of the characteristic function shifted by the frequency : where the characteristic function of the timing irregularity , characterized by the probability density function , is defined as 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: The summation term in (13) after substitution (5) is For , we have due to the presumed independence: 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 After substitution in (14) and taking into account (15), (16), the summation (14) simplifies to 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 The final formula then becomes 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 , 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.
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] where is an integer, , is the mean heart period, and is a dimensionless modulating signal interpreted as the fractional change of the instantaneous heart rate relative to .
The sampling pattern can be written as a sequence of the Dirac pulses: The Fourier transform of (21) is referred to as the spectrum of counts (SPC) [12]: The spectrum of a sampled signal is related to the true spectrum by the frequency-domain convolution: 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 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: 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: where is the mean sampling frequency corresponding to the first carrier .
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 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 . Specifically, if an input signal has frequency , the aliased component appears at the frequency . If the spectrum is plotted in a range above , 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): 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 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
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: Using in (29), in (28), we obtain an approximate formula for the ARI in term of the relative timing standard deviation : 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), . The resulting formula for the ARI estimation, analogical to (28), becomes Since , and using expansion of : ARI can be approximated by 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 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 Combining (33) with (35) and (26) results in the formula: 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 instead of a conventionally used range .
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 . 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 . 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 ensures that randomly chosen proximate sampling times , fulfill the monotonicity condition 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)
(a)
(b)
(c)
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)
(a)
(b)
(c)
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)
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)
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.