  Properly separating and subtracting renormalons in the framework of the
operator product expansion (OPE) is a way to realize high precision computation
of QCD effects in high energy physics. We propose a new method (FTRS method),
which enables to subtract multiple renormalons simultaneously from a general
observable. It utilizes a property of Fourier transform, and the leading Wilson
coefficient is written in a one-parameter integral form whose integrand has
suppressed (or vanishing) renormalons. The renormalon subtraction scheme
coincides with the usual principal-value prescription at large orders. We
perform test analyses and subtract the ${\cal O}(\Lambda_{\rm QCD}^4)$
renormalon from the Adler function, the ${\cal O}(\Lambda_{\rm QCD}^2)$
renormalon from the $B\to X_ul\bar{\nu}$ decay width, and the ${\cal
O}(\Lambda_{\rm QCD})$ and ${\cal O}(\Lambda_{\rm QCD}^2)$ renormalons from the
$B,\,D$ meson masses. The analyses show good consistency with theoretical
expectations, such as improved convergence and scale dependence. In particular
we obtain $\bar{\Lambda}_{\rm FTRS}=0.495\pm0.053~\text{GeV}$ and$
(\mu_\pi^2)_{\rm FTRS}=-0.12\pm 0.23~\text{GeV}^2$ for the non-perturbative
parameters of HQET. We explain the formulation and analyses in detail.
