DPOAE GFs were measured from the right ears of 12 adult participants (9 female, 3 male; mean = 22.9 yr, SD = 1.23 yr, range = 20–25 yr). Participants reported a negative history of hearing problems and listening difficulties. Participants’ ear canals were free of cerumen accumulation, and eardrums were visually normal by otoscopic inspection. Middle ear function was normal as assessed by 226-Hz tympanometry. Participants passed a DPOAE screening (both primaries swept continuously in frequency, f2 = 0.75–5 kHz, f2/f1 = 1.22, L1 = 65 FPL, L2 = 55 dB FPL), with peaks in amplitude fine structure ≥ 10 dB SPL from 1–4 kHz. Audiometric thresholds were not tested, since the study was focused on participants with reasonably strong DPOAE magnitudes, rather than with clinically normal behavioral thresholds, per se. The research protocol was approved by the University of Iowa Institutional Review Board, and written informed consent was obtained from all participants.

Stimuli were presented and ear canal responses measured using an ER10X probe system (Etymotic Research). The probe system was connected to a 24-bit sound card (RME Fireface 802). Playback and recording were controlled by custom MATLAB software [31]. The sampling rate was 96 kHz, and recordings were streamed to computer disk for post-hoc analysis. A critical part of obtaining high-quality DPOAE recordings is ensuring probe stability. The ER10X probes were inserted into the canal and then stabilized by attaching the probe cord to a headband (3M H10A Peltor Optime earmuffs with the outer plastic and foam core removed to leave a donut shape surrounding the pinna). Thevenin source calibration [32] was completed for each ER10X probe prior to data collection. In each ear, an in-situ calibration routine calculated the characteristic and load impedances. The loudspeaker voltage drive was adjusted to deliver the stimulus at the desired forward pressure level (FPL) at each test frequency. For the longer, modified test paradigm (described below), the in-situ calibration routine was repeated every 3.6 min during the recording session, and the output voltage was adjusted to ensure stable stimulus pressures at the eardrum.

DPOAEs were elicited with L1 fixed at 65 dB FPL while L2 was swept continuously from −2 to 72 dB FPL [15]. Each tone sweep presentation was 9 s in length; therefore, L2 varied at a constant rate of 8.22 dB/s. These rate and duration parameters were chosen to minimize analysis edge effects while maximizing analyzed DPOAE SNR. Each test consisted of \(M=\) 24 nine-second presentations, for a total length of 216 s. DPOAE GFs were collected for each participant at a single f2 frequency between 1–4 kHz (mean = 1817.3 Hz; SD = 847.1 Hz; min = 1056.0 Hz; max = 3881.0 Hz). To maximize SNR and to avoid potential GF notches due to two-component interactions, the value of f2 was chosen from a peak in the DPOAE amplitude fine structure obtained during the hearing screening [33]. Primary tone ratio was fixed at f2/f1 = 1.22.

As described in the next section, each GF obtained from the swept stimulus paradigm was fitted with a smooth curve that captured the essential features of peak regions and predicted DPOAE levels below the noise floor. To test the accuracy of these predictions at low levels, the same 12 participants were tested using a second, modified test paradigm. This modified paradigm included focused averaging at L2 = 5 dB FPL (L1 = 65 dB FPL), a level which produced DPOAE amplitudes below the noise floor of the swept paradigm. This single “snapshot” at L2 = 5 dB FPL provided a representative test of low-level prediction accuracy. In the modified paradigm, two fixed L1 sweeps were presented, followed by 2 min of continuous averaging with L1 held constant at 65 dB FPL and L2 fixed at 5 dB FPL. The interleaving process was repeated 13 times for a total of 26 interleaved sweeps and 26 min of averaging time at the low level (~ 39 min of total recording time). After the recording was completed, the swept and steady portions were de-interleaved and analyzed separately. The collection of swept recordings was analyzed as described in detail below. The steady low-level portion was analyzed using the same basic signal processing methods. The long averaging at L2 = 5 dB FPL reduced the noise floor by approximately 21 dB, yielding DPOAEs with approximately 10 dB SNR in most participants. Probes were not removed from participants’ ears between tests.

Recordings were analyzed using ordinary least-squares regression applied to a series of partially overlapping time-windowed waveforms. Each analysis frame was 500 ms in length. Seventy-one analysis frames were considered for each sweep, with each frame centered at a temporal location corresponding to a target L2 FPL value. Target values ranged from 0–70 dB FPL in 1 dB steps (75% temporal overlap between adjacent frames).

Let the variable \(}\) represent a column vector of recorded ear canal pressures in an analysis frame (\(}\) of length n = 48,000 samples), multiplied by a Hann window. The overdetermined system of equations \(}}=}\) was considered, where

$$}=\left[\begin _ \\ _\\ \vdots \\ _\end \begin _ \\ _ \\ \vdots \\ _\end\begin_\\ _\\ \vdots \\ _\end\right], }=\left[\begin_\\ _\\ _\end\right] , }=\left[\begin\begin_\\ _\\ \vdots \end\\ _\end\right]$$

(1)

The first and second columns in \(}\) contained Hann-windowed sinusoids of frequency 2f1-f2 in cosine and -sine phase, respectively. The third column in \(}\) was a vector of ones, to account for any mean offset from zero. The vector \(}\) contained the three unknown coefficients corresponding to the three columns of \(}\). The best-fitting (in the least squares sense) coefficients, \(\widehat}}\), were found by minimizing the matrix norm

$$\widehat}}=\underset}}}}-}}\Vert }^$$

(2)

The estimate of DPOAE magnitude and phase in that analysis frame was expressed in complex rectangular form as

where \(i\) is the imaginary operator. In this paper, variables written in script font, such as \(\), indicate the variable is complex. The fitting procedure was repeated for each of the \(M=\) 24 stimulus sweeps in the same analysis frame, and the complex coefficients were stored in a vector \(}\) of size 24 × 1 (\(M=\) 26 in the modified test paradigm).

In calculating DPOAE signal and noise levels from \(}\), a weighting strategy was used to reduce or remove the influence of any elements contaminated by intermittent noise. For relatively long signals such as these nine-second sweeps, this strategy allows effective artifact control at the resolution of single analysis frames (500 ms), avoiding the need to discard entire sweeps. An iterative bisquares reweighting algorithm [34] was applied to the magnitudes of the 24 complex coefficients in \(}\) to obtain a column vector of corresponding weights \(}\). This method minimizes the weighted sum of squares, with the weight assigned to each data point dependent on the distance of that point from the weighted mean in the previous iteration. Data points near the weighted mean receive full weight, while data points farther from the weighted mean receive reduced weights. The process of reweighting is repeated iteratively until the weighted mean converges (i.e., the fractional change in the current weighted mean from the previous weighted mean reaches some criterion amount). The final weighted mean minimizes the effect of outliers. For ease of subsequent computation, the final weights were normalized by dividing each element of \(}\) by the sum of \(}\), so that the weights summed to one. The DPOAE signal pressure was then calculated as the weighted complex mean across sweeps:

$$\overline}=}}}^},$$

(4)

where superscript \(\mathsf\) is the transpose operator. For creating GFs, DPOAE (weighted) mean magnitude was calculated as

$$\tt=\sqrt\text}\right)}^+\text\text }\right)}^},$$

(5)

where \(\tt\) written in Latin font indicates the variable is strictly real. The DPOAE noise floor was calculated as the weighted estimator of the standard deviation of the mean (the “standard error of the mean”) across sweeps:

$$}_}=\sqrt^}^ }\right)\sum }\left(\left(}-\overline }\right)}-\overline }\right)}^}\right)},$$

(6)

where \(\overline }\) is the weighted complex mean of \(}\), and the asterisk in Eq. 6 indicates complex conjugation. \(}\) is the normalized vector of weights described above, and \(M\) is the number of sweeps.Footnote 1 In Eq. 6, the vectors are multiplied element-wise.

The analysis procedure described above was applied to the 71 analysis frames in the recorded sweep, centered at target values of L2 ranging from 0 to 70 dB FPL. In Eqs. 5 and 6, the variables \(\tt\) and \(}_}\) refer to values obtained from a single analysis frame. To be clear: \(}\) is a complex vector across sweeps at a single stimulus level, the scalar \(\tt\) is the mean magnitude of \(}\), and throughout the rest of this paper, \(}\) is a real vector composed of \(\tt\) at each analyzed stimulus level. Similarly, the real vector \(}_}\) is composed of \(}_}\) across stimulus levels. In other words, \(}\) and \(}}_}}\) are 71 × 1 vectors obtained by analyzing the entire recording. Where needed, subscripts are used to indicate single elements in \(}\) (e.g., \(}}_\)). The best point estimate of the noise floor was taken as the mean across \(}}_}}\), written \(\overline}_}}\). To create an GF, \(}\) was plotted as a function of f2 stimulus magnitude.

The GFs obtained from 12 participants using the standard swept paradigm are shown in Fig. 1. Figure 1a shows the data on log-log axes, which is helpful for visualizing the lower L2 levels. Figure 1b shows the same data on linear-linear axes, which is the form in which the data were analyzed and fitted. The linear scale clearly shows the compressive growth of the GF. In the fixed L1 paradigm, GFs grow to a peak and then decrease. The decrease in the GFs at higher L2 likely occurs because of two-tone suppression effects and/or changes in the way scattered wavelets combine inside the cochlea as L2 approaches L1 [14, 35]. In general, GFs appeared to decrease asymptotically above the peak. For some participants, the GF continued its downward trajectory for L2 > 65 dB FPL; however, for other participants, the GF became flatter, or in some cases began to increase slightly. The cause of these differences is unclear, and additional experiments with L2 going to higher levels are needed to understand these observations. For the remainder of this report, GF data are considered for L2 from 0–65 dB FPL.

DPOAE GFs for 12 participants. Panel a shows data on a log-log scale; Panel b shows the same data on a linear-linear scale. Both Panels: Mean DPOAE magnitudes (\(}\)) are shown by blue lines with dot markers. Noise floors, calculated as the standard error of the mean (\(}}_}}\)), are shown by black lines with downward pointing triangle markers. Stimulus levels directly measured at the probe tip can be inaccurate due to standing waves in the canal; therefore, in Panel a, measured DPOAE levels (Ldp) are expressed as a function of target L2 values in dB FPL instead of measured L2 values in dB SPL. In Panel b, the quantities are expressed as linear pressure magnitudes (\(_|\) and \(_|\))

It is of interest to characterize observed GFs with a statistical model. Models reduce the dimensionality of the data, allowing a succinct description using a relatively small number of parameters. They also facilitate quantifying GF characteristics, such as slopes and thresholds, which have been used to assess cochlear function. Under the assumption that DPOAEs are causally related to an underlying nonlinear cochlear vibration pattern, we hypothesized that the observed GFs could be characterized by a small number of parameters belonging to a compressive nonlinear growth function representing this vibration pattern. This nonlinear growth function is assumed to correspond to a physical reality (vibration of cochlear structures) that could, in principle, be measured but currently cannot be noninvasively observed in humans. We constructed a latent growth statistical model which related the set of unobservable latent variables (the non-linearity parameters) to observable variables (the GF and noise floor estimate) via a set of structural equations and a generalized linear regression. A block diagram of the model is shown in Fig. 2. The latent portion of the model is shown inside the ellipse. The stimuli (DPOAE primaries) were passed through the non-linearity, and the output of the non-linearity was analyzed in the frequency domain to obtain a resulting predicted GF. The optimization loop shows the process by which the values of the coefficients for the non-linearity were adjusted to produce the predicted GF that best fitted the observed GF data.

Block diagram of the structural model. The primary input to the model is the variable triple of parameter values for the latent non-linearity (\(\\)). A secondary input is the fixed pair of primary (stimulus) waveforms used to evoke DPOAEs (\(}}_,}}_\)). The primary waveforms are passed through the nonlinearity, resulting in a distorted waveform (\(}\)), which is analyzed to extract magnitude at the 2f1-f2 frequency (\(}\), shown in red font). The vector \(}\) is the output of the model and is referred to as the “predicted GF”. An iterative loop is used to adjust the nonlinearity parameter values to produce the optimal match between the predicted GF and the observed GF data measured from an individual ear (\(},\overline_}}\), shown in blue font). Optimization is based on maximum likelihood. A generalized linear regression model is used, wherein each GF value is conceptualized as a random draw from the Rician distribution. Rician parameter values are dependent on SNR, which varies across the GF. The Rician parameter values (\(\},}\}\), shown in purple font) associated with each element in the predicted GF (i.e., each element in vector \(}\)) are determined using link functions. The likelihood of \(\},}\}\), given the observed GF (\(}\)) is calculated. For an observed GF measured from an individual ear, the optimal values for the latent nonlinearity triple are those producing the predicted GF with corresponding Rician parameters having the maximum joint probability, given the observed GF

There are multiple sources of cochlear non-linearity that may contribute to DPOAE generation, including processes associated with mechanoelectrical transduction (e.g., [36]) and electrical-to-mechanical transduction (e.g., [37]). The non-linearity used in the structural model does not necessarily represent any of the specific non-linearities associated with DPOAE generation, but it may approximate a composite non-linearity that captures the effective contribution of individual mechanistic non-linearities. Previous work has described the conversion of sound stimuli into cochlear responses using a variety of saturating non-linear functions with sigmoidal shapes, such as arctangent, hyperbolic arctangent, logistic, first- or second-order Boltzmann, and third-order polynomials [21,22,23,24,25,26,27]. We used a simplified form of the generalized logistic function [38] as the non-linearity for the latent portion of the statistical model.

The generalized logistic function has six free parameters, providing great flexibility. In preliminary analyses, we found that we could accurately fit our observed GF data using a modified version of the generalized logistic function, wherein the number of free parameters was reduced to two (\(\beta\) and \(\gamma\)), plus one additional stimulus scaling parameter \(\kappa\). The simplified generalized logistic function had the form

$$}=-2+\frac^}}\right)}^}} ,$$

(7)

$$}=\kappa }}_+}}_ .$$

(8)

In Eq. 7, the scalar values −2 and 4 set the minimum and maximum outputs of the function at ± 2. These values were chosen for convenience in working with acoustic data (representing 2 Pa, or 100 dB SPL, a realistic range of sound pressures), but the goodness of the obtained fits is not dependent on these specific values. Output of the non-linearity was a function of instantaneous stimulus sound pressure (\(}\)). The vector \(}\) was the summed pressures of the two primaries, \(}}_\) and \(}}_\) (having frequencies \(_\), \(_\) and levels \(_\), \(_\), respectively). The third free parameter, \(\kappa\), scaled the lower frequency primary prior to summing with the higher frequency primary, as shown in Eq. 8. The waveform outputs of the instantaneous nonlinearity (\(}\); Eq. 7) were converted to the frequency domain, and the magnitudes at frequency 2f1– f2 were used to create a predicted GF (\(}\)). Figure 3 shows the effects of the three parameters on the nonlinearity (Panels a, c, and e) and the 2f1– f2 GF (Panels b, d, and f). With regard to the nonlinearity, the parameter \(\beta\) controls the steepness of growth, and the parameter \(\gamma\) controls where maximum growth occurs (\(\gamma\) >1 moves the maximum growth towards the upper asymptote).The parameter \(\kappa\) has no effect on the nonlinearity itself, and its main effect is to linearly increase or decrease the overall GF magnitude.

Effects of the three parameters on the shape of the non-linearity (Panels a, c, and e) and the shapes of the resulting predicted GFs (Panels b, d, and f). The effects of each parameter are shown in pairs of panels. Panels a, b show the effects of changing \(\beta\) with the other parameters held constant. Likewise, Panels c, d and Panels e, f show the effects of changing \(\gamma\) and \(\kappa\), respectively, with the other parameters held constant. When held constant, the values used were the mean of the fitted values from the 12 participants in this study (Table 1). In all panels (except e), curves are shown for 20 selected values of the parameter of interest. Increasing values of the parameters are indicated by colors changing from blue to red, and the main directions of GF change are indicated by arrows. Panels a, b show the effect of increasing the value of \(\beta\) is to steepen the slope of the nonlinearity (a) and move the peak of the GF up and to the left (b). Panels c, d show the effect of increasing the value of \(\gamma\) is to move the location of maximum growth towards the upper asymptote of the nonlinearity (c) and move the peak of the GF down and to the right (d). Panels e, f show that \(\kappa\) has no effect on the nonlinearity since it does not appear in Eq. 7 (e); however, the effect of increasing the value of \(\kappa\) is to move the entire GF up (f), i.e., expressed on a log–log scale, the GF simply starts and ends at a higher y-axis value with increasing \(\kappa\), while the shape remains constant

Producing the predicted GF (\(}\)) that best matched an observed GF (\(}\)) required optimizing the parameter triplet \(\\). The parameter values were estimated using likelihood and the probability density function (PDF) of the Rician distribution. GFs are generally obtained in such a way (i.e., the averaging time is held constant for each stimulus level) that the DPOAE magnitude varies while the underlying noise floor remains stationary. As a result, SNR varies along the x-axis. If all data points in the GF have an SNR greater than approximately 10 dB, the signal + noise magnitude can be modeled statistically using the Normal distribution. However, when DPOAEs have less than approximately 10 dB SNR, the associated signal + noise magnitude is no longer normally distributed. In the limit of no DPOAE present, the signal + noise magnitude is Rayleigh distributed. At low SNRs (i.e., \(\text\text\text\lesssim 10 \text\text)\), the signal + noise magnitude varies systematically between the Rayleigh and Normal distributions, and the Rician distribution quantifies this changing distribution. The Rician distribution is well known in various fields, including the signal processing of noisy images, such as those obtained by MRI. It has been occasionally used to study otoacoustic emissions (e.g., [39]).

Formally, the Rician distribution describes the magnitudes of circular-symmetric bivariate normal random variables, including those with non-zero means. Expressed in terms of the variables used in this paper, the PDF of the Rician distribution is

$$f\left(\tt;\nu,\sigma \right)=\frac}^}^}^+^\right)/2^\right)}_\left(\frac\nu}^}\right),$$

(9)

where \(f\left(\tt;\nu,\sigma \right)\) is the probability density of the DPOAE magnitude (\(\tt\)), given the Rician parameters \(\nu\) and \(\sigma\). \(_\) is the modified Bessel function of the first kind with order zero. The distribution has two parameters: \(\nu\), which controls the location, and \(\sigma\), which controls the scale (spread). To use the Rician distribution in the context of fitting predicted GFs to observed data, additional functions are needed to link each predicted GF value (each element of \(}\)) to its Rician parameter values, given the expected SNR. In other words, the linking functions were needed because the parameter values of the Rician distribution depend on SNR, which in this study varied as a function of L2. Given the predicted GF value (“signal”) and the measured noise floor (“noise”) at a given L2, the linking functions returned the appropriate SNR-based Rician parameter values. We defined two functions mapping SNR to the parameter values of the Rician distribution. Details are given in Appendix A. Applying these functions to \(}\) yielded two vectors \(}\) and \(}\), with elements which correspond both to \(}\) and the associated vector of primary levels. Maximum likelihood estimation was used to determine the most likely values of vectors \(}\) and \(}\), given an observed GF. The likelihood function was

$$\mathcal\left(},};}\right)=\prod\limits_^f(}_;_,_),$$

(10)

where \(\mathcal\left(},};}\right)\) is the likelihood of the parameter vectors \(}\) and \(}\), given the observed GF data in the magnitude vector \(}\), and \(N\) is the number of analysis frames. Note that although the likelihood is computed in terms of \(}\) and \(}\), these values were linked deterministically to the predicted GF via SNR, so that the likelihood of the set of coefficients \(\\) was indirectly determined.

Maximum likelihood estimation is an optimization problem, wherein the values of the parameters that maximize the joint probability of the data are sought. Numerical solutions can be found using a variety of optimization procedures. In this study, we used MATLAB’s constrained nonlinear optimization function fmincon to minimize the negative log likelihood, which is equivalent to finding the maximum likelihood estimate (MLE). The MLE for the parameter values producing the predicted GF was the set of coefficients mapped to the pair of vectors \(}\) and \(}\) that produced the maximum likelihood value, i.e.,

$$\begin\begin\text\text\text=\underset}}}& \mathcal\left(},};}\right)\mapsto \end \left\,\widehat,\widehat\right\}\\ \end.$$

(11)

In summary, when magnitude GFs are the measurement of interest, the signal + noise distribution varies at SNRs below approximately 10 dB, violating the homoskedasticity assumption made by standard regression models. Inclusion of low-SNR data in the GF therefore necessitates using the Rician distribution to model the noise. Each element of the predicted GF was linked to its set of SNR-dependent Rician parameter values, and maximum likelihood estimation was used to find the predicted GF that best fit the observed GF data.

The structural equation approach described above yielded estimates for the latent variables associated with the non-linearity. Ideally, the values of these latent parameters (\(\widehat,\widehat,\widehat\)) are key points of interest for inferences and predictions about measures of clinical interest, such as cochlear nonlinearity, behavioral audiometric thresholds, loudness perception, and detection of subclinical damage. The strength of these relationships remains to be quantified using suitably large data sets, as well as compared with the performance of traditional descriptions of GFs, such as peak levels and slopes. To facilitate comparisons with traditional measures, each predicted GF was quantified at four points: 1) peak, 2) inflection point below the peak, 3) inflection point above the peak, and 4) maximum slope. We found these values by first fitting the predicted GF with a cubic spline (on a linear–linear scale), from which derivatives and their zero crossings could easily be computed.

Some previous work with GFs obtained from other paradigms has estimated a “knee point” or inflection point associated with the onset of cochlear compression (e.g., [40, 41]), defined as the point where the GF departs from approximate linear growth. (Note that some reports focus on GF compression in log–log space, as opposed to linear–linear space.) In the fixed L1 paradigm, the slope changes continuously throughout the GF (Fig. 1b); however, it does not change at a constant rate. We defined inflection points where the rate of change in slope was maximal (i.e., the location of the minimum and maximum of the second derivative of the spline fit in linear space). The inflection point below the peak was the L2 value at which the slope was most rapidly changing towards maximum compression. The inflection point above the peak was the L2 value at which the slope was most rapidly changing away from maximum compression.

In previous work, GFs have sometimes been used to estimate “thresholds”. Estimates have been made in log-log, semi-log, and linear-linear space, with thresholds defined as the lowest stimulus level at which DPOAEs are measurable or the level at which the DPOAEs theoretically become zero based on a mathematical function fit to the data (e.g., [42, 43]). For the fixed L1 paradigm, this latter concept of threshold does not appear to make sense, since a non-zero GF magnitude is expected for any f2 stimulus magnitude > 0 Pa, so long as there is some amount of non-linearity associated with cochlear vibration near the f2 place. Defining “threshold” as the L2 value where the GF reaches the noise floor is also undesirable, because such a measure conflates the signal and the noise levels. Noise floor varies as a function of frequency, hardware used, and averaging time. Noise floor also varies across individual subjects even when accounting for these other variables. In our sample of 12 participants, for example, noise floors showed a 10.7 dB range (Fig. 1a). Previous studies have rightly considered SNR when fitting GF magnitudes, for example including only data points with SNR > 6 dB ([e.g., 42, 43]). Such SNR-based criteria, like the Rician-based method described in this paper, improve the robustness of estimation procedures.

Based on theoretical expectations as well as human and animal measurements made in our labs, subjects with pristine hearing have steeper maximum GF slopes (on a linear-linear plot), while those with hearing loss have shallower slopes. In the fixed L1 paradigm, the steepest slope appears to be found in the limit as f2 stimulus magnitude (\(_|\)) approaches 0 Pa. Support for this assertion will be presented in a later section. We therefore defined GF maximum slope as the first derivative of the spline fit evaluated at \(_|\) = 0 Pa. Although this GF maximum slope is not a “threshold” in the traditional sense, it may correlate with behavioral audiometric thresholds in a manner similar to GF “thresholds” reported by previous work. Though it is not possible for a threshold measure to be truly independent of noise floor (since noise will always affect the precision of the DPOAE magnitude estimate), the slope of the predicted GF in the limit as \(_|\) approaches 0 Pa is otherwise independent of noise floor.

Confidence intervals (CIs) may be obtained in several ways. When evaluating the present GF data, we found that the commonly-used profile likelihood approach underestimated the actual variability, likely due to heteroskedastic residuals and relatively modest sample sizes. To obtain accurate interval estimates, we computed bootstrap confidence intervals. The paired bootstrap and residual bootstrap [44] are techniques commonly applied to regression models, but these are inappropriate in the present case due to nonconstant spacing on the x-axis and the presence of heteroskedasticity. In such cases, the wild bootstrap is an option [45, 46]. We implemented a modified version of the wild bootstrap which involved adding random samples of noise to the model fit, rather than adding reweighted versions of the residuals themselves. For each resample, \(}}_=}_+}_\), where the tilde accents indicate resampled values. \(}_\) was a complex number generated by the sum of two random draws from a normal distribution centered at zero with a standard deviation of \(\overline}_}}\sqrt\), with the second draw multiplied by \(\sqrt\) to make \(}_\) complex. This resampling was done for \(M\) sweeps, so the resulting matrix of resampled values could be processed as described by Eqs. 4–6. The resampled vectors \(\widetilde}}\) and \(\widetilde}}_}}}\) were fit as described above to obtain resampled MLE estimates \(\,\widetilde,\widetilde\}\). Using these resampled coefficients, predicted GF curves were generated, and their associated slopes, peaks, and inflection points were found. One hundred bootstrap resamples were obtained, and 95% confidence intervals (CIs) were computed from the sorted resamples for each measurement of interest. One hundred resamples is generally too small for determining statistical significance; however, in the present case the purpose of the resamples was to give a sense of the precision associated with the fits and measurements for each participant, for which 100 resamples was sufficient.