IV. Frequency domain
A. Spectral analysis
Based upon the above discussion, we can see that time domain methods are extremely useful, if and only if the periodicity of interest has a relatively high amplitude (i.e., represents a large proportion of the variance of the signal). The physiological response systems that psychophysiologists study seldom meet this criterion. We are, therefore, forced to identify other statistical methods which might be more amenable to our needs. An alternative approach to the study of periodicities is to apply frequency domain methods.
In contrast to the time domain techniques, frequency domain techniques are those based upon the spectral density function. Spectral analysis is a mathematical technique which decomposes a time series into constituent frequencies or periodicities. The amplitude or variance associated with each frequency component is known as the spectral density estimate. There is a mathematical relationship between the time domain autocorrelation procedures and spectral analysis. The spectral density function is the Fourier transform of the autocovariance (unstandardized autocorrelation) function. The Fourier transform is an algebraic method of decomposing any time series into a set of pure sine waves of different frequencies, with a particular amplitude and phase angle for each frequency. The algebraic sum of the sinusoidal components adjusted for phase angle will accurately reproduce the original time series. Moreover, the inverse Fourier transform of the spectral densities will reproduce the autocovariances.
In its strictest sense, however, the Fourier analysis of a time series is the decomposition of the series into the sum of sine and cosine terms. This decomposition is called the periodogram. The periodogram is the finite Fourier transform of the autocovariance function. Although the periodogram is relatively easy to calculate, it does not provide consistent estimators for the spectral densities. Thus, methods are required to smooth the periodogram to provide spectral density estimates which have less bias and variance. Spectral analysis includes methods of smoothing or windowing to enhance the description of stable estimators with stochastic processes. Recall that deterministic processes are totally predictable and that physiological processes contain probabilistic components and are, therefore, stochastic. This means that the parameters of the periodic process (i.e., amplitude, phase, and frequency) can not be totally predicted by the past history of the signal.
Figure 1c illustrates the power spectrum of the sine wave illustrated in Figure 1a. Note the dominant peak at 0.2 Hz in the spectrum of power densities. This means that a major source of variance in the time series may be characterized by a sine wave with a period of 5 seconds (i.e., 0.2 Hz). Note that the power spectrum attributes spectral densities to frequencies adjacent to the dominant frequency. If spectral analysis accurately decomposed the time series of a pure sine wave, why does the spectrum illustrated in Figure 1c have spectral density (i.e., variance) estimates at frequencies other than 0.2 Hz?
The answer is based upon an understanding of the application of spectral technology to stochastic processes. If we performed a Fourier transform on the deterministic sine wave plotted in Figure 1a, then the corresponding spectrum would distribute all the variance at 0.2 Hz. However, Figure 1c represents the results of a technology which has been developed to deal with stochastic processes. This analysis includes procedures to deal with the stochastic nature of the data. One of these procedures is called windowing.
Windowing is a technique for smoothing the spectral density function derived from the Fourier transform. The objective is to minimize the bias and variance of the spectral density estimates. Windowing may be implemented in either the time or frequency domain. In the frequency domain this may be accomplished by summing the weighted spectral density estimates on both sides of a specific frequency of interest. By stepping this weighted function across the spectrum (i.e., periodogram), abrupt shifts in the spectrum will be smoothed which may have represented unstable estimates. The number of frequencies included in the window is often determined by the investigator and the coefficients are a function of the specific window selected. A simple frequency domain window may be generated by calculating a moving average, two points to the left and two points to the right of a target frequency and stepping this through the spectrum. This unweighted moving average is called the Daniell window. The expected value of the new spectral density estimate becomes less biased as the duration of the time series increases. As more points are added to the moving average in the frequency domain, the variance decreases. Other frequency domain windows such as the Tukey-Hanning, Hamming, and Parzen use more complex weighing coefficients. Windowing may also occur in the time domain. For example, the Bartlett window is often used with the fast Fourier transform (FFT). This procedure consists of performing the FFT on segments of the data and summing the spectra to provide one smoothed spectrum. More detailed discussion of windows may be obtained in time series textbooks (see Chatfield, 1975; Got tman, 1981).
There are a variety of methods to calculate spectral density estimates. A common technique is called the fast Fourier transform (FFT). This technique greatly reduces the time required to perform a Fourier analysis on a computer. A limiting feature of most FFT algorithms is that they require data sets to be in powers of two (e.g., 128, 256, 512, 1024, etc). A discrete spectral analysis can be calculated with any number of data points by transforming the autocovariance function. The statistical properties of both the discrete Fourier transform and the FFT can be enhanced through frequency domain windowing. For a more mathematically based discussion of these issues see Bohrer and Porges (1982) or a recent time-series text such as Chatfield (1975) or Gottman (1981).
The value of spectral analysis becomes obvious when it is used to describe the parameters of periodicities in a complex wave. Spectral analysis provides a method to extract information about periodicities, even when they are embedded in a complex pattern and may represent only a small proportion of the total variance of the time series. Figure 5a illustrates a complex pattern composed of the sum of three sine waves of equal amplitude. The sinusoidal components of this wave form are illustrated in Figure 5b. The frequencies in this example have been selected to be similar to those found in the heart rate pattern and have been associated with thermoregulation (0.04 Hz), blood pressure (0.1 Hz), and respiratory influences (0.26 Hz). For our data analyses we have sampled this complex signal every 500 msec. Figure 5c illustrates the confusing autocorrelogram of this complex process. In our example a "lag" is equivalent to the 500 msec sampling interval and the three sine waves in the signal have periods of 25, 10, and approximately 4 seconds. If the autocorrelogram is of use in detecting periodicities in a complex signal it should exhibit periodicity every 50 lags for the 0.04 Hz sine wave, every 20 lags for the 0.1 Hz sine wave and approximately every 8 lags for the 0.26 Hz sine wave. Although the autocorrelogram exhibits periodicities, it is difficult to accurately identify underlying periodicities in these data. The corresponding spectral analysis illustrated in Figure 5d decomposes the variance of the process into constituent frequency components between 0.0 and 1.0 Hz. Note the easily identified frequencies and the similarity of the spectral density estimates of the three peaks reflecting the similar amplitudes of the component sine waves. The effect of frequency domain smoothing is also evident in the width of the peaks and the overlapping of the spectral density estimates for the 0.04 and 0.1 Hz signals. Recall that the frequency domain windows are employed to deal with the stochastic characteristics of the signal, even though in this example we used a deterministic simulation.
A more appealing example is provided by the spectral analysis of the data illustrated in Figure 6a. As illustrated in Figure 6b, these data are composed of the three frequencies used in Figure 5 with different amplitudes. The 0.04 Hz sine wave has an amplitude of 100, the 0.1 Hz sine wave has an amplitude of 50, and the 0.26 Hz sine wave has an amplitude of 25. Figure 6c illustrates the autocorrelogram which is dominated by a periodicity of approximately 50 lags reflecting the impact of the high amplitude 0.04 Hz sine wave which has a period of 25 seconds or 50 lags.
Figure 6d demonstrates the effectiveness of the spectral procedure in extracting estimates of the amplitude and frequency of the three frequency components even when one component represents a small proportion of the variance. Closer inspection illustrates that the power density estimates for each of the three frequencies reflect their relative variances. Recall that amplitude squared divided by two is equivalent to variance. Thus, amplitudes in the ratio of 4:2:1 would produce variances in the ratio of 16:4:1. Since power density estimates are a linear transformation of variance, the three peaks in the power density spectrum reflect this ratio.
We can now inspect the spectral decomposition of the respiration and heart period data presented in Figures 2c, 3c and 4c. Figure 2c illustrates the spectrum for the respiration signal. The spectrum is characterized by one dominant frequency at approximately 0.18 Hz (i.e., approximately 5.6 seconds per breath or 11.8 breaths per minute). Very little of the variance from the respiration recording is attributed to other frequencies. In contrast when we inspect the spectral densities for the heart period patterns a different picture emerges. Figure 3a illustrates the heart period pattern recorded simultaneously with the respiration data illustrated in Figure 2a. Inspection of Figure 3c identifies a prominent peak associated with the respiratory frequency, although most of the variance (i.e., spectral densities) of the heart period process has been attributed to frequencies slower than respiration. Figure 4c illustrates the spectrum for the heart period pattern displayed in Figure 4a. Note that most of the heart period variance is associated with extremely low frequencies, while the peak at the respiratory frequency is low.
As illustrated in Figure 8a, the heart period data and the respiration data illustrated in Figures 2a and 3a are superimposed to emphasize the covariation of the two processes. The impact of the respiratory system on heart period (i.e., respiratory sinus arrhythmia) is characterized generally by a decrease in heart period (i.e., heart rate acceleration) with inhalation and increase in heart period with exhalation. In our example, respiration was monitored by quantifying the changing chest circumference associated with breathing. Inspiration produces an increase in chest circumference and an increase in our measure of respiration amplitude. When the two processes are superimposed, the rhythmic respiratory process is clearly manifested as respiratory sinus arrhythmia in the heart period process. The heart period is systematically attenuated during inhalation (i.e., increases in chest circumference or respiration amplitude in our example) and increased during exhalation. In our example, the amplitude of respiratory sinus arrhythmia appears to be more variable than the relatively constant amplitude of the respiratory process.
Although the spectral decomposition theorem is quite straight forward, its application to physiological data is extremely complicated. The spectral decomposition function decomposes all time series into constituent sine waves, each sine wave with an amplitude proportional to the variance the sine wave accounts for in the process. With pure sine waves the spectral decomposition is not ambiguous and functions to perfection. When time series are the sum of pure sine waves whose periodicities are hard or impossible to decipher from a graphic representation of the time series or from the autocorrelogram, the Fourier transform perfectly identifies the frequency and amplitude of the periodic components. Although spectral analyses will decompose any time series into constituent sinusoidal components, physiological activity is a composite of aperiodic and periodic components. Moreover the periodic activity is not sinusoidal and is seldom stationary. Therefore, spectral decomposition will distribute variance at low frequencies for the trends and slow aperiodic activity and to high frequencies for the harmonics of the periodic processes and high frequency components of the aperiodic process. Thus, the peak frequencies identified in the spectral analysis of physiological data may not correspond to periodic physiological process.
The above paragraph describes a situation in which spectral analysis does not provide consistent spectral density estimates. When the data are not stationary, no inference can be made regarding the characteristics of the spectrum. If, however, the time series is stationary, it is possible to obtain consistent estimates of the processes from a single finite sample of data. This important point has been mathematically proven by the ergodic theorems, which show that, for most stationary processes, the sample moments (i.e., mean, variance, autocovariance, etc.) of an observed record converge to the corresponding population moments.
The above characteristics of physiological time series make the application of spectral technology to psychophysiological research difficult. Figures 2c, 3c, 4c, and 8b provide appropriate examples to discuss three major concerns about spectral analysis of heart period patterns.
First, when the greatest power density occurs at the lowest frequency of the spectrum the data are non-stationary and violate a primary assumption underlying spectral analysis. Spectral analysis always works. It always decomposes the variance of the process into sine wave components. Thus, regardless of potential violations of assumptions, spectral analysis will always provide descriptive information. However, the designation of peak in the spectrum as being statistically significant implies that the underlying assumptions have not been violated. The spectral decomposition of an aperiodic non-stationary process will distribute a large amount of variance in the lowest frequency band. To deal with this problem, many researchers detrend the original data with low order (e.g., linear) regression fits across the entire data set prior to analysis. This method functions well, if and only if, the data have a linear or low order trend. Although these trends often account for a significant proportion of the variance, the residual series are seldom stationary. In the section on Filters we will describe methods developed to remove the non-stationary components of time series. The only other sources for variance in the lowest frequency are: one, if the researcher failed to remove the mean level from the data prior to analysis; or two, if there is a real periodicity at a frequency slower than the slowest defined frequency band. If the latter is the case, there may still be problems in the interpretation of faster periodicities in the spectrum. Since physiological processes are not perfect sine waves, extremely slow periodic processes might contribute variances to faster frequencies in the spectrum.
Second, the spectra of physiological processes do not always reflect the multiple periodic processes theoretically assumed to be present. For example, in the literature it has repeatedly been stated that there are three periodicities in the heart rate spectrum. These rhythms have been theoretically associated with respiratory influences (i.e., respiratory sinus arrhythmia) at approximately 0.12 to 0.40 Hz, blood pressure feedback (i.e., Traube-Hering-Mayer wave) at approximately 0.08 to 0.10 Hz, and thermoregulation at approximately 0.04 Hz (see Kitney & Rompelman, 1980). It is hard to identify any peaks in Figure 4c. In Figure 3c the 0.04 Hz and 0.1 Hz periodicities are not represented. However, with some imagination one can identify a periodicity of approximately 10 to 12 seconds and another periodicity of approximately 20 to 30 seconds with careful visual inspection of Figure 3a and Figure 4a. The ability to study these rhythms is greatly influenced by the aperiodic characteristics of the heart period pattern. Appropriate filtering might enhance the possibility of confirming whether these signals exist and the ability to study their behavior during psychophysiological experimentation. However, there are instances when these rhythms do not exist. These rhythms which reflect neural modulation of the heart may be greatly attenuated or blocked due to individual differences in neural control, pharmacological treatments, disease states or physiological stress (Porges, 1986; Porges, in press; Porges, McCabe, & Yongue, 1982).
Third, the dominant frequency in one physiological response spectrum is not necessarily the dominant frequency in another physiological response spectrum, even when both variables have common mechanisms. In our example in Figure 8b, the peak of the respiration spectrum is at a slightly faster frequency than the peak of the heart period spectrum. Yet, our observations of the two processes in Figure 8a, would lead us to believe that the rhythmicity of the breathing and the manifested respiratory sinus arrhythmia would be at the same frequency. The peaks may differ because the periodicities in the heart period pattern are less sinusoidal than the periodicities in the respiratory data. The heart period pattern includes more low frequency and aperiodic influences. It is possible that the "respiratory" peak in the heart period spectrum is, in part, influenced by variances associated with low frequency and aperiodic influences and is, thus, not reflecting the frequency and amplitude associated with an accurate measure of respiratory sinus arrhythmia. To assess whether this is truly happening, in a subsequent section we will attempt to filter from the heart period pattern all aperiodic influences and periodicities slower than breathing. If our speculation is correct, the appropriately filtered data will have a peak frequency more similar to the peak of the respiration series.
B. Bivariate time series
Often there is interest in evaluating the covariation of two time series. Psychophysiological constructs such as response fractionation (Lacey, 1967) and cardiac-somatic coupling (Obrist, Webb, Sutterer, & Howard, 1970) have been defined in terms of the statistical relationship among physiological response systems.. Similarly, arousal theory has been tested by assessing multiple variables. However, it is possible that the descriptive statistics of mean and variance used to define these constructs are not sensitive to all the underlying organizational characteristics of the nervous system. Specifically, mean and variance statistics are insensitive to rhythmic co-occurrence (i.e., coherence). Thus, the conclusions of earlier research (e.g., Lacey, 1967) may have been a function of an insensitive statistical methodology rather than of the underlying principles of neural and behavioral organization.
The cross-correlation function is a bivariate time domain method which is useful in describing co-variation of two response systems when each is characterized by the same predominant periodicity. In Figure 7a the cross-correlogram is illustrated for the simultaneously recorded respiration and heart period illustrated in Figure 8a. Note that the periodic nature of the respiratory rhythm is reflected in the cross- correlogram. The magnitude of these correlations oscillate and dampen from approximately +/- .50 to approximately +/- .25 at lags representing multiples of a time displacement of approximately 12 lags or 6 seconds.
Since physiological time series are seldom characterized by a single periodicity, the bivariate time domain methods like the cross-correlation suffer from the same weaknesses attributed to the autocorrelation. A more effective method of describing the co-variation of two physiological processes is the frequency domain technique of cross-spectral analysis.
Refer back to our example of simultaneously recorded respiration and heart period illustrated in Figure 8a. Note that there appears to be a synchronous pattern of decreases in heart period associated with inhalation and increases in heart period associated with exhalation. In Figure 8b the spectra are superimposed and there is a peak in heart period spectrum approximately at the peak of the respiration spectrum. Figure 8b does not provide the relevant information regarding covariation of the two processes. The figure merely describes their similarity in the spectral decomposition. It is possible to have similar spectra for both processes during states in which the processes are totally unrelated. Cross- spectral analysis is the appropriate frequency domain method of assessing the covariation of respiration and heart period at the breathing frequencies.
Cross-spectral analysis generates a coherence function which is a measure of the best linear association of each observed rhythm in one variable on the same rhythm in the second variable. Coherence is the square of the correlation between the two processes at a specific frequency. The coherence at any specific frequency is the square of the cross-spectral density divided by the product of the spectral densities of each series at the specific frequency. Note the similarity of this relationship with the calculation of the squared correlation coefficient. The cross-spectral density parallels the squared cross-products and the spectral densities parallel the variances. Conceptually, the coherence may be thought of as a time-series analogue of eta or the proportion of variance accounted for by the influence of one series on the other at each specific frequency. For a more detailed discussion of coherence and a relevant modification of the procedure for psychophysiological data see Porges et al. (1980).
Figure 8c illustrates the coherence spectrum for the heart period and respiration data illustrated in Figure 8a. Coherence scores are generated for all frequencies in the spectrum. Coherence, similar to the squared correlation coefficient, varies from O to 1.0. Note that the coherence is not constant across all frequencies and there is a peak in the coherence spectrum in the general area that the dominant respiratory frequencies occur. Careful inspection of Figure 8c suggests that the coherence spectrum is exhibiting peaks at the dominant frequency of breathing and at integer multiples of this frequency (i.e., harmonics). These harmonics peaks do not represent important or interpretable data. Rather they are caused by the fact that both respiration and heart period processes are not perfectly sinusoidal and the sine wave fit to the periodicity leaves variances at higher frequencies. These small variances (note in Figure 8b that the spectral densities are extremely small at frequencies faster than respiration) are associated with the respiratory process and its manifestation on heart period and therefore should exhibit a high degree of covariation. The relatively high coherence at the harmonic frequencies confirms this point.
When the coherence is high, it is possible to interpret the phase relationship between the two series. Figure 9a illustrates the phase spectrum for the cross-spectral analyses conducted on the respiration and heart period data illustrated in Figure 8a. The frequency in the phase spectrum associated with the high coherence between respiration and heart period is identified by an arrow. Note that the phase varies between +/- Pi (i.e., with Pi equal to 180 degrees or one-half cycle). The phase at the frequency of peak coherence suggests that the heart period oscillations are approximately one-half cycle delayed from the respiration. However, inspection of Figure 8a demonstrates that the oscillations are virtually in phase. This contradiction is caused by the scaling of variables. When respiration amplitude increases during inspiration, heart period decreases. If we were to plot respiration amplitude and heart period on similar axes (i.e., larger numbers higher on the ordinate), the two series would appear to be temporally displaced by approximately one-half cycle and thus, consistent with the analysis.
Since physiological rhythms are never a perfect sine wave and occur over a band of frequencies, the calculation of a summary statistic that describes the proportion of shared variance of the two systems would be useful. In previous publications (see Porges, et al.,1980; Bohrer & Porges, 1982) we have described a method of modifying the coherence function which would enable one number to describe the general relationship or coupling between two time series. In our example, we have selected a frequency band between 0.12 and 0.4 Hz to calculate the summary statistic for the coherence spectrum, the weighted coherence (Cw).
If the spectral densities were equally distributed, an unweighted mean coherence (for all frequencies) would accurately described the relationship. For any other situation, it is necessary to calculate a weighted coherence, which provides an exact measurement of the proportion of the variance of one series that is shared between the two processes.
In our example, the Cw between heart period and respiration has been defined as the proportion of total variance of the heart period process that is shared with the respiration process within the frequency band defined between 0.12 Hz and 0.4 Hz. This frequency band has been selected because it is most representative of the spontaneous respiratory rhythms in human adults. During various experimental conditions it might be functional to shift this band to accurately describe the respiratory behavior of the subjects. In our example, we have weighted each of the coherence estimates within the 0.12 Hz to 0.40 Hz frequency band by the spectral density for the heart period process at that frequency. This procedure assumes that respiration is driving the heart period process. If we assumed that heart period was driving respiration, we would weight the same coherence estimates by the spectral densities for the respiration process. In our example Cw is .76.
C. Vulnerabilities of spectral analysis
Spectral analysis decomposes the variance of the time series into constituent sine waves of various frequencies and amplitude. For spectral analysis to appropriately function it is necessary that the data are stationary and that the data are not aliased. If the data do possess these two properties, the spectral analysis may not be meaningful.
Stationarity implies that the expected value of the mean and variance are independent of the segment sampled and that the autocovariance function is solely a function of lag. Stationarity assumes that statistical properties of the time series do not change over time. Implicit in the property of stationarity is the assumption that time of initiation of an experimental manipulation is irrelevant. Thus, the application of spectral analysis which assumes stationarity would be inappropriate to many psychophysiological manipulations in which the experimental manipulation violates the stationarity property by creating a change in mean heart rate or heart rate variance.
Aliasing is a problem associated with selecting an inappropriate sampling rate. Spectral analysis decomposes the signal into component variances associated with frequencies from zero to one-half the sampling rate. Thus, if respiration were sampled every 500 msec (i.e., 2.0 Hz), spectral analysis as illustrated in Figure 2c would decompose this signal into frequencies from zero to 1.0 Hz. The fastest frequency we can detect in any spectrum is called the Nyquist frequency. The Nyquist frequency for any data set, assuming equal sampling intervals, is the frequency of a sinusoid whose period is twice the time interval between successive cycles. In our examples in Figures 1c, 2c, 3c, and 4c, the data were sampled at 2.0 Hz (i.e., sampling interval = 500 msec) and the Nyquist frequency is 1.0 Hz (i.e., period = 1000 msec). If there are oscillations faster than 1.0 Hz, then sampling at 2.0 Hz is too slow and the spectral analysis of the time series would be insensitive to these faster periodicities. The consequences of sampling too slow is two-fold. First, it would be impossible to visually identify rhythms whose period is less than the distance between time sampled observations. Second, the variance of the faster process is not totally lost or "filtered" from the data, but is "folded back" on the slower frequencies in the spectrum. Thus, it is possible for high frequencies to be aliased as low frequencies.
Figure 10a illustrates a sine wave and Figure 10b illustrates the time series derived from the sine wave when the data are sampled too slowly. Assume that the Figure 10a represents a 3 Hz sine wave and Figure 10b illustrates data sampled from the 10a sine wave at the rate of 5 Hz. Note the difference in number of oscillations between the two figures: the original sine wave has three oscillations and the sampled data exhibit two oscillations. In this example, the Nyquist frequency is 2.5 Hz and the variance associated with the 3 Hz sine wave would be aliased as frequencies below 2.5 Hz.
If we know the frequencies embedded in the physiological signal, we can avoid the problems associated with aliasing by increasing the sampling rate. Unfortunately, in many electrophysiological response systems we are not aware of the entire range of periodic activity. In fact, one application of spectral technology has been to describe this range. Without an appropriate a priori physiological knowledge base, this approach might provide inexplicable epiphenomenon. In addition, inappropriate electrode application may, in itself, produce high frequency artifact in the signal which may be aliased as slower frequencies.
