The electrocardiogram (ECG) contains well-known PQRST(U) pattern. Its morphology has been extensively studied since the technique of recording was introduced in the beginning of the twentieth century. The use of computers for ECG recording enabled studies of dynamic properties of the ECG. In particular, two important characteristics have been revealed: amplitude and frequency modulation.

The modulation of the cardiac cycle is extensively studied based on the signal that is generated from intervals between consecutives R-peaks, known as heart rate variability (HRV). Numerous approaches have been proposed and the signal have been decomposed in various ways—analyzing their spectral properties by frequency and time-frequency methods such as Fourier transform [2], autoregressive spectral estimation [11], wavelet transform [9, 17], or their scaling properties [6].

The amplitude characteristics and amplitude modulation received less attention, especially because the notion of the ECG amplitude is ambiguous: the heart as a source of the signal is moving and the signal is collected from the surface of the body that is moving as well. In addition, it is usually recorded as a 2-D projection. However, it has long been recognized that one of the sources of the amplitude modulation—the respiratory related movement of the thorax can be usefully used to extract the information about the respiration.

The influence of respiration on recording cardiac potentials was reported as early as in 1967; Flaherty et al. [5] and Moody et al. [13] were probably the first to discuss the possibility to derive the respiratory signal from the ECG. Since then several algorithms were proposed to reconstruct respiration by demodulating the ECG signal (e.g. [1, 4, 8, 10]).

The ECG-derived respiration signal either based on changes in heart rate or based on direct effect of respiration on the ECG waveform has been particularly used in studies of obstructive sleep apnoea [3, 16].

The problem to reconstruct waveform from the discrete Fourier analysis is that the biological signal has complicated properties; being non-stationary wave, including random noise and having interactive components. However, though people rarely notice, the essential problem is that there is no way to gather impartial and statistical samples. That is why samples from an unique sampling frequency are biased for spectral analysis even if fine sampling.

The purpose of this study is to introduce an algorithm to classify the multi-modulating waveforms on ECG signals. The algorithm is composed by the three particular scale transforms and the statistical Fourier transform (SFT) with impartial and statistical samples [12]. As a result, the multi-modulating source is decomposed by the SFT and is reconstructed into ECG, respiration and the other signals simultaneously by particular transformations.

The principle of this algorithm is the selection of three particular scale transforms. Two scales convert modulating waveforms into the periodic signal and the other converts the periodic components into a few terms. The role of the SFT is to hold sample information for any transforms.

Three scale transforms a_j(t), u_j(t), o_j(t) are adopted to the periodicity and to the reduction of Fourier components.

a_j = a_j(t) is a scale to match time intervals in the standard’s intervals and is called adaptive scale and u_j = u_j(τ_j) is a scale to match the magnitude spans in the standard’s magnitudes. Those scales convert the source x′(t) into the standard x_j(t) as x_j(t) = x′(a_j)/u_j(a_j). a_j(t) and u_j(a_j), then, are defined as a join of each a_kj(t) and u_kj(a_kj) (k = 1,..., k_j·) when the standard repeats k_j times in a window T_j.

o_j = o_j(t), called inherent scale, is a timescale to convert the standard x_j(t) into the Fourier series of r_j components passing through 2r_j poles of x_j(t). The scale transform τ_j(t) is defined as the composite transform as follows (Appendix 2):

1

The transformed signal (x′(τ_j) − ε_j(τ_j))/u_j(τ_j) occupies r_j instantaneous components of the skipped ik_jth term in a window T_jincluding k_j modulations. While, instantaneous components of the ε_j(τ_j)/u_j(τ_j) scatters into broad components from the 1st to the nth terms. As a result, (x′(τ_j) − ε_j(τ_j))/u_j(τ_j) fits the scales τ_j, u_j and ε_j(τ_j)/u_j(τ_j) unfits scales τ_j, u_j, and then, the source x′(t) is decomposed into two class modulations of x′(τ_j) − ε_j(τ_j) and ε_j(τ_j) (Appendix 3).

That is, x′(τ_j) − ε_j(τ_j) is extracted by v_j(τ_j) = Σu_j(τ_j)h_ij(τ_j) from the ik_jth (i = 1,.., r_j) components and ε_j(τ_j) is extracted by the other ε_j(τ_j) = Σu_j(τ_j)h_ij(τ_j) except the ik_jth components when n instantaneous components of x′(τ_j)/u_j(τ_j) are defined as h_ij(τ_j). The source x′(t), therefore, is decomposed into x′(t) − ε_j(t) and ε_j(t) by inverse scale transform t = t(τ_j) under an adequate u_j(τ_j).

The magnitude scale u_kj(τ_j) is defined as an expression of the magnitude pattern in each window T_kj when the signals v_j(τ_j)(=x′(τ_j) − ε_j(τ_j)), ε_j(τ_j) and u_j(τ_j) are expressed by {v_kj(τ_j)}_k, {ε_kj(τ_j)}_k and {u_kj(τ_j)}_k (k = 1,.., k_j), and each changes its waveform independently. Especially, v_j(τ_j) is expressed by a constant value u_j(τ_j) = u_kj in each window T_kj when its magnitude modulation takes a similar figure with the standard (Appendix 4).

When x′(t) is m classes of Fourier series with terms of r₁, r₂,..., r_m (n = r₁ + r₂ +   + r_m), it expresses Eq. 2.

2

In Eq. 2, t = t(τ_j) is the inverse scale transform of the jth scale τ_j to scale t, v_j(t) is a Fourier series of the jth class modulation, and ε_m(t) is the mth residual signal. Equation 2 shows that x′(t) is reconstructed completely by periodic Fourier series v_j(τ_j) on scale τ_j (j = 1,2,..., m).

Since the probability of components on the jth class modulation is to be p_j = r_j/n, the probability to encounter x′(t) with m modulations is p(r₁,r₂,..., r_m) = p₁^r1, p₂^r2,..., p_m^rm which relates to the geometrical mean . Equation is as follows (Appendix 5).

3

It clearly shows that the quantity of information on multi-modulating signal is m itself.

When h_ij(τ_j) is the ith instantaneous Fourier component of standard x_j(τ_j) in scale τ_j and h_ij(t) is its modulating component in scale t, the morphological relation is given by the following Formula (4):

4

In Eq. 4, τ_j/t is the Jacobian. Therefore, formula (4) shows that the modulation of every component at scale t is expressed by scale’s modulations on its magnitude u_j(τ_j) and its time τ_j(t) = o_j(a(t)). Especially, the adaptive scale a_j(t) describes composite transform a_j(t) = c₁(c₂( c_n(t))) and that of exchanged orders a_j(t) = c_n( c₂ (c₁ (t))), so the sth scale c_s(t) will express its specific modulation from the standard such as modulations of PQ, QRS, QT and RR.

The analysis of scale classification of the three modulating signals on V4 lead was carried out by using the statistical Fourier analysis with maximum terms of 256 with A/D converted 16-bit samples of the sampling frequency 20 kHz over successive 50 s. To gather independent samples with 256 terms, samples needed to be accumulated by the high speed 20 kHz of A/D converter with a large dynamic range of 16 bit to avoid the deform of a tiny respiratory signal (Fig. 1).

Fig. 1 Synchronous frequency samplings on scale t and τ

The ECG standard x_j(t) was selected from the source x′(t) without breathing. On the other hand, x′(t) with natural breathing and the respiratory flow are recorded simultaneously by using Number 3 Fleisch pneumotachograph. Both the ECG and the respiratory scale transforms were the same as in Fig. 2a, b of Appendix 2. The ECG x′(t) − ε_j(t) and its baseline wander ε_j(t) were analyzed when u_j(τ_j) = 1 and the respiratory standard x_j+1(t) was selected from stationary sequences of ε_j(t) in 16 respiratory cycle. (ε_j(τ_j+1) − ε_j+1(τ_j+1))/u_j+1(τ_j+1) and ε_j+1(τ_j+1)/u_j+1(τ_j+1) were analyzed when u_j+1(τ_j) of the formula (5) in Appendix 4 and lung ventilation V_j+1(t) and its residual ε_j+1(t) were determined.

Fig. 2 ECG and respiratory standard waves and specific points on scale t. a Cardiac standard, b respiratory standard, and c detecting dots on a part of the source

Signals corresponding to every process on the scale classification of m = 2 under constant magnitude scale of u_j = 1 are illustrated in Fig. 3 with five cardiac periods. That is, intervals among specific points in source x′(t) in Fig. 3a are transformed into the same intervals of the standard signal with the adaptive scale a_j(t) to x′(a_j(t)) of Fig. 3b and with the inherent scale τ_j(t) to x′(o_j(a_j(t))) of Fig. 3c. The transformed signal x′(τ_j) is decomposed into statistical Fourier components on a fundamental period of 53 cardiac intervals. The transformed signal x_j′(τ_j) in Fig. 3d is a reconstructed signal from the 53rd, 106th, 159th and 212th components, and its residual ε_j(τ_j) in Fig. 3e from all other components. The extracted ECG x_j′(t) in Fig. 3f and the base line wander ε_j(t) in Fig. 3g are the inverse scale-transformed signals from scale τ_j to scale t. When tiny-leak signals on the ε_j(t) synchronized with intervals of QT are replaced with linear trends, the figure of pattern ε_j(t) takes that of the lung ventilation in Fig. 3i. The smoothing signal, V_j(t) in Fig. 3h is a filtered signal of ε_j(t) with low-pass filter (LPF) under 30 terms (0.6 Hz). V_j(t) is illustrated with enlarged four times. The signal V(t) in Fig. 3i is the lung ventilation curve which is an integral signal of the direct measurement of the respiratory flow

Fig. 3 Classified process of modulating signal of the algorithm

Any differences of P ~ U waveforms between the source x′(t) and reconstructed ECG x′_j(t) bear getting information, and also the differences between ventilation V(t) and baseline wander V_j(t) are small when frequencies lower than breathing are neglected. x′_j(t) provides time-interval information on Minnesota codes better than that of x′(t) clearly in visual.

Correlation coefficients ρ between the signal V_j(t) and V(t) are calculated by each ventilation with scatter plots in Fig. 4a, and are classified by P ~ U waves of 53 cardiac signals in Table 1a and by respiratory phases of 16 ventilations in Table 1b. Total ρ takes even values (ρ ≈ 0.82) and has less interference between P ~ U waves in cardiac cycles. However, it takes low values (ρ ≤ 0.5) at the transient regions of I, III between the inspiration and the expiration corresponding to loops in Fig. 4a. The inspiration and the expiration except transient regions are statistically dominant at higher value (ρ ≈ 0.95). It is clear that the scale classification extracts the electrical activity of the heart with constant magnitude scale of u_j(τ_j) = 1 in hole window of T_j. The extracted baseline wander is constructed with lung ventilation and other biological information.

Fig. 4 X–Y loop among modulating Vj(t), Vj+1(t) and V(t)

Table 1 Correlation coefficients ρi among Vj(t), Vj+1(t) and V(t) on cardiac and respiratory stages

The magnitude scale transform is applied to the scale classification of the base line wander V_j(t) with a respiratory standard in Fig. 2b. Scale classification of m = 4 is illustrated in Fig. 5, the source x′(t) in Fig. 5a, the extracted ECG x_j(t) in Fig. 5b, the extracted respiratory signal V_j+1(t) in Fig. 5c, the residual ε_j+1(t) in Fig. 5d, the leak ECG ε(t) over 0.6 Hz of ε_j(t) in Fig. 5e and the ECG of x_j(t) + ε(t) in Fig. 5f. The signal V(t) in Fig. 5g is 16 measured respiratory ventilations. The magnitude scale Δu_j+1(t) is illustrated in Fig. 5h.V_j+1(t) and ε_j+1(t) are displayed in the reverse because the extracted respiratory waveform is a visually reversed V(t). ε_j+1(t) is similar to base-line trends of V(t) which relate to the thoracic gas volume (TGV).

Fig. 5 Scale classification of the four modulating signals in an electro-cardiac signal

The scatter plots of 15 lung ventilations between V_j+1(t) and V(t) take the similar pattern as illustrated in Fig. 4b. Correlation coefficients ρ in Table 1c takes higher value than the others though transient regions change I, III into II, IV.

The frequency distributions of the amplitude components on |V(f)|, |V_j+1(f)|, |ε_j+1(f)| and that of phase shifts Δ(f) between V(f) and V_j+1(f) are illustrated in Fig. 6. On frequency ranges from 0.2 to 0.5 Hz, amplitude distribution in |V(f)| and |V_j+1(f)| takes similar pattern and those of phase shifts are invariant near 0°. Amplitude distribution in |ε_j+1(f)| has not only the lower frequency components under 0.06 Hz corresponding to the trends of TGV but also the leak components of the lung ventilation from 0.25 to 0.5 Hz. The components of |ε_j+1(f)|, moreover, has the harmonics among 0.04, 0.1, 0.18 and 0.28 Hz though those harmonics disappear in |V(f)| and |V_j+1(f)| at least.

Fig. 6 Frequency distributions of the measured V(f) and the extracted Vj+1(f), εj+1(f)

As a result, the time and magnitude scale is effective to classify modulating signals from the higher frequency V_j+1(t) to the lower frequency ε_j+1(t) because a residual ε_j+1(t) has less interference of the higher V_j+1(t) and successive sequential ε_j+1(t) will supply the signal for detecting other lower modulation.

A discrepancy between time domain and frequency domain analysis happens commonly by using biased samplings for discrete signal processing in various way such as Hilbert transform [2] and wavelet analysis [9] relating with an unique sampling frequency. Moreover, the mathematical concepts of one-to-one correspondence, n unknown amplitudes and phases are derived from 2n samples, hinders to orientate the statistics in the signal processing.

The SFT, however, has no contradiction between Fourier analysis and the statistics which resolute components with the same signal-to-nose ratio (SNR) for any signal decreasing amplitude components with order 1/f such as ECG signal, which is unsatisfied in the autoregressive spectral estimation [11]. The SFT component, moreover, is given independently by a time trace of instantaneous components at each instant as shown in Appendix 1. That is, its amplitude and phase are not calculated by the inner products between the source and trigonometric bases but by the extracted traces of instantaneous sin-wave themselves [12]. The concepts of signal arrangements at series {t_k} is another thing on the SFT because the replacement between t_i and t_j needs the replacements of all the synchronous frequency samplings between t_i and t_jto hold the same instantaneous components. The summation of replacements, then, is about . In other words, the SFT supplies a signal ensemble arranging samples of the series {t_k}. In this study, discrete samples from A/D converter of 20 kHz satisfy the independency as shown in Appendix 6. As a result, the random noise including a signal is not analyzed as a white noise but as inherent components characterizing arrangements. Therefore, the reconstructed waveform of a baseline wander holds a prototype even if it includes many random noises.

The necessary condition of the spectral analysis known as Nyquist’s condition [14, 15] becomes precise on the SFT as shown in Appendix 3, which confirms its condition statistically. The potentiality to scatter the aliased components impartially discriminates the waveforms under the unsatisfied Nyquist’s condition. In this algorithm, while the inherent scale transform converges to the standard x_j(τ_j) with four skipping components, others ε_j(τ_j), either continuous or random, disperse components into 256 terms.

Since the scale-transformed samples on the SFT are gathered from a whole window at any instant, those samples reflect its scale transform better than any local one. The random and leak signal on the extracted ECG x_j′(t) are reduced into the tiny level under a stable behavior of the random noise. The 256 terms is enough to analyze scale-transformed signal in a window under 53 cardiac cycles because transformed components with the ranges, 8 ≤ k_j ≤ 53, takes almost the same figures. The high-resolution analysis, therefore, is accomplished by the relative high SNR from those instantaneous components.

The problem on the scale transforms is to be adequate scales, not to be expressing the true biological modulation. a_j+1(t) is, however, selected a traditional scale for non-stationary sequences in each respiratory cycle which carry information stored in Fourier phases as already indicated in [6] and a_j(t) is selected from cardiac cycles without breathing in addition to the non-stationary sequences. As a result, o_j(a_j) and o_j+1(a_j+1) make those mutual relationship of phases more simple relations. Taking own standard characterized with number of poles is also reasonable to determine the waveform reconstruction and to discriminate the own HRV. Moreover, the lost waveform is scattered into the other components without losing the morphologic information such as the leak ECG waveform.

It is well known the volume conductivity of the thorax takes inequality and changes the electro cardiac potentials on the thorax’s surface with respiratory movements [1, 5].

In the scale classification of m = 4, it indicates that a constant magnitude scale u_kj = 1 produces unchangeable heart-activity x_j′(t) so that the signal ε_j(t) appears in respiratory modulations from their changes of the conductivity. It suggests that the standard without breathing follows the formula (2) of the signal classification and the source follows the formula (4) which expresses the modulating standard with Jacobian, τ_j/t. A composite scale c_s(t), therefore, is expected to express the inherent modulation from the standard such as PQ, QT and RR. That is, the spectral analysis of c_s(t) on RR provides the distribution of the HRV without a tachogram which is expressed in other ways by others [2, 9, 11]. The extracted ECG has the potentiality to detect R–R modulations of the heart itself.

The algorithm extracts the lung ventilation curves {V_j+1(t)}_k from V4 lead as changing the thoracic gas volume [1, 4]. The difference of amplitude distribution between V(t) and V_j+1(t) from spectral analysis takes the same pattern around 0.3 Hz [8]. The phase difference between V(t) and V_j+1(t), however, occurs in the transient regions between the inspiration and the expiration which deduce from inadequate respiratory standards, false respiratory peaks and the gas compression in the thorax [7]. It is effective to determine specific components, such as harmonics of ε_j+1(t) around (0.04, 0.1, 0.18 and 0.28 Hz) in cardiovascular responses [2, 17].

Information of the partial ventilation is also expected from other 12 leads [8]. The respiratory Jacobian τ_j+1/t, like that of ECG, has the potentiality to detect respiratory modulation itself as a part of the HRV and u_j+1(t) is expected to detect the amplitude modulation such as sighs distribution. The residual ε_j+1(t) supplies information of waveforms not only of the thoracic gas volume but also of the blood flow relating to cardiovascular responses [2, 11, 17].

It can be concluded that this algorithm reveals the simultaneous measurements of ECG and respiration. The comparison with the frequency distribution and the waveform evaluates new modulations clearly, such as u_kj+1 for detecting changes of the tidal volume (TV) and for finding TGV and cardiovascular responses. This contribution of the algorithm indicates diagnostic potentials on the sleep apnoea from ECG recordings [16].

Instantaneous Fourier components on the SFT at any instant t_k are derived from samples of the synchronous frequency samplings from the source, having maximum components number of n, as illustrated in Fig. 1a. Those samples at any instant t_k are independent of each other when the t_k is given by t_k = kT_j/s (k = 0,..., s − 1) in a window T_j for a prime number s of s > 2n. Although the sampling timing from the transformed source at scale τ_j is the same pattern in Fig. 1a, timing at scale t becomes as illustrated in Fig. 1b. Any source sample, therefore, is expressed by own components on each scale independently.

The ECG standard x_j(t) and its scale-transformed signal x_j(τ_j) are illustrated in Fig. 2a and those of the respiratory x_j+1(t) and x_j+1(τ_j+1) are illustrated in Fig. 2b.x_j(τ_j) is four terms of Fourier series passing eight poles of x_j(t). x_j+1(τ_j+1) is a sin wave passing two poles of x_j+1(t). Inherent scale τ_j = o_j(t) and τ_j+1 = o_j+1(t) are defined as transforms from x_j(t), x_j+1(t) into x_j(τ_j), x_j+1(τ_j).

The detected specific points for the source and the ECG standard are illustrated in Fig. 2a, c as a dot curve. These dots indicate start, pole and end of the P ~ U waves in each cardiac cycle. Specific points of lung ventilation are the poles of inspiration, expiration, midpoints and expiratory inflection point in each respiratory period. The midpoints represent respiratory maximum flows and expiratory inflection points distinguish expirations relating to pleural pressure. Points are determined to divide time interval of the specific-points into 6 ~ 16 to approximate true signals. Both the ECG and the respiratory standard are line curves of connected points. The adaptive scales a_j(t), a_j+1(t) are coordinated by divided points after adjusting detected specific points to those of the standard. The scale-transformed signal x′(τ_j) and that of the residual ε_j(τ_j) at scales τ_j, τ_j+1 are determined from formula (2), whereas the scale transforms are τ_j(t) = o_j(a_j(t)) and τ_j+1(t) = o_j+1(a_j+1(t)).

The inherent transform o_j(t) satisfies the necessary condition of the SFT analysis because Fourier components appear in the skipped ik_jth (i = 1,..,r_j) terms under n > r_jk_j when x′(τ_j)/u_j(τ_j) is periodic waves of k_j standard. Contrarily, instantaneous components are scattering into broadband of 1st ~ nth components when the source x′(t) unfits scales τ_j, u_j.

The necessary condition of the SFT analysis is that the mean m_i(t_k)(i = n + 1,..., 2n) over the (n + 1)th of synchronous frequency sampling takes the same value at any instant t_k as m_n+1(t_k) = m_n+2(t_k) =   = m_2n(t_k) which has the same properties with Nyquist’s condition. Each instantaneous component includes errors of aliased components under the unsatisfied condition, but the summation of all the instantaneous components comes back to a source sample completely at any instant and the SNR of each component is improved with an adequate n value because the aliased errors are scattered into n components monotonously than those of the unique-interval samplings.

The Δu_kj(=u_kj − 1) is given by the following formula (5) when x_k′(τ_j) = v_kj(τ_j) + ε_kj(τ_j) under u_kj(τ_j) = 1.

5

The constant component of each window T_kj is compensated under the continuous condition with an end expiratory point at each edge of window T_kj.

The cases of modulations take when the Fourier components of n = r₁, + r₂, +   + r_m are distributed into m cases. The probability p(r₁, r₂,..., r_m) to encounter x′(t) with m modulations is given by Starling’s formula as below.

When the random signal in Fig. 7A-b is a permuted source in Fig. 7A-a, the amplitude distributions on the SFT are illustrated in Fig. 7B and on the direct Fourier transform (DFT) in Fig. 7C, respectively. While amplitude distributions of the source take the same forms between the SFT in Fig. 7B-a and the DFT in Fig. 7C-a, those of the random signal are different in Fig. 7B-b and C-b. That is, the distribution takes an exponential form on the SFT but takes an even form on the DFT which is well known as the form of a white noise.

As amplitude and phase components of the SFT reconstruct the random signal in Fig. 7A-b completely, the inverse rearrangement of the reconstructed random noise comes back to the source in Fig. 7A-c. Contrarily, the waveform from the DFT is still a random in Fig. 7A-d not the source. As a result, the distribution of samples is kept in SFT but lost in DFT. The SFT supplies true instantaneous Fourier components for any scale transform including arrangements. In other words, a signal with any random noise holds far more stable information in the SFT than in any other integral transform.