electrophysiological brain signals are widely studied to get insights into the inner function of the brain .
the electro - encephalogram ( eeg ) , as an example , has been analyzed for decades and is particularly popular because of its noninvasiveness , wide availability , relatively small cost and excellent temporal resolution which enables capturing the fast neural dynamics . because it is known that electrophysiological brain signals exhibit important spectral characteristics , frequency transforms are often applied .
however , since in general brain signals do not possess the statistical property of stationarity , time - frequency transforms are of special interest .
such methods are able to represent a given signal jointly in the time - frequency domain , called its time - frequency representation ( tfr ) .
thereby the signal 's spectral components can be analyzed in relation to its temporal dynamics . in a wide sense
an important feature of biological signals , and particularly brain signals , is their inter- and intraindividual variability .
that is , under fixed experimental conditions , the obtained signals exhibit heterogeneousness not only between groups of subjects , but also between subjects within the same experimental group and even within the same subject between multiple experimental trials .
existing signal analysis techniques either do not take this issue into account adequately or usually treat it by defining frequency bands of interest rather than single frequencies , and similarly time intervals instead of sharp instants .
however , this strategy requires a priori knowledge about the variability to appropriately set the interval widths , and it is imprecise because it blindly includes all information contained in that time / frequency region .
a good example is the time - frequency coherence analysis of two given input signals , for instance , by means of the cross short - time fourier transform , the cross wigner - ville distribution , or the wavelet coherence [ 1 , 3 ] . all of these methods relate both signals at fixed time / frequency ( or scale ) locations .
thus , if signal a exhibits the same neural activation as signal b , but signal a 's pattern is shifted in frequency just a little , none of the abovementioned techniques will be able to find the strong similarity of a and b. although coherence estimation and other rigid strategies have been successfully applied  for more than 30 years 
, this issue has inspired us to develop a general flexible method of pattern analysis and corresponding feature extraction in electrophysiological tfr data . by abstracting from the tfr images and working with the representation of a tfr pattern ,
tfr patterns reduce dimensionality by representing neural activity in a wide spectrotemporal region by comparably few quantities .
pattern - based outlier detection has the potential to become a useful tool for data quality assurance . in ongoing studies
we are employing the presented method , for instance , to estimate functional brain connectivity by means of a pattern - based approach . in the following ,
we present the developed neuroinspired interpretable model which is able to capture general time - frequency patterns .
we use solely eeg data for demonstrations here , but our method is applicable to general electrophysiological signals or even to other signals showing similar behavior .  
algorithms for robustly fitting the novel model to time - frequency representations are presented . a strategy for finding an appropriate model order
is given , and distance functions are defined which quantify ( dis-)similarity of two given models .
these methods are tested in section 3 , where real as well as simulated data are used to demonstrate our technique 's functionality and robustness .
while our technique is not restricted to specific time - frequency distributions , we employ the smoothed pseudo wigner - ville distribution   in this work . this is a quadratic transform estimating signal power in the time - frequency domain , whereby all quantities in this work are real numbers .
the transform generates quite smooth tfrs , which means that neighboring pixel values are correlated .
although only positive values can be interpreted as signal power , the wigner - ville distribution introduces also negative values in general .
we will refer to time - frequency representations as mappings y : t  f   which estimate signal power for each point in the time - frequency domain .
of the numerous ways to quantify tfr patterns , we choose to fit a parametric surface to the data . because of the spatial correlation inherent in the data , traditional regression assumptions about independence of observations do not hold here .
this absence of strong gradients in the tfr images also invalidates most image feature extraction techniques , which are often based on edges and texture . our model ,
however , is especially designed for spatially correlated data ; furthermore its parameters are interpretable .
these quantities are useful features which embody important information about the underlying signal and thereby considerably reduce data dimensionality . using our method ,
we propose an extension of the well - known gaussian model for tfr analysis . the gaussian model for multivariate data x   is defined by 
 ( 1)y(g)(x)=c+a  exp(12(x)t1(x ) ) 
 with c   being a constant additive offset , a   the amplitude relative to the offset ,    the constant n - dimensional mean vector , and  denoting a n  n symmetric positive definite matrix .
positive definiteness ensures that the argument to the exponential function is always negative ; additionally we know that e is bounded by zero and one for negative x. therefore , the exponential factor scales the final amplitude between 0 and a relative to the offset c. the term ( x)(x   ) is also known as the squared mahalanobis distance of x with respect to  and  . because in our context this function represents arbitrary data in contrast to statistical distributions ,  will also be called the position vector , and  is the spread matrix , its entries ij are denoted spread parameters .
firstly , the model will be shaped like a peak for all possible parameter values by imposing the constraint that  ( and thus also  ) is symmetric positive definite .
thereby , the model will never be able to completely  degenerate .  because the model is not flexible enough to fit small local variations of an expected pattern , the gaussian model is relatively insensitive to local data outliers and is also unsusceptible to overfitting .
an additional aspect of robustness is that extreme peak deformations are directly reflected in extreme parameter values .
thereby degenerated models can be easily detected or may even be prevented by imposing parameter constraints .
the gaussian model is well - suited to extract bivariate peaks from brain signals ' tfr data , reflecting short intervals of neural excitement in a specific frequency range .
an instance of the above described surface , in the bivariate case , is fully identified by its parameter vector 
 ( 2)p=(c , a,1,2,11,12,22)t . 
 the absolute peak height c + a , the peak position ( 1 , 2 ) , and the peak orientation can be derived from the parameter vector .
further relevant quantities are the temporal peak onset , peak offset , and peak duration ( as the difference of the previous two ) .
as already mentioned , the gaussian model 's robustness comes at the cost of inflexibility .
while some local effects in tfr data can be appropriately explained by ( 1 ) , more general patterns of activation do not follow peak - like shapes , as will be shown later .
thus a generalization of the gaussian model would be desirable , especially concerning the ability to represent patterns of activation rather than just independent  events  in the spectrotemporal domain . at the same time , a generalized method should maintain maximum robustness in order not to degenerate easily and to prevent overfitting .
the so - called gaussian mixture modeling ( gmm ) is a straightforward extension , but this method still assumes ( multiple- ) peak - shaped data . in the following , the smooth natural gaussian extension ( snage ) model is presented as a flexible extension of the multivariate gaussian model . before giving a formal definition , we explain the idea in an intuitive way , guided by figure 1 . the n - variate gaussian model can be described by its ( n + 1)-dimensional peak point p = ( 1 ,  , n , a ) and its spread parameters controlling the exponential flattening relative to p along the independent variables ' dimensions . now the idea is to not use only one , but k peak points p   , i = 1  k , interpolated by a smooth ( n + 1)-dimensional curve of peaks .
a way to think of the surface in figure 1 modeled by snage is to  shape  it by sliding an n - variate gaussian model along the curve , its peak point being connected to the curve and thus varying in height ( a ) and position ( vector  ) .
thereby complex smoothly  bent  patterns of data with varying amplitude ( dependent variable ) can be captured .
analogously to the gaussian model , an n  n spread matrix determines the model 's shape , which is a surface for 
the tradeoff between robustness and flexibility can be controlled by choosing the number of peak points k. using many peak points will allow for good fits to complex patterns but will also increase the danger of overfitting .
choosing a small k yields a robust model , but its ability to capture complex patterns will be limited . by setting k = 1
, snage reduces to the traditional gaussian model as a special case . regarding the standard gaussian model , each function value y(x ) is fully determined by the mahalanobis distance of the point x to the unique mean point . but since the snage model offers infinitely many mean points
this issue will be addressed in the next section , where a formal definition is given .
let the  number of peak points  be denoted by k   , k  1 .
let  : [ 1 , k] denote a smooth n - dimensional curve of 
means ,  and let a : [ 1 , k] be a smooth one - dimensional curve of  amplitudes  along . let further the  offset  c   , and let    be a symmetric positive definite matrix
. define 
 ( 3)y~(snage):n[1,k]y~(snage)(x , u)c+a(u )  exp(12(x(u))t  1(x(u ) ) ) . 

note that y~(snage ) is a family of traditional gaussian models , parameterized by the curve parameter u. in order to construct a function which is independent of u , we define 
 ( 4)u(x)argmaxu[1,k]|y~(snage)(x , u)c| . 
 the snage model is then given by 
 ( 5)y(snage):ny(snage)(x)y~(snage)(x , u(x ) ) .
the function u*(x ) defines that traditional gaussian model which assigns to x the largest absolute amplitude relative to the offset c among all members of the gaussian family y~(snage ) .
this is necessary to cope with positive as well as negative a(u ) . in tfr analysis ,
a(u ) can be restricted to only positive values ( see section  2.3.1 ) , in which case ( 5 ) in fact simplifies to 
 ( 6)y(snage)+(x)maxu[1,k]y~(snage)(x , u ) . 
 for y as well as y the max (  ) function
curve , in the sense that ||(u1 )  (u2)|| is small , but |a(u1 )  a(u2)| is large .  
we assume that the diagonal entries of the spread matrix  , that is , the spread along each dimension , are sufficient to control a tfr pattern 's  width . 
therefore , we fix offdiagonal entries to zero for the sake of robustness . thereby , the mahalanobis distance in ( 3 ) reduces to the weighted euclidean distance .
the two curves a(u ) and (u ) are yet to be defined in terms of discrete parameter values . for good parameter interpretability
, we choose to form both curves by interpolating k points p   , i = 1  k by b - splines of degree 3 , which yields the curve : 
 ( 7)p:[1,k]n+1 ,  p(u)=(p1(u)  pn(u)pn+1(u))=(1(u) 
, we obtain a sufficiently smooth model which inherits both the splines ' flexibility and the gaussian standard model 's robustness .
our model further inherits the b - splines ' local control property ; that is , varying a p will affect the model only in the p 's vicinity .
additionally , the degree of flexibility can be adapted to the data at hand by varying k from 1 ( single gaussian peak ) to arbitrary flexibility with k > 1 .  
an instance of the n - variate snage model with k points ( or of order k ) is fully represented by the parameter vector of length 1 + n + ( n + 1)k : 
 ( 8)p=(c,11,22,  ,nn , p1(1), 
while the offset c and the spread parameters are mainly responsible for the prediction of data values , the curve p(u ) interpolating the p is directly interpretable as it models the main path of peaks in the data . in the next section
we show how to fit the model to data in a robust manner . given a time - frequency representation y(x ) , x  t  f , we aim to find a parameter vector p so that the respective model y fits the data  best ,  in the sense that it minimizes a cost function .
we use the sum of squared differences of the data and the modeled surface : 
 ( 9)sse(p)=x=(t , f)ttt , ff(y(tfr)(x)y(snage)(x))2 . 

this quantity is also called the  sum of squares due to error  which is zero if the model perfectly fits the data .
sse is a nonlinear function dependent on p. using squared differences pronounces outliers , but these are not expected to occur frequently in our smooth tfr data . in order to find a locally minimum solution of sse , a nonlinear least squares algorithm implemented in the matlab optimization toolbox , lsqnonlin ,
given an initial parameter vector and the cost function sse , this optimizer produces a sequence of models pi .
the iteration hopefully converges to a p * with minimum cost , that is , best resemblance between model and data .
we attempt to provide advantageous starting conditions for the optimizer by preprocessing the tfr and by obtaining an initial parameter vector p0 which is expected to be close to the optimum with respect to sse .
moreover , an iterative refinement scheme is proposed to be able to robustly fit models of high order .
tfr data are badly scaled , showing differences of several orders of magnitude in values of time , frequency , and signal power , which affects optimization performance . to address this problem
, lsqnonlin offers a way to take into account typical values for each dimension for gradient estimation . also concerning this issue , any euclidean distance operating on tfr data in our algorithms
is weighted appropriately . furthermore , smooth objective functions are desirable so that the low - order taylor approximations used during optimization resemble the cost function in a relatively large neighborhood around the current point . to this end ,
the tfr images are smoothed and subsampled , which has the additional benefit of faster cost function evaluations . finally , since negative values in the time - frequency domain are not interpretable , they are usually set to zero .
since lsqnonlin performs local optimization , a start parameter vector should be chosen such that it already lies in the basin of the cost function 's ( unknown ) global minimum , that is , a vector whose cost is already low . to this end ,
the constant offset and the spread parameters are initialized to c = 0 , 11 = ( tmax  tmin)/5 and 22 = ( fmax  fmin)/5 if the underlying tfr 's time and frequency axes are bounded by tmin , tmax and fmin , fmax , respectively .
one may choose any sensible values alike , but the spreads should not be initialized too small in order to obtain a generalizable model . yet , more thought has to be put in choosing the number and coordinates of the peak points p. in the following we propose a way to compute initial estimates of the time , and frequency coordinates of all p directly from the data
. a reasonable approximation to the unknown optimal curve of peaks can be found by tracing a path through the tfr image y from left to right , which runs through areas of high pixel intensity .
more precisely , the sum of all pixel intensities along this path should be as high as possible .
this is an optimization problem in turn , yet its solution can be computed in quadratic time complexity ( provided that the path 's slope is bounded ) by a dynamic programming algorithm , see . because a global optimum is guaranteed to be found , this strategy is insensitive to local outliers and noise . to this end , a similar approach as described in   is employed .
following the notation therein , we define our energy function to be equal to the tfr values themselves , that is , e(i ) = i. a horizontal path s * ( called seam in ) which maximizes ( this is in contrast to , where minimum energy seams are computed ) this simple cost function is found by dynamic programming . see figure 3 for an example .
additionally , the paths are constrained by imposing an upper bound k on their slopes .
this value depends on the time - frequency resolution here and once more represents a compromise between robustness and flexibility .
given the found path s * and the desired model order , k evenly spaced samples are subsequently drawn from a smooth approximating curve to obtain estimates of the first two coordinates of the p , i = 1  k. we choose to empirically set the ps ' last components , interpretable as amplitudes relative to the initial constant offset c = 0 , to max(y )  c = max(y ) .
once a parametric representation of the data is available , its accuracy can be improved in a step - wise manner , as is presented in the following section .
as already stated , the number of points p controls the model 's robustness which complements its ability to resemble complex patterns .
therefore , the demand for a near - optimal initial parameter vector increases with the model order k. employing the optimal image path method described in the previous section yields a 
reasonable  estimation , but sampling k equidistant points p ,  i = 1  k from the resulting curve is a simplification .
in fact , it can be observed that the optimizer tends to concentrate the p in time - frequency regions of high signal variability . for low model order
k , this shortcoming of our initial parameter estimation algorithm can be compensated easily by the optimization algorithm , but it may become a problem for increasingly flexible models . for this reason we propose an iterative scheme.find initial parameters
p0 by means of an optimal path ( see section  2.3.2 ) for a first , robust model of low order k0 .
let i0.fit the model to the data to obtain optimal parameters pi*.construct the optimal curve of peaks p(u ) by interpolation of the peak points ( see section  2.2.1).obtain the ki+1 = ki + 1 peak points for a refined model by uniformly sampling ( with respect to the spline 's sites uj ) the curve computed in the previous step.construct the parameter vector pi+1 of the refined model from pi by replacing the old ki peak points with the ki+1 new ones.if ki
< kmax , let ii + 1 and continue with step 2 . find initial parameters p0 by means of an optimal path ( see section  2.3.2 ) for a first , robust model of low order k0 .
fit the model to the data to obtain optimal parameters pi*. construct the optimal curve of peaks p(u ) by interpolation of the peak points ( see section  2.2.1 ) . obtain the ki+1 = ki + 1 peak points for a refined model by uniformly sampling ( with respect to the spline 's sites uj ) the curve computed in the previous step . construct the parameter vector pi+1 of the refined model from pi by replacing the old ki peak points with the ki+1 new ones .
the curve found in steps 2 and 3 will exhibit smaller gradient magnitude , that is , traversal speed , in areas of high signal variability than in other regions .
we aim at maintaining the curve - defining points ' optimum distribution found by the fitting algorithm and at enhancing the model 's flexibility mainly in these areas .
it turns out that simply by uniformly sampling the fitted curve ( step 4 ) we obtain a new interpolated curve which retains these properties .
an application of this algorithm is demonstrated in section  3.1 , where the resulting sequence of nested models is evaluated . in this section
distance measures are necessary , for instance , to quantify how well the data exhibit an expected pattern .
distance functions which are based solely on the curve of peaks p(u ) were found to be quite effective .
other possibilities include parameter vector distances and pixel - wise differences of signal power of the models ' generated data . by comparison ,
curve - based distance functions have the advantage of being able to interrelate models of different orders .
additionally , they are not influenced by the less informative parameters ( offset c and the entries of  ) . a popular distance measure for parametric curves
is the frchet distance . in the continuous case , the frchet distance of two parametric curves p1(u ) and p2(u )
is defined by 
 ( 10)df_max(p1,p2)=inf,  maxu  d  ( p1((u)),p2((u ) ) ) . 

here , (t ) and (t ) are monotone reparameterizations of the two curves , and d (  ) denotes ( weighted ) euclidean distance . in words ,
we search for those reparameterizations which make the curves the most similar with respect to maximum point - wise euclidean distance along the curves .
this maximum for these reparameterizations is returned as the two curves ' continuous frchet distance . in practice ,
the discrete frchet distance is frequently applied , whose computation is based on dynamic programming once more . in the discrete case , an additional distance function df_sum(p1 , p2 ) can be obtained by replacing the max function with a sum over u. that way , df_sum represents an average distance , being less prone to outliers in the curves .
we demonstrate the workflow to determine the appropriate model order k by fitting a tfr of real eeg data in this section .
typically we determine the necessary model complexity by fitting data with good signal to noise ratio ( snr ) in order to prevent the overestimation of k. for example , one possibility to achieve sufficient data quality is to average several tfrs which are expected to show similar patterns . the averaged tfr of real eeg data shown in figure 5 will guide the following explanations .
the depicted brain signals located in the lower frequency bands were recorded from the temporal brain region during a face recognition experiment .
these data exhibit a pattern of activity which is too complex to be captured by a traditional gaussian model . the iterative scheme described in section  2.3.3
an optimal path ( see section  2.3.2 ) estimates the initial parameters for the first , least flexible model of order k0 .
since the appropriate k is still unknown , a sufficiently large number kmax = 7 is chosen for the refinement . at each refinement stage the respective model
is evaluated , and in the end the most suitable 
 ( 11)k{k0,  ,kmax } ,  k0kkmax 
 is chosen as the model order for future fittings on lower - quality data .
these are the cost function value ( sse , see section  2.3 ) , the coefficient of determination r and its adjusted version radj .
although the use of quantities based on the coefficient of determination is discouraged for nonlinear models , they are applied here nonetheless for two reasons .
they are found to perform well for our purposes , and the proposed alternatives ( aic   and bic ) are not easily applicable here .
this is because the assumption of normally distributed residuals often does not hold , which is supported by a highly significant shapiro - wilk test   at  = 5% yielding p < 10 for this experiment .
the results show that for the data at hand a model of order k * = 4 is sufficient to capture the variability .
having determined the maximum model complexity on high - quality data , such a model can now be fitted to the rest of the data .
if the tfrs are not expected to vary substantially , like when fitting a model to signals from several nearby sensors , a previous fit may serve as the initial model .
however , if , for instance , multiple data segments of the same sensor should be fitted , the tfrs ' patterns may vary strongly . in this case , initial parameters should be chosen depending on the data by using the method of optimal paths described in section  2.3.2 .
since in this example k * = 4 is a quite moderate number , the iterative refinement may also be skipped . however , in general we would start with k0 = 3 or k0 = 4 and refine up to the determined k * , as proposed in section  2.3.3 .
we want to assess our model 's robustness by simulating data and measuring how strongly the model is affected by additive gaussian noise . to this end
, artificial data are created in the time domain , and their tfrs are computed to which our model will be fitted .  
we created a signal consisting of three consecutive oscillations , representing an alpha - theta - alpha eeg pattern at 10  hz/4  hz/10  hz respectively over a time span of 2.5 seconds .
the simulated sample rate is 250  hz . a plot is shown in figure 9 .
these data are quite challenging for our model because three distinct peaks emerge in the tfr which could be more appropriately modeled by a mixture of independent gaussian peaks .
however , we want to demonstrate the flexibility of the snage model which should also be able to cope with patterns of this form . for the following experiments we chose to start the fitting with a model of order k = 5 to account for the pattern 's complexity and perform one refinement step .
initial parameters are estimated by finding optimal paths , which means that no a - priori information about the known optimal model is passed to the fitting procedure other than the number of p to use .
we define the optimal model by fitting the noise - free simulation in the same way .
the distance measures from section  2.4 are used to determine how well the simulated pattern is found . in this experiment we added gaussian noise , which is appropriately filtered with respect to the sampling frequency , to the simulated data in the time domain .
signals exhibiting signal to noise ratios of 15  db up to + 10  db were generated in steps of 2.5  db . at each snr , ten distinct noise realizations are created to obtain representative results .
this independent noise in the time domain will produce correlated noise in the time - frequency domain due to smoothing .
this experiment serves to assess how strongly our algorithm is affected by pattern variability , respectively , to investigate its robustness .
small pattern distortions should ideally only slightly alter the optimal model , reflecting its robustness and avoidance of overfitting .
we further want to find out to what degree our model is able to find the simulated pattern at all .
we note here that adding noise increases the tfrs ' maximum amplitudes exponentially which strongly affects the comparability of different models . without normalization
but this would merely reflect the decreasing data amplitudes and contain no information about the quality of fit .
normalizing maximum data values are not an appropriate option either , because , for negative snrs , this would keep the noise constant while exponentially shrinking the pattern 's pixel intensities . even if the optimal model was perfectly recovered from the noisy simulation , high distances
only if signal power is excluded from model distance estimation , the returned values are useful representatives of how well the pattern was found .
the snage 's robustness to noise with respect to the pattern 's power is therefore not regarded here .
this is done by setting the third dimension of the path of peaks p(u ) to zero during frchet distance computation . 

convergence to the optimal model seems to require high signal to noise ratios . at 7.5  db ,
the distance measures ' variances fall off , reflecting the point of reliable pattern extraction .
apparently , the noise and interferences introduced in this experiment considerably impair the fitting process . in figure 11 , this issue is exemplarily investigated . at the positive snr of 5  db , where distance variances across the noise realizations are still high
the pattern was in fact found , but only in a different way than was expected .
nevertheless , this example shows that the impact different kinds of noise may have on the fitting process . to get a better feel for the average ability to fit the pattern under the influence of noise , see figure 12 . at each noise level ,
shown is a sequence of mean models which progressively look more similar to the true pattern .
in fact , the average fitting capability concerning both the positioning of the peak points in the time - frequency domain and the estimation of surface values is better than expected after having studied figure 10 . apparently , although the mean distances are still decreasing at negative signal to noise ratios , they are already small enough for successful pattern extraction on average .
an example is the subplot corresponding to snr = 0  db in figure 12 , which already clearly resembles the simulated pattern.this experiment shows that interferences between the desired signal and additive noise affect the fitting process quite strongly in the worst case .
positive signal to noise ratios of at least 7.5  db are found to be necessary for reliable pattern extraction in this investigation
. however , successful data modeling is also possible at lower snrs , as is seen in the average case .
the snage model is especially suited for time - frequency representations of electrophysiological signals because of their ( expected ) nonnegativity , smoothness and their patterns following a path of peaks .
however , our robust model is able to cope with data which do not exactly meet these requirements .
in order to retain robustness , we imposed several restrictions on our model , like neglecting offdiagonal spread parameters and holding the spread matrix constant over the curve of peaks .
an interesting question remains how the model 's flexibility and robustness would be affected if these constraints were dropped .
little effort would be necessary to include the stated extensions . as is typical for nonlinear optimization problems ,
therefore , a - priori knowledge about the optimal model can be incorporated by starting the optimization with a model which was previously fitted to similar data . moreover
, an algorithm based on an optimal path is developed to estimate initial model parameters directly from the data .
however , this method is limited to positive peak polarity by trying to maximize the path 's average amplitude .
additionally , the technique will not be able to find initial models which exhibit multiple contemporary components .
in such a case , the nonlinear optimization algorithm , which was found to work well , must compensate .
further strategies for the estimation of initial parameters would be desirable . in particular , the extraction of the curve interpolation points p from the optimal path possesses potential for improvement .
an open question is how we should deal with the spatial correlation of both the dependent variable and the residuals in a statistical inferential context .
further work is necessary to facilitate statistical testing , for instance , to assess the null hypothesis that an expected pattern is not contained in the data .
concerning the presented measures of model distance , the frchet distances were found to be very useful to assess model similarity in our experiments involving simulated noise .
their distinct advantage is their independence of model order and the disregard of the less interpretable parameters . on the other hand , spurious high distances
this can be attributed to the fact that the simulated data exhibit three independent peaks , which is a violation of the snage 's assumption of a connected path of peaks .
therefore , a combination of frchet values and a pixel - wise distance function based on the models ' generated data seems advantageous .
when applied to noisy time signals , the snage model adapts too well to the corresponding smooth time - frequency representations .
because the optimized cost function does not take into account information about the expected pattern , the model simply tries to capture the tfr data as accurately as possible .
data preprocessing and tfr interference suppression are therefore extremely important . adding penalty terms to the cost function and/or providing explicit initial parameters are ways to point the optimizer in the right direction .
however , even without specifying a - priori knowledge , the model was able to find the simulated pattern for low signal to noise ratios in the average case .
the analysis of time - frequency representations of electrophysiological signals calls for flexible methods accounting for inter- and intraindividual data variability .
we present the flexible , robust , and interpretable model snage , which extends the established gaussian model .
its ability to extract 3d features from time - frequency representations of electrophysiological data is demonstrated
. however , the model applies to general multivariate data which exhibit similar behavior . in this work , several techniques to improve the model fitting performance are described .
we show how to estimate start parameters directly from the data . an iterative scheme to refine optimized models
experiments with real as well as simulated data demonstrate the snage model 's robustness and flexibility . under the influence of severe noise ,
the developed technique is best suited for patterns which are too complex to be appropriately captured by a gaussian model , but still simple enough to facilitate robust fits . to summarize , due to its robustness and flexibility the snage model possesses the potential to become a beneficial tool for practical eeg / meg analysis , including functional brain connectivity analysis , outlier detection , time - frequency denoising , and feature extraction .
 electrophysiological signals such as the eeg , meg , or lfps have been extensively studied over the last decades , and elaborate signal processing algorithms have been developed for their analysis . 
 many of these methods are based on time - frequency decomposition to account for the signals ' spectral properties while maintaining their temporal dynamics . 
 however , the data typically exhibit intra- and interindividual variability . 
 existing algorithms often do not take into account this variability , for instance by using fixed frequency bands . 
 this shortcoming has inspired us to develop a new robust and flexible method for time - frequency analysis and signal feature extraction using the novel smooth natural gaussian extension ( snage ) model . 
 the model is nonlinear , and its parameters are interpretable . 
 we propose an algorithm to derive initial parameters based on dynamic programming for nonlinear fitting and describe an iterative refinement scheme to robustly fit high - order models . 
 we further present distance functions to be able to compare different instances of our model . 
 the method 's functionality and robustness are demonstrated using simulated as well as real data . 
 the snage model is a general tool allowing for a wide range of applications in biomedical data analysis . 
1. Introduction
2. Methods
3. Results
4. Discussion
5. Conclusion