A data-driven approach is proposed to construct anisotropic damage models with a minimal number of internal variables from phase field - crack propagation simulations on Representative Volume Elements (RVEs) of quasi-brittle materials. The approach resorts in particular to a harmonic analysis of damage. The orientation distribution functions of two elastic moduli are numerically determined while accounting for the effects of the nucleation and propagation of microcracks by the phase field method. Given these two functions, the effective elastic tensor of a material without or with microcracks is uniquely determined. The expansions into two Fourier series of the relative variations of these two functions related to an undamaged reference state and to a damage state make appear damage internal variables naturally. The number and natures of these variables can be optimized by truncating the Fourier series according to the degree of approximation desired and using POD (Proper Orthogonal Decomposition). Thus, 2D and 3D anisotropic damage models can be constructed without resorting to usual assumptions made in damage mechanics. This construction holds for complex microstructures including image-based ones and for arbitrary loading history. 2D and 3D applications are provided to evaluate the accuracy of the damage models constructed and to show the potential of the approach proposed.