Metabolic profiling of natural products is used to map correlated concentration variances of known and unknown secondary metabolites in extracts. NMR-spectroscopy is in this respect regarded as convenient and reproducible technique with the ability to detect a wide range of small organic compounds. Two-dimensional J-resolved NMR-spectra are used in this context to resolve overlapping signals by separating the effect of J-coupling from the effect of chemical shifts. Often one-dimensional projections of these data are used as input for standard multivariate statistical methods and only the intensity variances along the chemical shift axis are taken into account. Here, we describe the use of parallel factor analysis (PARAFAC) as a tool to preprocess a set of two-dimensional J-resolved spectra with the aim of keeping the J-coupling information intact. PARAFAC is a mathematical decomposition method that fits three-way experimental data to a model whose parameters in this case reflect concentrations and individual components spectrum along the chemical shift axis and corresponding profiles along the J-coupling axis. A set of saffron samples, directly extracted with methanol-d4, were used as a model system to evaluate the feasibility and merits of the method. To successfully use PARAFAC the two-dimensional spectra (n = 96) had to be aligned and processed in narrow windows (0.04 ppm wide) along the chemical shift axis. Selection of windows and number of components for each PARAFAC-model was done automatically by evaluating amount of explained variance and core consistency values. Score plots showing the distribution of objects in relation to each other, and loading plots in the form of two-dimensional pseudo-spectra with the same appearance as the original J-resolved spectra but with positive and negative contributions are presented. Loadings are interpreted not only in terms of signals with different chemical shifts, but also the associated J-coupling profiles.