A method of imaging at least two chemical species using magnetic resonance imaging with signal separation for two chemical species resulting in separate signal datasets for these two chemical species, the method comprising: - acquiring first and second echo datasets at different echo times resulting in first and second acquired complex datasets, - modelling the first and second acquired datasets by employing a multi-peak spectral signal model of at least one of the chemical species, said modelling resulting in first and second modelled complex datasets, said first and second modelled datasets comprising first and second phase errors and the separate signal datasets for the two chemical species and the modelling of the first and second datasets employing a linear combination of the separate signal datasets for the two chemical species, multiplied by first and second phasors, the first and second phasors comprising the first and second phase errors, wherein the weights for the linear combination are derived from the spectral signal model of the chemical species, wherein the determination of the separate signal datasets for the two chemical species comprises: - determining magnitudes of the echo data of the first and second acquired datasets by - deriving a bi-quadratic equation resulting from the system of two nonlinear equations obtained by the combination of the magnitudes of the echo data of the first and second acquired datasets with the first and second modelled datasets, - obtaining two pairs of associated true and false solutions for the separate signal datasets for each of the two chemical species by analytically solving the bi-quadratic equation, - obtaining true and false phasor candidates based on the difference between the phase errors of the first and second modelled complex datasets given the two pairs of associated true and false solutions for the separate signal datasets for each of the two chemical species, - extracting the true phasor from the true and false phasor candidates by performing a regional iterative phasor extraction and - determining the true solution for the separate signal for each of the two chemical species based on the extracted true phasor. - acquiring first and second echo datasets at different echo times resulting in first and second acquired complex datasets, - modelling the first and second acquired datasets by employing a multi-peak spectral signal model of at least one of the chemical species, said modelling resulting in first and second modelled complex datasets, said first and second modelled datasets comprising first and second phase errors and the separate signal datasets for the two chemical species and the modelling of the first and second datasets employing a linear combination of the separate signal datasets for the two chemical species, multiplied by first and second phasors, the first and second phasors comprising the first and second phase errors, wherein the weights for the linear combination are derived from the spectral signal model of the chemical species, wherein the determination of the separate signal datasets for the two chemical species comprises: - determining magnitudes of the echo data of the first and second acquired datasets by - deriving a bi-quadratic equation resulting from the system of two nonlinear equations obtained by the combination of the magnitudes of the echo data of the first and second acquired datasets with the first and second modelled datasets, - obtaining two pairs of associated true and false solutions for the separate signal datasets for each of the two chemical species by analytically solving the bi-quadratic equation, - obtaining true and false phasor candidates based on the difference between the phase errors of the first and second modelled complex datasets given the two pairs of associated true and false solutions for the separate signal datasets for each of the two chemical species, - extracting the true phasor from the true and false phasor candidates by performing a regional iterative phasor extraction and - determining the true solution for the separate signal for each of the two chemical species based on the extracted true phasor. - determining magnitudes of the echo data of the first and second acquired datasets by - deriving a bi-quadratic equation resulting from the system of two nonlinear equations obtained by the combination of the magnitudes of the echo data of the first and second acquired datasets with the first and second modelled datasets, - obtaining two pairs of associated true and false solutions for the separate signal datasets for each of the two chemical species by analytically solving the bi-quadratic equation, - obtaining true and false phasor candidates based on the difference between the phase errors of the first and second modelled complex datasets given the two pairs of associated true and false solutions for the separate signal datasets for each of the two chemical species, - extracting the true phasor from the true and false phasor candidates by performing a regional iterative phasor extraction and - determining the true solution for the separate signal for each of the two chemical species based on the extracted true phasor.