Detecting the Harmonics of Oscillations with Time-variable Frequencies

A method is introduced for the spectral analysis of complex noisy signals containing several frequency components. It enables components that are independent to be distinguished from the harmonics of nonsinusoidal oscillatory processes of lower frequency. The method is based on mutual information and surrogate testing combined with the wavelet transform, and it is applicable to relatively short time series containing frequencies that are time variable. Where the fundamental frequency and harmonics of a process can be identified, the characteristic shape of the corresponding oscillation can be determined, enabling adaptive filtering to remove other components and nonoscillatory noise from the signal. Thus the total bandwidth of the signal can be correctly partitioned and the power associated with each component then can be quantified more accurately. The method is first demonstrated on numerical examples. It is then used to identify the higher harmonics of oscillations in human skin blood flow, both spontaneous and associated with periodic iontophoresis of a vasodilatory agent. The method should be equally relevant to all situations where signals of comparable complexity are encountered, including applications in astrophysics, engineering, and electrical circuits, as well as in other areas of physiology and biology.


I. INTRODUCTION
The identification of periodic processes and their harmonics in noisy signals is an important data processing task in a variety of fields, and a number of numerical approaches to the problem have been developed in recent years, particularly in astronomy [1][2][3][4] and electrical engineering [5][6][7][8][9].One reason for continued interest in the problem is that useful information about the dynamics of a nonlinear system often can be obtained from a study of the harmonics in its response to an oscillatory input [10][11][12][13][14][15].It is an approach that has proven very useful in the analysis of physiological time series, where the underlying system is highly complex, nonlinear, and resistant to modeling approaches.
Physiological oscillatory processes arise in, e.g., electrocardiograms (ECGs), electroencephalograms (EEGs), blood pressure (Traube-Herring-Mayer waves), and intracranial pressure monitoring (Lund waves).The profile of the pressure wave in the circulatory system produced by each heartbeat provides useful information about the state of the heart and the cardiovascular system, as does the form of the ECG in the case of atrial arrhythmias [16], atrial fibrillation [17], and ventricular fibrillation [18,19].Typically such wave forms are described in the time domain, but identification of the major Fourier components could provide equivalent information that might be quantified more easily.Harmonicbased quantitative analyses of the wave forms associated with such oscillations can provide useful information about the underlying physiology: The existence of a nonlinear response function can be diagnosed by the detection of harmonics of the fundamental frequency, a principle that has been applied to electroretinography [10,20], to electroencephalography [21], and to the cardiovascular system [15].
In electrical engineering, methods based on spectral analysis are commonly used to identify the sources of such * aneta@lancaster.ac.uk harmonics so that they can be eliminated [22].Much the same problem arises in impedance spectroscopy [23].In such cases the fundamental driving frequency is known and constant and the primary problem is that of spectral estimation by using Fourier [6] or wavelet [9] analysis of the noisy data.An allied task is the estimation of the fundamental frequency from limited noisy data [5,24,25], which also occurs in mechanical vibration analysis [26].
In this paper we propose and validate an approach to the problem that allows identification and characterization of one or more nonsinusoidal oscillatory processes from a short time series of noisy data-data of the kind that often arise in practice, e.g., in the cases of vibrating machinery, electrical circuits, astronomical bodies, and physiological time series.Unlike earlier methodologies, which generally assumed a single underlying oscillatory process of fixed fundamental frequency, and/or that all spectral peaks that can be reliably identified in the signal belong to the same oscillatory process, the method proposed is also applicable to those cases where there may be two or more fundamental frequencies that are not necessarily constant.
Section II A sets the scene by reviewing briefly the methods hitherto available for identifying harmonics in cases where the underlying oscillatory process occurs at a fixed frequency.In Sec.II B we introduce our proposed approach, and it is then described in detail in Sec.III.The method is tested in Sec.IV, first on two different numerically generated time series in Secs.IV A and IV B, respectively, and then, in Sec.IV C, on physiological signals to demonstrate its robustness in the face of the noise and nonidealities present in real experimental time series.The results are discussed and summarized, and conclusions are drawn, in Sec.V.
where the fundamental frequency may be unknown, and so methods focus on accurate curve fitting or estimation of Fourier components [1,2], or on identifying a minimum in the entropy of the signal for a given trial period [3].A wavelet-based approach has been applied to the similar problem of spectral estimation that arises in the case of sunspot recordings [4], and this method is applicable to cases where the fundamental frequency or frequencies are time variable.Again, it is generally assumed that only a single independent oscillatory component is present, or that several such components can be distinguished easily.No method exists for distinguishing the harmonics of a time-variable fundamental frequency where unrelated processes at or near the harmonic frequency are known to exist (or are suspected).

B. An approach applicable when the period is variable
Where a time series contains only a few periods of oscillation and is subject to noise, where several significant oscillatory processes are superimposed with spectra that may overlap, and where the frequency and shape of the processes may be time variable, a different approach is required.The method of finding harmonics that we propose here is based on mutual information, by analogy with the Shannon entropybased sychronization index used by Tass et al. [27] to identify 1:n synchronization in bivariate data.Our comparison of phase time series differs from the approach used by Stogbauer et al. [28] to identify admixtures of signals in that, here, we specifically seek to demonstrate a relationship between different spectral components in order to better understand the underlying process responsible for the harmonics.
The Shannon entropy provides a measure of the unpredictability or information content of samples from a discrete distribution [29].The distribution (which may be conditional on another variable) of values has high entropy if it is uniform and corresponds to an unpredictable variable, and low if the distribution is sharply peaked and the variable can be predicted.The mutual information is a measure of the entropy that is "missing" from a conditional distribution once information about some underlying variable has been taken into account.
Cincotta et al. [3] assume monotonic phase growth of the signal fundamental, and test for significant mutual information between the distribution of signal values and a set of potential periodic functions.Rather than assuming a fixed signal period, for the case of oscillators with nonlinear phase growth, we instead introduce the wavelet transform to determine the phase at each point in time of each possible fundamental frequency in the data.After appropriate binning of the phase data (see below) to generate symbol sequences made up of discretized phase values, the mutual information of every possible pair of phase time series in the wavelet transform can be calculated.Cincotta et al. [3] have shown that the Shannon entropy of randomly and uniformly distributed signal values is expected to be a normally distributed variable for data of finite length.Within a short period of observation, however, because the phase time series extracted from the wavelet transform are intrinsically oscillatory, the phase series may have insufficient time to become uniformly distributed.Surrogate testing therefore must be applied instead to check for significant cross correlations.
The method of surrogate testing for mutual information measures has been developed by Paluš and Vejmellka [30] for the case of bivariate data representing two oscillatory processes that may, possibly, be causally related.Here we generate surrogate data in order to determine the distribution of mutual information values for independent processes with frequencies of near integer ratio, in order to confidently reject the null hypothesis when the actual value obtained reaches some threshold determined from the surrogate distribution.
We demonstrate the results of this algorithm on numerical examples incorporating time-variable frequencies and noise: first a single periodic process with harmonics, and then a signal incorporating two independent periodic processes with harmonics.Finally we apply the method to two physiological time series incorporating several independent periodic processes with harmonics, owing to a periodic perturbation in the first case and occurring spontaneously in the second.

A. Extracting phases
We consider a finite, noisy signal including one or more oscillatory processes with characteristic (but time-variable) frequencies and a nonsinusoidal profile, i.e., a shape giving rise to harmonics.We introduce the harmonic finder method to identify the harmonics amid spectral peaks owing to noise or unrelated processes.First the signal is detrended by use of a moving average to remove the mean and the spectral content outside the band of interest, and then a Morlet wavelet transform is applied [31].In practice, the transform used here will not provide reliable results for oscillation time scales longer than ≈1/6 the length of the time series, or shorter than twice the sampling period by the Nyquist criterion.The Morlet mother wavelet is a complex plane wave convolved with a Gaussian envelope function, For f 0 = 1 the characteristic frequency of the wavelet is 1 σ .The wavelet transform of a time series is the convolution of the complex wavelet [Eq.( 2)] with the time series at each scale σ , We use a wavelet of center frequency unity, rescaled in increments of 5% between 0.0025 and 2.5 Hz.The center frequency of the mother wavelet determines the time and frequency resolution of the transform at a given scale.The width of the wavelet in frequency space is proportional to the rescaled frequency ( f 0 σ ) of the wavelet, and the time scale is proportional to the rescaled period (σ ), thus giving logarithmic resolution.
This isolates the sinusoidal components of the signal and produces a matrix of complex values representing the phase of each frequency component at each point in time.In general, the phase growth is nonlinear.As the period of the underlying oscillation varies, the higher-frequency components vary in lockstep in order to preserve the shape of the wave form.Thus the phase time series corresponding to the fundamental or first harmonic will provide enough information to predict the phase time series corresponding to higher harmonics.Where the time series is short, low-frequency components in integer ratio will always have nonzero mutual information.For the limited time interval available, a spurious association may be found between phase values, because insufficient time is available for the (arbitrary) initial phase relationship found at the beginning of the time series to change.We test for a statistically significant excess of mutual information in the phase time series for different frequency components.

B. Finding the mutual information
The phase time series for each frequency component is extracted and discretized by equipartition into 24 bins, and then all the phase time series are compared pairwise.For each pair of phase time series the Shannon entropy [29] of the phase distribution p (φ 1 ) of the higher-frequency component is determined, and then the mean entropy of the distribution given the corresponding phase of the lower-frequency component for each moment in time, p (φ 1 φ 2 ).The number of bins was chosen such that peaks in p (φ 1 φ 2 ) would each fall into a single bin for frequency ratios of 1:2, 1:3, and 1:4.By Nyquist's criterion, the maximum ratio that remains detectable in theory is 1:12, but the presence of noise in the data generally prevents this limit from being attained.Because the distribution obtained is subject to noise, use of a larger number of bins may not provide higher resolution and can become counterproductive owing to the smaller number of samples available to determine the density of the distribution in each bin.Thus the choice of 24 for the number of bins represents a compromise.
The entropy of the discretized phase values is and the mean entropy of the conditional distribution is The mutual information content of the signals is equal to the difference The proportion of mutual information for φ (ω 1 ,t) and φ (ω 2 ,t) approaches 1 as ω 1 approaches ω 2 .We also applied the k-nearest-neighbor estimator of mutual information described by Kraskov, Stogbauer, and Grassberger [32] to our data and found that it gave results comparable to those obtained from the method of binning.In both cases surrogate testing is required because the mutual information values are strongly biased by the correlations present in the relatively short phase time series, not just by the discretization into bins.In practice, we preferred to use binning because of its greater computational speed, which was useful when evaluating many surrogates.

C. Generating surrogates
In order to test for statistical significance the (much lower) 2:1, 3:1, . . ., n:1 mutual information values that we obtain, we need a surrogate distribution of mutual information for time series that have the same spectral content, but in which the oscillations at different frequencies are known to preserve no particular phase relationship.The use of surrogates avoids the problem of having to calculate the probability of a given M (φ 1 ,φ 2 ) value for phase time series with known long autocorrelation times, rather than the random distribution assumed by Cincotta et al. [3].
We can generate a set of surrogate data by randomizing the phases of the Fourier components of the original signal and transforming back into the time domain.In practice, the distribution of signal values may be non-Gaussian, in which case this procedure may substantially alter the distribution of the data.To avoid this, we sort and then map the signal values onto an equivalent number of sorted values drawn from a normal distribution.Then we Fourier transform this normalized signal, phase randomize the components, and inverse transform to recover a surrogate with a Gaussian distribution, which we map back onto the values from the original time series.These are referred to as amplitude-adjusted Fourier transform (AAFT) surrogates [33].To generate a surrogate mutual information value for each frequency pairing, we determine the phase series mutual information of the putative harmonic frequency drawn from the surrogate data, and the putative fundamental frequency drawn from the original data.
The wavelet transform is a function of time incorporating all the Fourier components around each wavelet frequency.In order for this surrogate approach to be effective, any spectral peaks in the data must be sufficiently broad that phase randomization of the Fourier components does not merely result in a fixed phase shift, but a detectable change in the phase dynamics as determined by the wavelet transform.There must be enough information in the phase growth time series to randomize.
Kralemann et al. [34] have described the importance of considering the phase dynamics of a system of nonlinear oscillators separately from the wave form associated with one projection of the system, and the possibly nonlinear growth of a "protophase" extracted from such a variable.The process of filtering the signal by use of the wavelet transform, and the subsequent equipartitioning of the phase time series, largely avoids this issue.The procedure for checking for one or more fundamental oscillations and their associated harmonics is predicated on the assumption that the shape of the wave form itself, and not just the phase dynamics, may provide physically or physiologically relevant information about the system.

D. Testing significance
Finally, a local maximum in the mutual information for a pairing of phase time series is deemed to indicate a harmonic if it occurs some number of standard deviations above the expectation value of the mutual information distribution for the uncorrelated surrogates.In practice, the number of standard deviations required should be adapted to take into account both the acceptable rate of false positives (typically no more than 5%) and the effect of multiple testing on this false positive rate.
The wavelet transform corresponds to a filter in the frequency domain, and the lower the center frequency of the mother wavelet (f 0 ), the lower the frequency resolution, resulting in a broader range ω 1 -ω 2 for which the neighboring phase time series are dominated by a single component.Conversely, greater frequency resolution results in fewer independent samples of phase and poorer statistics.

IV. TESTS AND APPLICATIONS
Before applying the method to real physiological data, we test it on numerically generated time series of known harmonic content.

A. A numerically generated square wave
To begin, we generate a random-walk noise signal of length 2000 s, sampled at 100 Hz.This signal B is the cumulative sum of 200000 randomly selected values w(τ ) from a normal distribution N of unit variance: where The random-walk signal B is detrended by subtraction of a moving 400-s average to give stationarity over long time scales.This Brown-noise time series has power at all frequencies above that corresponding to the detrending time scale and serves to obscure our test signals and their harmonics from easy detection.Note that the use of additive white Gaussian noise would be inadequate for this application as a simple smoothing-resampling of the time series would be enough to reveal the underlying test signal.
As a first test of the harmonic finder method, we generate an artificial, approximately periodic signal incorporating harmonics, and hide it under the Brown noise.Figure 1(a) shows the noise time series superimposed on a square wave S of amplitude five units, which itself has a variable period.The phase of the square wave φ grows at a mean rate of 1 cps, plus a variable term equivalent to a 10% standard deviation, generated from another detrended Brown-noise time series (as above) rescaled by mapping onto samples from a normal distribution N(1,0.01).This produces slow changes in the rate of phase growth throughout the time period of observation, and corresponding changes in the period of the square wave: The square wave and its harmonics are not obvious in the Fourier transform of the noisy signal shown in Fig. 1(b), but we can see a possible harmonic at 3 Hz in the wavelet transform [Fig.1(c)].We wish to determine rigorously whether this oscillation is a harmonic of the 1-Hz signal or an independent higher-frequency oscillation.Figure 1(d) shows the mutual information of each phase time series pairing: Note that this plot is of course dominated by the 1:1 relationship between the phases of wavelet components of close or identical frequencies.To test the mutual information values obtained at frequency ratios of 2:1, 3:1, etc., we generate a set of AAFT surrogates for the original signal, as described above.
This destroys the phase relationships between the harmonics of the square wave, while leaving the spectrum largely unchanged.Note that the existence of such surrogates demonstrates that the identification of harmonics by checking for spectral peaks is inadequate, as the surrogates no longer include a square-wave component despite maintaining the same spectrum as the original signal.The mean mutual information values for the phase time series of the original signal and the phase time series of 1000 such surrogates are shown in Fig. 1(e).We see high but spurious values at spectral peaks and for low frequencies, as expected owing to the lower number of independent measurements of phase when autocorrelation times are long compared to the time scale of simulation and/or measurement.
Figure 1(f) shows the actual mutual information values relative to the surrogate distribution, which we will refer to as an M-plot.Both the third and fifth harmonics are clearly identifiable.A threshold of five standard deviations above the mean of the surrogate distribution is marked in blue, and the local maxima found at 3 and 5 Hz are marked by pluses (+).Five standard deviations were chosen as the criterion to ensure that the 1596 tests performed (comparing all possible pairs of phase time series to their surrogate distributions) would not be expected to produce false positives for a normal distribution of surrogate mutual information values.The method has evidently succeeded in this particular case.
Note that, if the square-wave signal frequency was without variability, no meaningful calculation of mutual information could be performed because the phase relationship (and the complete phase time series) would be determined completely by the initial phases of the two components.A random-walktype Brown-noise signal was used to introduce this variability, ensuring that the variability was significant over long time scales (longer than several times the averaging time associated with the Morlet wavelet envelope).In practice, this is a realistic pattern for the medium-term variability of physiological data, as presented below.
Low integer ratios between harmonics are more easily detected than high ratios, with the result that, in the case of a square wave, the M plot may correctly find significant associations between the third and ninth harmonic (not shown), but not the first and ninth harmonic.

B. Two numerically generated square waves
Next, the method is applied to a more complex signal, consisting of a Brown-noise time series as before, plus a pair of uncorrelated square waves with amplitudes of ten units and mean frequencies at 1 and 2 Hz, as plotted in Fig. 2(a).The aim here is to test the ability of the method to detect independent periodic processes and correctly attribute their harmonics to them.The 1-and 2-Hz components apparent in the Fourier transform might be interpreted naively as the fundamental and harmonic of a single oscillatory process.In the wavelet transform, more components become visible, but the relationship between them is far from obvious.After the mutual information of the phase time series for each pair of frequencies has been calculated and compared to the distribution obtained with AAFT surrogates, the components of two independent oscillations become distinguishable.Crucially, there is no detection at (1 Hz, 2 Hz) in the M plot of Fig. 2(f), indicating that these components are independent.A maximum at (1 Hz, 3 Hz) is found as in case I, together with the maxima at (2 Hz, 6 Hz) corresponding to the second independent square wave.The (1 Hz, 5 Hz) and (2 Hz, 10 Hz) maxima fall below the 5σ threshold for detection in this complicated signal, but as expected they reappear if more data are included or if a larger square-wave amplitude is used.
Having identified two independent fundamental frequencies in the data, we can reconstruct the approximate shape of each.First we detrend the signal to remove energy owing to noise below the fundamental frequency.Then, after equipartition of the fundamental phase values, we average over the signal values found in each bin to remove energy owing to noise from above the fundamental.The error in the mean wave-form value for each phase bin is dominated by the error in the phase discretization procedure, which tends to smooth the wave form by assigning neighboring but distinct values to the same bin, or by simply incorrectly determining the phase value associated with a certain signal value where the signal is dominated by noise.Figure 3(a) shows 2 s of the signal after detrending (in which the square-wave oscillations are not at all apparent) and Figs.3(b) and 3(c) show the reconstructions obtained by using the phase time series of the two fundamental frequencies, in which we see that both square waves can be independently reconstructed from the data with approximately the correct amplitudes.Accuracy is not improved by the addition of more data owing to the systematic misattribution of phase values caused by the noise: Only a reduction of the effective noise in the phase time series, e.g., by use of a mother wavelet with a higher center frequency, could improve the accuracy of the reconstruction for this data.

C. Skin blood flow
In human subjects, noninvasive measurements of skin blood perfusion by laser Doppler flowmetry (LDF) enable the investigation of several linked processes in the cardiovascular system.Each has its own characteristic time scale of variability, including cardiac, respiratory, myogenic, neurogenic, and endothelial processes [35].The characteristic frequencies of these processes are known approximately, but vary from subject to subject and in time.Harmonics of low-frequency oscillations may be difficult to identify and, consequently, it may be difficult to determine correctly the signal energy associated with each process.The signal is also subject to noise.This represents a useful test case in which to apply our method of statistical testing for correlated harmonic content.
A DRT4 laser-Doppler flowmeter with DP1T-V2 probe (Moor Instruments, Axminster, Devon, U.K.) was used to record subcutaneous blood flow noninvasively.Flow is expressed here in (arbitrary) standard perfusion units (SPU), proportional to the density of red blood cells in the sample volume, and to the mean velocity of the cells [36,37].
We consider 1800 s of laser Doppler skin blood-flow data, subject to both spontaneous oscillations and induced lowfrequency oscillations owing to the periodic administration (locally, by iontophoresis) of a vasodilatory drug: This is a standard procedure, and the resultant signal is typical of LDF measurements of this kind that require analysis.Acetylcholine was administered to the skin of the left arm for a period of 20 s, commencing every 260 s, over a 30-min period of measurement.Details of the procedure are given by Shiogai, Stefanovska, and McClintock [38].
For this case, the wavelet time resolution must be relatively good because the low-frequency components may have highly time-variable frequencies.The rescaling increment is chosen appropriately for the associated frequency resolution of the wavelet.Note that the number of cardiac oscillations in the time series is comparable to the number of square-wave periods in the numerical examples above.Without making use of prior knowledge of which frequency components may produce harmonics, we examine the wavelet transform across the whole spectrum from the cardiac frequency downward.
The raw signal during intophoresis for one subject is shown in Fig. 4(a).The period and shape of the iontophoretic oscillation in blood flow is apparent from the raw signal and its wavelet energy spectrum in Fig. 4(b).By visual inspection, we see that variability in the signal is dominated by the cardiac pulse oscillation and the periodicity of the iontophoretic protocol.The corresponding raw M plot identifies the harmonics of the protocol and also the harmonics of the natural cardiac oscillation as shown in Fig. 4(c).Because there is no intrinsic variability in the iontophoretic protocol administered by the machine, but there is intrinsic variability in the heart rate, comparison with surrogates reveals that only the maxima associated with the cardiac component are significant.At the iontophoretic frequency, the AAFT surrogates have the same uniform phase growth as the actual data, and thus the same mutual information with the harmonics of this process, illustrating an unavoidable limitation of this form of testing.
In fact, both the cardiac wave form and the wave form associated with periodic iontophoretic vasodilation can be reconstructed by using the binned wavelet phase series, as shown in Fig. 4(e).The cardiac wave form corresponds to the classic blood-flow wave form associated with the cardiac pressure wave traveling through the microvasculature.The wave form associated with the iontophoretic protocol matches the pattern visible to the eye in the raw data.
Finally, we apply the method to control measurements that were made on undisturbed skin as part of another study.This demonstrates the existence of low-frequency harmonic processes independent of any periodic stimulus.A sample of the signal is shown in Fig. 5(a), its wavelet energy spectrum is shown in Fig. 5(b), and the surrogate-adjusted M plot is shown in Fig. 5(c).We clearly see harmonics of the oscillatory processes in the cardiac, respiratory, and myogenic bands, although the fundamental frequencies are less clearly defined.Focusing on the second harmonic of each frequency, we see peaks in Fig. 5(c), indicating harmonics at the cardiac (1 Hz), respiratory (0.25 Hz), and myogenic (∼0.07 Hz fundamental) frequencies.The M plot is most complex around the myogenic frequencies, with two sets of two peaks at (0.07 Hz, 0.15 Hz), (0.08 Hz, 0.16 Hz) and (0.10 Hz, 0.20 Hz), (0.11 Hz, 0.22 Hz).We reconstruct the form of each oscillation in Fig. 5(d), taking the myogenic fundamental to be 0.07 Hz.
The presence of several overlapping sets of harmonic oscillations in the myogenic band indicates complicated oscillatory behavior, possibly including significant rate variability of the underlying process.The reconstructed wave form is sharply peaked, and this may provide a clue to the origin of the harmonics.Oscillations of the smooth muscle in arterioles are thought to alter the vascular resistance by widening and narrowing the vessels.By Poiseulle's law, the resultant flow is expected to vary as the fourth power of the radius, and this nonlinear relationship may account for some of the high harmonics corresponding to the nonsinusoidal nature of the myogenic oscillation profile.

V. SUMMARY AND CONCLUSION
The application of the method to both simulated data and real-time series demonstrates its usefulness in detecting the harmonics of nonsinusoidal oscillations in relatively short segments of noisy data.Where the characteristic profile of a particular oscillatory process is visible by the eye, the use of a statistical method with surrogate testing may help to validate the observation objectively.Where the profile is unclear this method may open up new ways of analyzing the data, if the source is otherwise considered as a "black box," and if it is unknown what relationships may exist between the spectral components of its output.
While the harmonic content of nonsinusoidal oscillations is of intrinsic interest, in other circumstances the existence of unrecognized harmonics in a time series thought to consist of independent oscillatory processes may cause statistical distortions.For example, when performing wavelet coherence analysis, an areawise test of coherence [39] may give false positives if a single large coherent patch is formed by spuriously coherent oscillations and their harmonics (which would be treated as independently testable oscillations by surrogating techniques, unless the harmonics had been identified).Identifying the harmonics in the data set enables us to eliminate this possibility.
The harmonic finder method introduced in this paper can be used to identify the harmonics of noisy periodic processes with time-variable frequencies, such as oscillations in skin blood flow (both stimulated and spontaneous).The use of the wavelet transform enables us to determine the phase of each component at each point in time, and use of the Shannon entropy enables us to identify correlations between putative harmonics.When frequencies are time variable, surrogate testing allows us to identify significant correlations.We can thus test for harmonics without the use of a fixed trial period and hence determine statistical significance while making minimal assumptions about the data.

FIG. 1 .
FIG. 1. (Color) Analysis of a numerically generated noisy squarewave signal.(a) The raw signal (arbitrary units), (b) its Fourier transform, (c) its wavelet transform, (d) the corresponding M plot, (e) the mean M plot for AAFT surrogates, and (f) the M plot relative to the mean and standard deviation of the surrogates.Regions more than 5σ above the surrogate mean are marked with blue, and local maxima are marked with +.

Figure 2 (
b) shows the Fourier transform, Fig. 2(c) shows the wavelet transform, and Figs.2(d)-2(f) show the M plots for this more complicated case.

FIG. 2 .
FIG. 2. (Color) Analysis of a numerically generated noisy time series containing two square-wave signals.(a) The raw signal (arbitrary units), (b) its Fourier transform, (c) its wavelet transform, (d) the corresponding M plot, (e) the mean M plot for AAFT surrogates, and (f) the M plot relative to the mean and standard deviation of the surrogates.Regions more than 5σ above the surrogate mean are marked with blue, and local maxima are marked with +.

FIG. 3 .
FIG. 3. (a) A 2-s sample of the noisy two-square-wave numerical data (arbitrary units) after detrending.(b) The reconstructed 1-Hz wave form (black curve) and the profile of the original square wave (gray curve).(c) The reconstructed 2-Hz wave form (black curve) and the profile of the original square wave (gray curve).Bin numbers should be multiplied by 2π/24 to recover the phase values in radians associated with each bin.

FIG. 4 .
FIG. 4. (Color) (a) A detrended blood-flow signal from human skin subjected to the periodic iontophoresis of a vasodilatory agent (in SPU).(b) The Morlet wavelet transform spectrum of the signal (in SPU 2 /s).(c) The corresponding M plot.(d) The M plot relative to the surrogate distribution.Regions more than 5σ above the surrogate mean are marked with blue.Because the iontophoresis period is fixed, the maxima at this frequency in the M plot also occur in the surrogates.(e) Reconstructed wave forms of the iontophoretic response (left) and the cardiac pulse (right).

FIG. 5 .
FIG. 5. (Color) (a) A detrended blood-flow signal from undisturbed human skin (in SPU).(b) Its Morlet wavelet transform spectrum (in SPU 2 /s).(c) The M plot relative to the surrogate distribution.Regions more than 5σ above the surrogate mean are marked with blue, and local maxima are marked with +.(d) Reconstructed wave forms of (left to right) the myogenic and respiratory oscillations, and the cardiac pulse.