1092
IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 44, NO. 11, NOVEMBER 1997
Spectral Decomposition in Multichannel Recordings Based on Multivariate Parametric Identification Giuseppe Baselli,* Alberto Porta, Ornella Rimoldi, Massimo Pagani, and Sergio Cerutti, Senior Member, IEEE Abstract—A method of spectral decomposition in multichannel recordings is proposed, which represents the results of multivariate (MV) parametric identification in terms of classification and quantification of different oscillating mechanisms. For this purpose, a class of MV dynamic adjustment (MDA) models in which a MV autoregressive (MAR) network of causal interactions is fed by uncorrelated autoregressive (AR) processes is defined. Poles relevant to the MAR network closed-loop interactions (clpoles) and poles relevant to each AR input are disentangled and accordingly classified. The autospectrum of each channel can be divided into partial spectra each relevant to an input. Each partial spectrum is affected by the cl-poles and by the poles of the corresponding input; consequently, it is decomposed into the relevant components by means of the residual method. Therefore, different oscillating mechanisms, even at similar frequencies, are classified by different poles and quantified by the corresponding components. The structure of MDA models is quite flexible and can be adapted to various sets of available signals and a priori hypotheses about the existing interactions; a graphical layout is proposed that emphasizes the oscillation sources and the corresponding closed-loop interactions. Application examples relevant to cardiovascular variability are briefly illustrated. Index Terms— Biomedical signal processing, cardiovascular variability signals, multivariate parametric models, multivariate spectral decomposition.
I. INTRODUCTION
T
HE ANALYSIS of biological signals often requires the comparison of multiple recordings which are differently affected by the same oscillation sources [1]–[4]. Parametric spectral analysis [5], [6] permits only the recognition and quantification of the oscillatory components in the single signals [7]; on the contrary, multivariate (MV) parametric identification [8]–[11] provides further information about the causal interactions among the signals [12], [13] and about cross-spectral patterns [14]–[16], but it is difficult to use it for a better comprehension of the rhythmogenetic mechanisms. So, it can be worth presenting the results of MV analysis in a way immediately contiguous to the more familiar one of spectral analysis, in order to understand whether a rhythmic component impinges on a specific channel and propagates to the other ones or whether several interacting channels are involved in its genesis and also to understand how the interactions between the different channels amplify or damp oscillations. Manuscript received March 5, 1995; revised April 16, 1997. Asterisk indicates corresponding author. *G. Baselli is with the Dip. di Elettronica per l’Automazione, Universit`a di Brescia, Via Branze 38, Brescia 25123 Italy (e-mail:
[email protected]). A. Porta and S. Cerutti are with the Dip. di Bioingegneria, Politecnico di Milano, Milan 20133 Italy. O. Rimoldi and M. Pagani are with Medicina II, Ospedale “L.Sacco,” Universit`a di Milano, Milan 20157 Italy. Publisher Item Identifier S 0018-9294(97)07417-X.
The present work proposes a method of spectral decomposition and of rhythm classification based on MV analysis. Rhythm sources are classified considering the closed-loops which can resonate and produce oscillations. Consequently, the spectral decomposition quantifies the classified oscillations in each signal. Therefore, the method is based on the identification of a quite general class of MV dynamic adjustment (MDA) models in which closed-loops of simple structure are evidenced both in their deterministic interactions and in their stochastic inputs. Signal interactions are described by causal finite impulse response (FIR) filters, forming a MV autoregressive (MAR) network which is driven by monovariate and uncorrelated autoregressive (AR) random inputs. A decomposition into spectral components relevant to the poles of the MAR network and to those of the single AR inputs is, consequently, derived and an useful classification of spectral components is obtained. This approach was developed and first applied for the study of the origin of rhythms in heart rate and arterial pressure beat-to-beat variability signals, in which low-frequency [(LF), around 0.1 Hz] and high-frequency [(HF), at respiratory frequency around 0.3 Hz] oscillations are expressive of the control action of the autonomic nervous system [17]–[19]. Nonetheless, in the present paper the class of models is generalized in order to introduce a method of wide application illustrated in Section II. The flexibility of the procedure is based on the possibility of choosing the model structure best suited to describe the signal interactions and consequently the rhythm sources in the particular problem under study. The MDA class of models includes monovariate AR, ARXAR, MAR models and more complicated structures with exogenous inputs, AR stochastic inputs and various interaction loops. Examples of simple model layouts and their graphical representation are illustrated in Section III-A. Section III-B presents an analysis example on cardiovascular variability data using three different model structures (separate monovariate AR, bivariate AR, and specific double-loop structure). This example is provided in order to illustrate the effects of model choice and show diagrams useful in presenting the analysis results. II. METHODS A. The Class of Models The general structure of the considered class of MDA models is that of a MAR network fed at the level of each signal, by single-channel AR processes, , which are supposed uncorrelated from one another. Therefore, it can be written
0018–9294/97$10.00 1997 IEEE
(1)
BASELLI et al.: MULTIVARIATE SPECTRAL DECOMPOSITION
1093
Fig. 1. General structure of a MDA model with X independent AR exogenous inputs and (L-X ) interacting signals. (L-X ) independent AR stochastic inputs u and the exogenous inputs feed the MAR network of interacting signals. So L+1 loops are evidenced which, when resonating, define different rhythm sources.
m
where is the vector of signals and is the matrix of the causal FIR filters which describe the interactions among the signals. Elements on the diagonal are of the type (2) i.e., are strictly causal FIR filters and represent the regression of on samples of its own past. The other elements are causal FIR filters which describe the effect (exhausted within samples) on a signal of a different signal As for elements outside the diagonal, the notation , equivalent to , is introduced (and used also in the graphical representations) by analogy of the regression of on in single output models
(3) The immediate (i.e., not delayed) effects described by the coefficient are present only in some interactions and determine the canonical form of the MV process (see below). Any element of is described as a single-channel AR process, therefore (4) is a strictly causal FIR filter and where of the vector of uncorrelated white noises
is an element
(5) Any can be supposed is the all-pole spectral factor of white and the relevant set to zero; when all are white the MDA model is reduced to a MAR one.
The joint process form of a MV model is and describes the signals as the product of the white noise vector filtered by the matrix of rational transfer function From (1) and (4), MDA models are characterized by (6) where (7) is the transfer function matrix of the MAR network and contains the spectral factors of the uncorrelated inputs. The general black-box structure of MDA models given by (1) and (4) can be easily adapted to particular applications and include a priori notions about existing interactions and rhythm sources; so gray-box models for specific signal analysis problems are obtained. The absence of any interaction is accounted for simply by substituting the corresponding element with a zero in matrix In other cases, an element of can be constrained to have a more simple structure (e.g., a pure gain) compared to the others, simply by limiting its order. In addition, the possibility of a given residual to introduce rhythms can be excluded omitting . Drawings, such as those shown in Fig. 2 and described in Section III-A, display the obtained layout and emphasize closed loops as oscillation sources. It can often be assumed that some of the signals are exogenous, i.e., they do affect the other signals without being affected. Let us arrange vector so that the last signals, , are exogenous and uncorrelated, while the first ones, , are
1094
IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 44, NO. 11, NOVEMBER 1997
(a)
(b)
(c)
(d)
(e)
(f)
Fig. 2. Examples of simple structures pertaining to the class of MDA models: (a) monovariate AR model; (b) ARXAR model; (c) bivariate MDA model with single loop and two AR stochastic inputs; (d) bivariate MAR model; (e) complete bivariate MDA model; (f) MDA model with double-loop interaction of two signals plus an exogenous input.
interacting signals, as represented in Fig. 1. So, the MAR network is described by an with an upper block-triangular structure with for all and
The canonical form of the joint process is fixed by imposing that the white noises are uncorrelated also at zero-lag, in accordance with a complete uncorrelation of the model inputs. So the variance matrix of is diagonal
(8)
(9)
where is the matrix that represents the MAR network of the interacting signals and is matrix. The block (with in fact, exogenous signals are described by the single-channel AR processes fed directly by the corresponding white noises The corresponding is accordingly omitted; so that and (in fact an AR block would be directly in series to ). As shown in Fig. 1, the general structure obtained after the definition of exogenous signals contains specific loops which generate the model poles. This characteristics will be used in Section II-C to classify the rhythm sources of the model and to quantify their contribution to each signal by the relevant spectral components. B. Model Identification and Validation Here, the joint process description of the MV model is used to define the canonical form imposed and introduce the consequent identification procedure. It is also shown that the validation of a specific model structure must satisfy the corresponding joint process description of all the signals whether they are supposed exogenous or not. In other words, identification should result into white and uncorrelated prediction errors; otherwise, some of the a priori assumptions about the existing interactions should be rejected.
where is the variance of In order to fulfill this constraint, the presence of some immediate effects, , is generally required, because zero-lag cross correlations among the signals have to be explained by the modeled interactions rather than by correlations among the inputs. Loops of immediate effects are not allowed; therefore, matrix (which contains the , with null diagonal) is triangular. Consequently, and are triangular as well, with diagonal elements equal to one. Which immediate effects have to be considered depends on physical considerations; e.g., if the th sample is measured in real time before , then is estimated, while is set to zero. A parametric identification of the joint process requires the computation of an estimate of by minimizing a scalar function of the prediction error variance matrix (e.g., its determinant, with the constraints imposed by the chosen canonical form. The prediction-error vector appears in the prediction form of the model (10) where it indicates the difference between and its one-step prediction After a correct identification of represents an estimate of the input white noises Identifiability requires that all branches of the closed loops are enlightened by the inputs, i.e., that all are persistently exciting [9]. The diagonality of permits us to approach the estimate by means of direct methods and minimize the single prediction
BASELLI et al.: MULTIVARIATE SPECTRAL DECOMPOSITION
1095
complex
error variances, , on the diagonal by means of single-output prediction models. In other words, the prediction of from the other signals is optimized separately, as if the loops were virtually opened. In MDA models this leads to the identification of single-output DA models (ARXAR) with inputs
variable moves on the unit circle (with ), from phase 0 to In order to further decompose each partial spectrum into components corresponding to its peaks, it is necessary to classify the poles of , which, according to (6), are the roots of the determinant of
(11)
(14)
(12)
As clearly shown by Fig. 1, plus the roots of each the poles of the exogenous signals are separated from those borne by
The minimization of , through a generalized least squares , of algorithm, yields directly a parametric estimate of (with ), and of [20], [21]. From (11) it is clear that a random input is actually derived as equation residual, i.e., by difference from the considered interactions. If is absent, the algorithm is simplified to a least squares one. In particular, each exogenous signal requires to identify only its AR parameters The prediction errors obtained with the estimated parameters have to be tested in order to verify the identification hypothesis: 1) whiteness of each ; 2) uncorrelation between each pair ; 3) uncorrelation at zero-lag or diagonality of A residual autocorrelation or crosscorrelation of the prediction errors indicates either an insufficient order or an inappropriate model structure (e.g., an incorrect assumption of exogenousness or the exclusion of interactions that actually exist); a nondiagonal indicates an incorrect choice of the immediate effects. The order of the model corresponds to the order of its FIR filters, which can be either unique (order as above) or change from filter to filter (the above formula can be readily extended to this case). Usually, the model order is determined after successive increments until the hypotheses are fulfilled and a minimum of Akaike’s figure of merit is reached [22]. In conclusion, the MDA model layout is set by deciding: 1) which causal effects, are present (hence, the eventual presence of exogenous inputs); 2) which immediate effects, are allowed; 3) which self loops, , are hypothesized; 4) which residual inputs, are colored (presence of and the relevant poles). Identification is performed on single-output ARXAR models (AR for exogenous signals), but identification hypotheses and model adequacy are verified on the joint process of prediction errors.
(15) In addition, the roots of characterizing an exogenous may appear only in the transfer functions starting from the corresponding Therefore, the spectrum of an interacting signal is decomposed into the sum of partial spectra due to stochastic (nonmeasured) inputs and due to exogenous (measured) ones shaped by rational transfer functions with form (16) and
(17) respectively. In fact, each white input is first colored by the poles of the corresponding AR filter, either in (16) or in (17). These poles are specifically related to the considered input and, therefore, are named -poles for stochastic inputs or -poles for exogenous inputs. Next, the input is filtered by the MAR network of interacting signals which generates spectral peaks corresponding to the poles borne by ; these poles are named closed-loop poles (cl-poles) because they are related to loops of causal interactions in accordance with Mason’s formula (closed-loops)
C. Multivariate Spectral Decomposition
(pairs-of-separate-closed-loops)
If the identification hypotheses are satisfied, each interacting signal is described by the MDA model as the sum of uncorrelated contributions, each driven by a white noise Uncorrelation permits the decomposition of the spectrum of into a sum of partial spectra due to the different (i.e., the spectra of the
(13) is the sampling period expressed in seconds; the where frequency scale ranges from 0 to 1/2 Hz, while the
(18) are related to the interferences between Polynomials different pathways and loops and do not generate any component but contribute to give different weights to the cl-pole components in the different partial spectra. In conclusion, classes of poles are defined which characterize different oscillation sources even at similar frequencies: cl-pole components due to the modeled interactions; -pole components which quantify rhythmogenetic mechanisms external to the considered interactions and carried by the colored stochastic inputs; -pole components of the exogenous inputs that are directly measured.
1096
IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 44, NO. 11, NOVEMBER 1997
The partial spectra have form similar to that of a singlechannel autoregressive moving average (ARMA) process (a white noise filtered by a rational transfer function) and are accordingly decomposed into components relevant to each pole (or pair of conjugate poles), thus, obtaining their central frequency and power [7] (see Appendix for details). So each partial spectrum is divided into cl-pole components and either -pole or -pole components (when is exogenous). The global size of a cl-pole component is evaluated by summing again its contributions extracted from the different partial spectra and is compared with the size of the other types of components. In a finer detail, the partial contributions by which the cl-pole component is formed indicate the extent to which the closed-loop resonance is excited by each input. Central frequencies and either absolute or percent powers of components can be represented as bargraphs on the frequency axis (Fig. 6). Each bar can be labeled according to the pole type and colored according to the input to which it is due (i.e., the partial spectrum from which it is derived). Consequently, -pole (random inputs) and -pole (exogenous signals) components have the color of the corresponding partial spectrum, ; while, cl-pole components (derived from all partial spectra) result in stacked bars in which the contribution of each input is distinguished by a segment of different color. III. APPLICATION EXAMPLES Given signals experimentally recorded or monitored, the general MDA model of Fig. 1 allows all possible interactions and all possible colored stochastic inputs. This quickly increases the number of loops, i.e., of rhythm sources. A priori physiological knowledge about the system constraints, or working hypothesis, or even the skills of the researcher can be easily included in the model application by omitting some interaction paths, thus, greatly simplifying the model and tailoring it to the specific biological problem. Some signals can be considered exogenous by omitting all directed to them; this can be done when it is clear that these signals affect the other ones but are not part of control loops; e.g., respiration in respect to cardiovascular signals, an electrooculographic or electromyographic noise in respect to electroencephalographic signals, an external stimulus in respect to neurographic or myographic signals, the large heart period variability in respect to the small QT duration variability in the electrocardiogram, etc. Some other interactions can also be excluded by omitting the corresponding Self loops can be excluded by setting to zero the corresponding ; this is correct when it can be hypothesized that the control loops affecting that signal also involve some other recorded signal and, thus, are explicited by other loops; e.g., heart period variability when arterial pressure variability is considered, efferent activity in neural loops when the corresponding afferent activity is recorded, muscular activity when the antagonist activities are all included in the model, etc. Other a priori features can be included by hypothesizing that some stochastic inputs are not colored (white), by limiting the order of specific blocks, and by including delays in the proper block. In order to show how the MDA class can be easily adapted to the particular application, simple structures will be discussed
in Section III-A; the layouts of the specific models are shown in Fig. 2 with graphical skills intended to emphasize the loops, i.e., rhythm sources. In Section III-B, some of them will be utilized to analyze the interactions among cardiovascular variability signals, in order to emphasize the different results obtained in relation to the model complexity and specificity. A. General Examples Note the particular graphic arrangement adopted in Fig. 2 to depict the model layouts. The represented blocks are exclusively FIR filters and never rational transfer functions. Therefore, each block corresponds to a specific regression of the identification algorithm; autoregressions of measured signals are indicated by blocks, regressions of a signal over another one by , and autoregressions of stochastic In addition, the absence of poles “inside” inputs by the blocks permits us to emphasize the loops where the model poles are originated, and, consequently, the different rhythm sources. Fig. 2(a) shows that the single-signal case reduces to a single-channel AR model; the way in which it is drawn emphasizes the feedback performed by which generates the poles and the relevant spectral components of the signal. Separate AR spectral analyses performed on several signals can detect the existing rhythm, but do not provide any classification in terms of interaction. Fig. 2(b) shows the case of two signals, one of which is exogenous ; The relevant model is a single-output ARXAR. In this case, three classes of poles are distinguished in relation to: 1) the output feedback (cl-poles), 2) the spectral content of the exogenous signal ( -poles), and 3) the spectral content of the residual ( -poles). When two signals express the activity of two opponent mechanisms which interact in a closed loop, the model of Fig. 2(c) can be applied. The analysis tries to disentangle the contribution of resonances created by the signal interaction (clpoles) from external rhythms impinging either on ( -poles) or ( -poles). Fig. 2(d) shows the case of a bivariate MAR model. All the blocks of are present, while the random inputs are considered as white. As a consequence, only a class of cl-poles related to the triple-loop (tl-poles) of the model is obtained. Often a MAR model is used for a black-box analysis without relating its blocks to particular interactions. In the model of Fig. 2(e), all possible rhythm sources and interactions permitted by an MDA model with are included. The model of Fig. 2(f) considers the closed-loop interaction of with and also a self loop on ; this can be the case of a variable which is controlled by the action of the other signal , but also by different feedback mechanisms which cannot be measured and are lumped into the self loop. The resonances of the double-loop structure are expressed by dl-poles. The exogenous signal bears -poles, while the residual inputs are colored and bear -poles and -poles, respectively. B. Cardiovascular Example In the following, analysis examples are given concerning cardiovascular variability. Namely, three beat-by-beat series
BASELLI et al.: MULTIVARIATE SPECTRAL DECOMPOSITION
(a)
1097
(b) (a)
(c)
(d)
Fig. 3. Beat-to-beat variability series of (a) the systolic arterial pressure s and (c) heart period duration t of a conscious dog during BCO. The monovariate s and t AR spectra [(b) Ss (f ) and (d) St (f ) (solid lines), respectively] are shown with their relevant LF and HF spectral components (dashed lines).
are considered: systolic arterial pressure variability measured inside the th cardiac cycle , duration of the th heart period variability measured from the R-R interval , respiratory signal sampled at the beginning of the th cardiac cycle ; mean values are subtracted from all series; with a number of samples typically 300. The measurement order allows “immediate” effects (i.e., effects within the th cardiac cycle) from to and from to both and A single data set recorded during bilateral carotid occlusion (BCO) in a conscious instrumented dog will be analyzed first using separate monovariate AR models [Fig. 2(a)], next with a bivariate MAR model [Fig. 2(d)] applied on and , and finally by means of a specific MDA model with double-loop structure of and interaction and as exogenous input [Fig. 2(f)]. BCO was considered because this maneuver is particularly critical in relation to the interpretation of the physiological mechanisms enhancing the LF activity [24]. Indeed, due to the complexity of this maneuver’s ability to induce a sympathetic excitatory response by depulsating the baroreceptors at the carotid sinus level, but not unloading those at the aortic sinus level, several different control mechanisms could be contemporaneously active in producing LF oscillations. This data record is part of an experimental protocol carried out on conscious dogs in which various physiological maneuvers were performed in order to enhance sympathetic activity and, hence, LF oscillations in different ways: nitroglycerin infusion, coronary artery occlusion, and BCO [23], [24]. Series and and the respective results of monovariate AR spectral decomposition are shown in Fig. 3(a) and (c), respectively, relevant to a case of BCO. The monovariate spectra of Fig. 3(b) and (d) (solid lines) displays a modest LF component (dotted lines) at 0.1 Hz (virtually absent in control condition; not shown). These spectra also display a large fraction (compared to control) of power at lower frequencies,
(b)
(c) Fig. 4. MV spectral decomposition of series s based on the application of the bivariate AR model of Fig. 2(d) to the series s = y1 and t = y2 of Fig. 3. The total spectrum Ss (f ) (a) is decomposed into the sum of the partial spectra Ss=ws (f ) (b) and Ss=wt (f ) (c) driven at s and t level, respectively. A single type of components relevant to tl-poles is obtained in both partial spectra.
but the monovariate analysis was not able to resolve any specific rhythm in this band. A decreased (compared to control condition) yet marked effect of respiration on variability is documented by two large HF components (at close frequencies, Fig. 3(b) and (d); dotted lines), due to a small drift in respiratory frequency; respiration (not shown) also displays the same HF peaks. In Fig. 4, a bivariate analysis with the MAR model of Fig. 2(d) is performed on the same and data of Fig. 3. The bivariate analysis considers both signals at the same time; accordingly the total autospectra are slightly different and are based on more information. The total spectrum of computed by the model, [Fig. 4(a)], is decomposed into two partial spectra: [Fig. 4(b)], which describes the amount of variability driven at level (arterial pressure level), and [Fig. 4(c)], driven at level (sinus node level). A similar decomposition of into and (not shown) was obtained relevant to This model can be useful when characteristics such as gains and delays of and interactions are addressed; e.g., the ratio of over assesses the amount of variability driven at level and reflected on ; i.e., the gain of cardiac baroreflex mecha-
1098
IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 44, NO. 11, NOVEMBER 1997
(a)
(b)
(c)
(d)
=
=
Fig. 5. MV spectral decomposition of series s based on the application of the MDA model of Fig. 2(f) to series s y1 ; t y2 (same series of Fig. 3) and y3 (not shown). Ss (f ) is decomposed into the partial spectra driven at s and t level and by r : Ss=ws (f ), Ss=wt (f ) and Ss=wr (f ), respectively. All partial spectra contain components relevant to double-loop interactions (dl-poles) and each one contains the components relevant to the poles of the respective inputs: us ; ut ; and r: The different colors of the partial spectra (white, gray, and black, respectively) are used in the bargraphs of Fig. 6.
r
=
nisms described by the transfer function The further decomposition of the partial spectra into spectral components reveals by which amount each rhythm is driven at or at level, but does not provide any rhythm source classification, as the model considers only the triple-loop poles (tl-poles) excited by two white noises. The small component near 0.1 Hz is described as being excited almost entirely by The large amount of power at lower frequencies is not resolved by the model. Importantly, the lack of as exogenous signal compels the model to describe the HF components by means of fictional signal interactions; the prevalence of HF in indicates a strong direct modulation of on (a situation different from basal condition in the dog [18]). Therefore, the bivariate AR model structure allows us, with respect to the monovariate analysis, to evaluate the amount of power exchanged between the and series and to quantify the transfer function describing the - interactions. However, a rhythm classification cannot be carried out, because all possible rhythm sources are lumped together in the triple-loop structure. A more detailed analysis of rhythm sources and of respiration effects is performed through the model of Fig. 2(f). The double-loop network (dl-poles) represents the interaction between and elicited by cardiac baroreflex , by mechanical effects of heart rate variability on arterial pressure , and all other pressure regulation feedbacks The spectral content of and ( -poles and -poles, respectively) describe rhythms external to the considered interactions; the former ones can be related to rhythmic vasomotor activity while the latter to rhythmic neural activity modulating the sinus node. All these mechanisms were shown to be involved under different stimuli in the generation of LF waves [18], [19]. Respiration is described as a colored ( -poles) exogenous input impinging on both vari-
abilities, and is mainly involved in the description of the HF components. The model considers simultaneously the and series in order to extract the parameters of the various blocks and, next, to compute the relevant MV spectral decomposition. Fig. 5 displays the decomposition of signal obtained after the analysis of the same BCO record; with this model, the total spectrum [Fig. 5(a)] is decomposed into the sum of three partial spectra , and [Fig. 5(b)–(d), respectively]. In all the three partial spectra a LF dl-pole component appears, indicating that the closedloop interactions are involved in the generation of the LF rhythm by means of a resonance and frequency-selective amplification. This component appears at a frequency of 0.05 Hz and corresponds to waves lasting more than 30 beats visible in and even better in series [Fig. 3(a) and (c)]; so, this model, thanks to its specific structure, was able to clearly distinguish an oscillation even at this lower frequency. All the inputs do excite the dl-pole resonance, with a prevalence of Also, a small -component is detected at 0.12 Hz. The fraction of variability induced by respiration is described by [Fig. 5(d)], in which the breathing rhythm is represented by the HF -pole components; the LF dl-pole components in indicates that also slow respiratory modulations can excite the closed-loop resonance [18], [19]. The pole diagram of Fig. 6(a) displays the poles and indicates their frequency (phase in the -plane), effectiveness (closeness to the unit circle) and, also, their classification provided by the MV approach; this case evidences a dl-pole and a -pole at LF and two -poles at HF, which correspond to the components described above. The bargraphs of Fig. 6(b) and (c) indicate the percent powers and central frequencies of the spectral components resulting from the double-loop model identification; the former [Fig. 6(b)] bargraph summarizes the
BASELLI et al.: MULTIVARIATE SPECTRAL DECOMPOSITION
1099
(a)
(b)
(c)
Fig. 6. Summary of the spectral decomposition of s shown in Fig. 5 and of the results relevant to t simultaneously obtained. (a) Pole diagram with the identified dl-poles, ut -poles (at LF), us -poles (absent), and r -poles (at HF). Bargraphs representing the relative power and central frequency of the components of (b) s and (c) t. The partial spectrum considered and the relevant driving input (ws ; wt ; and wr ) are distinguished by the color of the bar (white, gray, and black, respectively); dl-pole components, which are present in all the three partial spectra, result in stacked bars and the proportions of the three colors indicate the relative importance of the inputs in exciting the closed-loop resonances.
decomposition of spectrum shown in Fig. 5, while the latter [Fig. 6(c)] is relevant to and is derived in a similar way. Considering Fig. 6(b), the components relevant to the dl-pole result into a stacked bar at LF; the gray part, relevant to the effect of , is prevalent over those due to (black) and to (white). This indicates that the LF rhythm is synchronized by the closed-loop regulation of , but is mainly driven by broadband modulations of The small -pole component is gray as excited only by ; -pole components are black; while a -pole component (not present) would be white. The decomposition bargraph of [Fig. 6(c)] is similar with an even more marked effect of in exciting the LF dl-pole. The LF enhancement observed in the series and, less remarkably on the series [24], is explained by the doubleloop model as a result of the resonance of the - double-loop interactions (dl-pole). However, the presence of an -pole suggests that other control mechanisms, external to the double loop and affecting primarily the heart period, are present. They appear to be able, first, to modulate the heart period duration and, next, to cause systolic arterial pressure changes as pointed out by the large LF percentage of the variability driven by changes [Fig. 6(b), gray part of the stacked bar]. The importance of the - closed-loop resonance at a frequency slower than the usual 0.1 Hz and the presence of different mechanisms external to the closed-loop interaction have been confirmed [19] on a small animal group (four dogs). It is remarkable that in the presence of a decreased gain of baroreceptive response (correctly identified by the model [19]), an increased pressure regulation resonance is found. Considering the usual simple control model, it can be hypothesized that the impairment of baroreceptive mechanisms can cause an increase in feedback delays, thus, decreasing the control stability with a tendency
to larger and slower oscillations. Conversely, considering the simultaneous presence in the neural cardiovascular control system of both inhibitory (negative feedback) and excitatory (positive feedback) mechanisms, the data can be interpreted as a consequence of a prevailing excitatory sympathetic circuit (operating at a slower frequency). Accordingly, the prevalent drive to this slow oscillation exerted through the heart period can be hypothesized to be an index of the concurrent activation of central cardiovascular control mechanisms induced by the disappearance of such an important inhibitory input as that represented by the afferences from the carotid sinuses. The BCO stimulus causes a large LF increase (compared to control) as expected [24]. However, the shift toward lower frequency renders this component difficult to resolve by the monovariate AR model (see Fig. 3) and also by the bivariate AR model (see Fig. 4). It is remarkable that the better-suited and more-complex double-loop model is able to improve the spectral resolution in addition to correctly identify the rhythm sources. IV. CONCLUSION The present method is proposed as a tool to be applied in fields of biological signal analysis, in which rhythms are classically studied by means of monovariate spectral analysis. It appears particularly suited where a limited number of signals is analyzed and sufficient a priori knowledge can be considered in the design of the MV model layout, so that the closed loops involved in the genesis of poles and of the relevant spectral components can be given specific physiological meanings. The flexibility of the model building procedure permits to explore different combinations of signals and different working hypotheses.
1100
IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 44, NO. 11, NOVEMBER 1997
The presentation of MV parametric analysis results provided by the proposed MV spectral decomposition is closely related to the familiar plots of autospectral analysis; in addition, the class of MDA models is suited to classify different rhythm sources. So, once the model structure is fixed, the complexity of MV parametric identification can be rendered transparent to final users even if more detailed information than in singlechannel analysis is provided.
The residuals for the computation of (A.3) for subsequent lags are derived from those at zero lag by multiplying by a factor (A.6) is divided into the sum of components Therefore, which are exponential series with base
APPENDIX ARMA SPECTRAL DECOMPOSITION This section shows a spectral decomposition method based on pole residuals. It was introduced for the analysis of ARMA spectra [7], but it is useful to decompose the addressed MV partial spectra. The general form of a monovariate ARMA spectrum (and also of the partial spectra introduced above) is (A.1) where ARMA
is the complex spectral function. transfer function
is the
(A.2) where the are the poles of the process with bears the zeros of the process which need not to be computed. is one for the canonical form of monovariate ARMA processes, but this is generally not true (apart on the main diagonal) for the transfer functions of a MV process; in any case, spectral decomposition is not affected by this feature. The characteristics of spectral decomposition can easily be shown by considering the inverse transform of , which yields the autocorrelation function (ACF)
(A.3) The line integral (A.3) can be obtained by summing components computed as residuals of the argument in respect to the poles inside the unit circle, i.e., the poles of ; e.g., the zero-lag sample is the sum of the residuals of , which are defined as
(A.4) and easily computed
(A.5)
,
(A.7) Finally, the desired spectral components are obtained by considering the bilateral -transform of the ACF components (A.8) Components relevant to real poles are real as well and produce spectral components centered either at 0 Hz or at Nyquist frequency, , depending on the pole sign. The combination of a pair of conjugate poles produces a real ACF component that is a damped sinusoid with frequency and damping time constant (i.e., damping is slower when the pole is closer to the unit circle). For the ACF symmetry, the relevant spectral component is real as well, , and has the desired characteristics of being bell-shaped and centered around the pole frequency, , with a peak growing higher and sharper when the pole approaches the unit circle. The skewness of the components is highly influenced by the zeros. The component power is for real poles, and for conjugate pole pairs; in fact the sum of the gives the process variance which represents its total power. The above result can be derived directly in -domain by Heaviside decomposition of in simple fractions relevant to all its poles; i.e., both the inside the unit circle and their reciprocals outside, for The simple , i.e., fractions are weighted by the relevant residuals of and , respectively. A spectral component relevant to a pole is derived from its simple fraction and that of its reciprocal (A.9) In fact, the former simple fraction in (A.9) is the third addendum in (A.8); while, by considering that it can be readily verified that the latter simple fraction in (A.9) corresponds to the first and the second addenda in (A.8). REFERENCES [1] A. Cohen, Biomedical Signal Processing. Boca Raton, FL: CRC, 1986. [2] R. Weitkunat, Digital Biosignal Processing. Amsterdam, the Netherlands: Elsevier, 1991. [3] D. O. Walter and W. R. Adey, “Analysis of brainwave generators as multiple statistical time series,” IEEE Trans. Biomed. Eng., vol. 12, pp. 8–13, 1965. [4] W. Gersch and J. Yonemoto, “Parametric time series models for multivariate EEG analysis,” Comput. Biol. Res., vol. 10, pp. 113–125, 1977.
BASELLI et al.: MULTIVARIATE SPECTRAL DECOMPOSITION
[5] S. M. Kay, Modern Spectral Analysis: Theory and Application. Englewood Cliffs, NJ: Prentice-Hall, 1988. [6] S. L. Marple, Digital Spectral Analysis with Applications. Englewood Cliffs, NJ: Prentice-Hall, 1987. [7] L. H. Zetterberg, “Estimation of parameters for a linear difference equation with application to EEG analysis,” Math. Biosci., pp. 227–275, 1979. [8] L. Ljung, “System identification,” Theory for the User. Englewood Cliffs, NJ: Prentice-Hall, 1987. [9] T. S¨oderstr¨om and P. Stoica, System identification. London, U.K.: Prentice-Hall Int, 1989. [10] P. J. Brockwell and R. A. Davis, Time series: Theory and methods. Berlin, Germany: Springer-Verlag, 1987. [11] P. E. Caines, Stochastic Linear Systems. New York: Wiley, 1987. [12] W. Gersch, “Causality and driving in electrophysiological signal analysis,” Math. Biosci., vol. 14, pp. 177–196, 1972. [13] J. P. Saul, R. D. Berger, M. H. Chen, and R. J. Cohen, “Transfer function analysis of autonomic regulation. II: Respiratory sinus arrhythmia,” Amer. J. Physiol., vol. 256, pp. H153–H161, 1989. [14] P. J. Franaszuck, K. J. Blinowska, and M. Kowalczyk, “The application of parametric multichannel spectral estimates in the study of electrical brain activity,” Biol. Cybern., vol. 51, pp. 239–247, 1985. [15] M. Morf, A. Vieira, D. T. L. Lee, and T. Kailath, “Recursive multichannel maximum entropy spectral estimates in the study of electrical brain activity,” Biol. Cybern., vol. 51, pp. 239–247, 1985. [16] J. C. Shaw, “Correlation and coherence analysis of the EEG: A selective tutorial review,” Int. J. Psycho-Physiol., vol. 1, pp. 255–266, 1984. [17] G. Baselli, S. Cerutti, S. Civardi, A. Malliani, and M. Pagani, “Cardiovascular variability signals: Toward the identification of a closed-loop model of the neural control mechanisms,” IEEE Trans. Biomed. Eng. vol. 35, pp. 1033–1046, 1988. [18] G. Baselli, S. Cerutti, F. Badilini, L. Biancardi, A. Porta, M. Pagani, F. Lombardi, O. Rimoldi, R. Furlan, and A. Malliani, “Model for the assessment of heart period and arterial pressure variability interactions and of respiratory influences,” Med., Biol. Eng., Comput., vol. 32, pp. 143–152, 1994. [19] G. Baselli, A. Porta, G. Ferrari, S. Cerutti, M. Pagani, and O. Rimoldi, “Multi-variate identification and spectral decomposition for the assessment of cardiovascular control,” in Computer Analysis of Blood Pressure and Heart Rate Signals, M. Di Rienzo, G. Mancia, G. Parati, A. Pedotti, A. Zanchetti, Eds. Amsterdam, the Netherlands: IOS Press, 1995. [20] T. S¨odestr¨om, “On the convergence properties of the generalized least squares identification method,” Automatica, vol. 10, pp. 617–626, 1974. [21] J. Gustavsson, L. Ljung, and T. S¨oderstr¨om, “Survey paper, Identification of processes in closed loop: Identifiability and accuracy aspects,” Automatica, vol. 13, pp. 59–75, 1977. [22] H. Akaike, “A new look at the statistical novel identification,” IEEE Trans. Automat. Contr., vol. 19, pp. 716–723, 1974. [23] O. Rimoldi, S. Pierini, A. Ferrari, S. Cerutti, M. Pagani, and A. Malliani, “Analysis of the short term oscillations of RR and arterial pressure in conscious dogs,” Amer. J. Physiol., vol. 258, pp. H967–H976, 1990. [24] O. Rimoldi, M. R. Pagani, S. Piazza, M. Pagani, and A. Malliani, “Restraining effects of captopril on sympathetic excitatory responses in dogs: A spectral analysis approach,” Amer. J. Physiol., vol. 267, pp. H1608–H1618, 1994.
Giuseppe Baselli received the Italian degree in electronic engineering from the Polytechnic University of Milano, Milan, Italy, cum laude, in 1983. He worked as a Research Fellow at the Department of Electronic Engineering of the same University. Since 1986 he is Research Fellow at the Department of Electronics for Automation, University of Brescia, Brescia, Italy, where he teaches Automatic Control and Biomedical Instrumentation. He also teaches a postgraduate course of the School of Endocrinology in the Medical Faculty, University of Brescia. His research activity is in the field of biomedical signal processing and physiological modeling. He deals mainly with the assessment of cardiovascular control mechanisms.
1101
Alberto Porta received the Italian degree in electronic engineering at the Polytechnic University of Milano, Milan, Italy, in 1989. Since 1996, he has been studying towards the Ph.D. degree in Biomedical Engineering at the same university He has been a Research Fellow on Automatic Control and System Theory at the University of Brescia, Department of Electronic Engineering. His primary interests include time series analysis, signal processing of nonlinear dynamics, and system identification and modeling applied to cardiovascular control mechanisms.
Ornella Rimoldi was born in 1954 and received the MD degree from the University of Milano, Milan, Italy, magna cum Laude. She is specialized in cardiology and clinical pharmacology. Since 1988 she is a Research Fellow with the Italian Research Council. In 1994 and 1995, she was a Visiting Lecturer at Harvard Medical School, Cambridge, MA. Since 1996, she has been an Honorary Registrar and Research Fellow with the Royal Postgraduate Medical School Hammersmith Hospital, London, U.K. She is involved in the study of neural cardiovascular control and computer application to blood pressure and heart rate variability analysis in experimental animals.
Massimo Pagani was born 1944 and graduated in medicine at the University of Milano, Milan, Italy, in 1969. From 1975–1977, he studied at Harvard Medical School, Cambridge, MA, under a Fogarty Fellowship, in the Department of Medicine. He is currently Professor of Medicine at the University of Milano. His main interest is the application of bioinformatics to the neural control and modeling, as well as the introduction of behavioral principles to clinical medicine. Dr. Pagani is Fellow of the American College of Cardiology and Fellow of the Circulation Council of the American Heart association, and member of many other professional societies.
Sergio Cerutti (M’81–SM’97) received the Laurea degree in electronic engineering from the Polytechnic University, Milano, Milan, Italy, in 1971. He spent over a year as a Visiting Professor at the Massachusetts Institute of Technology (MIT) and Harvard School of Public Health, Boston, MA. From 1982–1990, he was an Associate Professor of Biomedical Engineering at the Department of Biomedical Engineering of the same Polytechnic. From 1990–1994, he was Professor of Biomedical Engineering at the Department of Computer and System Sciences of the University of Rome “La Sapienza,” Italy. Since 1994, he has been a Professor of Biomedical Signal and Data Processing. His research activity is mainly dedicated to various aspects of biomedical signal processing and modeling related to the cardiovascular system and in the field of neurosciences. He is the author of various papers and books on these topics published in international scientific literature. Dr. Cerutti has been member of AdCom of IEEE-EMBS (1992–1996) as Representative of Region 8. He is Chairman of the Biomedical Engineering Group of Italian AEI (Association of Electrical Engineering) and is a member of IEEE-AEI, IFMBE-AIIMB, ESEM, IEC-CEI, and other international and national scientific associations.
