A method of image processing of magnetic resonance (MR) images for creating de-noised MR images, comprises the steps of providing image data sets including multiple complex MR images (S7), subjecting the MR images to a wavelet decomposition (S12) for creating coefficient data sets of wavelet coefficients (Sn,m) representing the MR images in a wavelet frequency domain, calculating normalized coefficient data sets of wavelet coefficients Formula (I) (S17), wherein the coefficient data sets are normalized with a quantitative amount of variation, in particular standard deviation Formula (II), of noise contributions included in the coefficient data sets (Sn,m), averaging the wavelet coefficients of each coefficient data set (S18) for providing averaged wavelet coefficients Formula (III) of the coefficient data sets, calculating phase difference maps (Δφn,m) for all coefficient data sets (S20), wherein the phase difference maps provide phase differences between the phase of each wavelet coefficient and the phase of the averaged wavelet coefficients Formula (III), calculating scaled averaged coefficient data sets of wavelet coefficients by scaling the averaged wavelet coefficients Formula (III) with scaling factors (Cn,m), which are obtained by comparing parts of the normalized wavelet coefficients of the normalized coefficient data sets Formula (I) that are in phase with the averaged wavelet coefficients Formula (III) (S22), calculating rescaled coefficient data sets of wavelet coefficients Formula (IV) (S24) by applying a transfer function Formula (V) on the coefficient data sets (Sn,m) and on the scaled averaged coefficient data sets, wherein the transfer function includes combined amplitude and phase filters, each depending on the normalized coefficient data sets Formula (I) and me phase difference maps (Δφn,m), resp., and subjecting the rescaled coefficient data sets to a wavelet reconstruction Formula (IV) (S25) for providing the denoised MR images.