# Observability of Spectral Components beyond Nyquist Limit in Nonuniformly Sampled Signals

Authors: Jozef Púčik; Oldřich Ondráček; Elena Cocherová
Journal: ISRN Signal Processing
 (2012)

---

## 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.


---

## Body

## 1. Introduction

Aliasing is a well-known phenomenon in signal theory that occurs if a continuous-time signal is sampled at a frequencyfs below twice the maximum frequency fm of the sampled signal, that is, fs<2fm. The sampling frequency halved, fN=fs/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-fs/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/(2tmin), where tmin is the smallest time interval between data, while others infer the Nyquist frequency from average sampling times tmean and  defines it as 1/(2tmean). 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. Section2 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)P(ω)=1N|∑n=0N-1x(tn)e-jωtn|2,
where x(tn) are data samples taken at the time instants tn. Notice that roots of a Fourier-type analysis with complex exponentials reside in a least-squares fitting of data to the complex sinusoid ejωt. Lomb [2] formulated the periodogram by means of the least-squares fitting of observed real-valued data to the model of the form:
(2)Acos⁡⁡(ω(ti-τ))+Bsin⁡⁡(ω(ti-τ)).
The redundant time delay variable τ, computed as
(3)tan(2ωτ)=∑n=0N-1sin⁡⁡(2ωtn)∑n=0N-1cos⁡⁡(2ωtn),
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]
(4)PLS(ω)=12{(∑n=0N-1x(tn)cos⁡⁡(ω(tn-τ)))2∑n=0N-1cos⁡2⁡(ω(tn-τ))PLS+(∑n=0N-1x(tn)sin⁡⁡(ω(tn-τ)))2∑n=0N-1sin⁡2⁡(ω(tn-τ))}.
Provided that the summed sine and cosine squared terms in (4) are both equal and approach to N/2, and x 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 amplitudeC, frequency ω0, sampled at the time instants tn:
(5)x(tn)=Cejω0tn.
The sampling instants are supposed to be the random numbers obeying the following model:
(6)tn=nT0+ξn,
where T0 is the mean sampling period, n=0,…,N-1, and ξn 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):
(7)X-(ω)=E[X(ω)]=E[∑n=0N-1x(tn)e-jωtn].
After substitution (5), (6) into (7) and rearrangement, we obtain
(8)E[X(ω)]=C∑n=0N-1e-j(ω-ω0)nT0E[e-j(ω-ω0)ξn].
In (8), we can recognize FT of a uniformly sampled sinusoid expressible by the periodic Dirichlet kernel:
(9)NDN((ω-ω0)T0),
where the kernel is
(10)DN(x)=1Nsin⁡⁡(Nx/2)sin⁡⁡(x/2).
The mean FT (8) can be rewritten in terms of the characteristic function shifted by the frequency ω0:
(11)X-(ω)=CNDN((ω-ω0)T0)T(ω-ω0),
where the characteristic function T of the timing irregularity ξ, characterized by the probability density function pξ, is defined as
(12)T(ω)=E[e-jωξ]=∫-∞∞pξ(ξ)e-jωξdξ.
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:
(13)P-(ω)=E[P(ω)]=1NE[(∑n=0N-1x(tn)e-jωtn)(∑k=0N-1x(tk)e-jωtk)*]=1NE∑n=0N-1∑k=0N-1x(tn)x*(tk)e-jω(tn-tk).
The summation term in (13) after substitution (5) is
(14)C2NE∑n=0N-1∑k=0N-1ejω0(n-k)T0e-jω(n-k)T0e-j(ω-ω0)ξnej(ω-ω0)ξk.
For n≠k, we have due to the presumed independence:
(15)E[ejωξne-jωξk]=T(-ω)T(ω)=T2(ω),
where we have assumed a symmetrical probability function in (12) and thus a real-valued characteristic function.In the specific case,n=k, expectation (15) is 1 and it can be formally expanded as
(16)E[ej0]=T2(ω-ω0)+[1-T2(ω-ω0)].
After substitution m=n-k in (14) and taking into account (15), (16), the summation (14) simplifies to
(17)P-(ω)=C2[1-T2(ω-ω0)]+C2T2(ω-ω0)∑m=-N+1N-1(1-|m|N)ejω0mT0e-jωmT0.
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
(18)SN(ω-ω0)=NDN2((ω-ω0)T0).
The final formula then becomes
(19)P-(ω)=T2(ω-ω0)C2SN(ω-ω0)+C2[1-T2(ω-ω0)].
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 T(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 (Figure2). The beat occurrence times (representing sampling times) tk generated by the IPFM model can be defined as [12]
(20)k=∫t0tk1+m(t)T0dt,
where k is an integer, t0=0, T0 is the mean heart period, and m(t) is a dimensionless modulating signal interpreted as the fractional change of the instantaneous heart rate relative to 1/T0.Figure 2
IPFM model for generating the beat occurrence times.The sampling pattern can be written as a sequence of the Dirac pulses:(21)spc(t)=∑k=-∞∞δ(t-tk).
The Fourier transform of (21) is referred to as the spectrum of counts (SPC) [12]:
(22)SPC(ω)=∑k=-∞∞e-jωtk.
The spectrum of a sampled signal Xs(ω) is related to the true spectrum X(ω) by the frequency-domain convolution:
(23)Xs(ω)=12πX(ω)*SPC(ω).
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
(24)SPC(ω)=1T0(2πδ(ω)+M(ω))*{δ(ω)+1π∑n0=1∞FT{cos⁡⁡(2πn0T0(t+∫-∞tm(τ)dτ))}}.
Interestingly, the spectrum of the modulating signal M(ω) is directly included in (24). Beside the Dirac pulse located at the frequency origin, infinite numbers of frequency modulated (FM) terms convolved with M(ω) 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:
(25)m(t)=M1cos⁡⁡(ωmt),
where ωm 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 J0, the argument of which is given by the gradually increased index of the frequency modulation:
(26)β=n0ωsM1ωm,
where ωs is the mean sampling frequency corresponding to the first carrier (n0=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 rangefs/2 to fs. 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-fs/2. Specifically, if an input signal has frequency fi,fN<fi<fs, the aliased component appears at the frequency fa=fs-fi. If the spectrum is plotted in a range above fs/2, both the spectral peaks, fi and fa, 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):
(27)ARI=Ai-AaAi,
where Ai is the spectral amplitude at the true (input) frequency fi, Aa is the spectral amplitude at the corresponding aliased frequency fa. 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 ω=+ωi=ω0, while the component aliased into 0-fs/2 is originated from the negative input frequency -ωi, shifted by value ωs, that is, -ωi=ω0 and ω=-ωi+ωs. The ARI can be then written in terms of the characteristic function as
(28)ARIRND=T(0)-T(ωs)T(0).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:
(29)T(ω)≈1+12!T′′(0)ω2=1-12σ2ω2.
Using ωs=2π/T0 in (29), T(0)=1 in (28), we obtain an approximate formula for the ARI in term of the relative timing standard deviation σ/T0:
(30)ARIRND≈2π2(σT0)2,
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 +ωi and by -ωi. A dominant aliased term is associated with the 1st FM carrier, that appears, after -ωi shifting by ωs, at the frequency ωa=ωs-ωi. Its amplitude is determined by the 0th-order Bessel function of the first kind for β given by (26), n0=1. The resulting formula for the ARI estimation, analogical to (28), becomes
(31)ARIIPFM=J0(0)-J0(β)J0(0).
Since J0(0)=1, and using expansion of J0:
(32)J0(x)=∑k=0∞(-1)k(k!)2(14x2)k≈1-14x2,
ARI can be approximated by
(33)ARIIPFM≈14β2.
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(34)tk-kT0=-M1ωmsin⁡⁡(ωmtk).
The deviations of the sampling times from a regular scheme is thus sinusoidal, with amplitude M1/ωm, and the standard deviation estimated as the RMS value of a continuous-time sinusoid is
(35)σ=12M1ωm.
Combining (33) with (35) and (26) results in the formula:
(36)ARIIPFM≈2π2(σT0)2.
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 ratefs. 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.65fs which is aliased to the frequency 0.35fs. Conversely, an input frequency 0.35fs below the Nyquist limit 0.5fs will be imaged to 0.65fs 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-fs instead of a conventionally used range 0-fs/2.Figures3 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.05T0). 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.122T0. 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.3T0 ensures that randomly chosen proximate sampling times tk-1,tk, fulfill the monotonicity condition tk-1<tk 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.Spectral amplitudes of the LS spectra—nonuniformly sampled signal with the random sampling pattern. The input frequency 0.35fs lies below the Nyquist limit. Different levels of sampling irregularity: (a) low, (b) medium, (c) high.
(a)
(b)
(c)Spectral amplitudes of the LS spectra—nonuniformly sampled signal with the random sampling pattern. The input frequency 0.65fs exceeds the Nyquist limit. Different levels of sampling irregularity: (a) low, (b) medium, (c) high.
(a)
(b)
(c)A similar set of simulations (Figures5 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 M1 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).Spectral amplitudes of periodograms—nonuniformly sampled signal with IPFM generated sampling pattern. Input frequency 0.35fs lies below the Nyquist limit. Different levels of sampling irregularity: (a) low, (b) medium, (c) high.
(a)
(b)
(c)Spectral amplitudes of periodograms—nonuniformly sampled signal with IPFM generated sampling pattern. Input frequency0.65fs exceeds the Nyquist limit. Different levels of sampling irregularity: (a) low, (b) medium, (c) high.
(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 Figure7(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.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.
(a)
(b)
(c)The results of application of different methods for HRV spectrum computation are summarized in Figure8. 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.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.
(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.



---

*Source: 643563-2012-06-28.xml*