The multi-dimensional nature of neuroscience data has made the use of multi-way statistical analysis suitable in this field. Parallel Factor Analysis (PARAFAC) is a multidimensional generalization of PCA with the advantage of offering unique solutions. However, imposing physiologically acceptable constraints would improve the interpretation of this type of analysis. We propose a new algorithm called Alternating Penalized Least Squares to estimate PARAFAC solutions which allows the use of different kinds of constraints. Applied to semi-synthetic and real spontaneous EEG time-varying spectra, a wide range of sparse and smooth solutions can be found separately as well as with these two properties combined. Smoothness is usually desired in spectra and different sparse scenarios are observed in the temporal evolution of physiological intermittent phenomena. The degree of constraints can be tuned through the weighting parameters, whose optimal values can be chosen by cross-validation and other measures. Using physiologically founded constraints, PARAFAC analysis was able to recover the observed ‘waxing and waning’ temporal behavior and smoothness of the EEG spectra. This tool can also be applied to epilepsy, sleep data and other clinical applications in the recognition of neurophysiological patterns of interest.