Perturbative QCD gives divergent series due to renormalons, and theoretical predictions using such series are essentially ambiguous. By subtracting renormalons from the Wilson coefficients in the framework of the operator product expansion (OPE), we can achieve a precise calculation of the QCD effect. We propose a method for renormalon subtraction with systematic approximation accuracy by using the Fourier transform. It utilizes the properties of the Fourier transform, and the Wilson coefficient with the renormalons removed is presented as a one-parameter integral whose integrand has suppressed (or vanished) renormalons. In this talk, I will show an application to B and D meson masses as one of the first analyses. By subtracting the ambiguities of $O(¥Lambda_{QCD})$ and $O(¥Lambda_{QCD}^2/m)$, the non-perturbative parameters of HQET are determined with high accuracy. The results are consistent with theoretical expectations, and improvements in convergence and scale dependence are confirmed.