A Gauss quadrature method in which the wavelet is used as a weighting function is developed for wavelet BEM. Non-orthogonal spline wavelets that can change the order of vanishing moments as well as the order of polynomial are considered in BE analysis. Although the increase in the order of vanishing moments leads to the increase in the sparseness of matrices, that also increases the number of intervals in which the wavelet is described by a certain polynomial. The proposed quadrature method does not need to divide the support of wavelets in the calculation of matrix coefficients, while the Gauss-Legendre formula obliges us to divide the support into several intervals. Consequently the proposed method allows to reduce the computational work for generation of matrices. Estimation of the error in the numerical integration is also attempted in order to decide the number of integral points for a specified tolerance. Numerical experiments axe carried out, and the feasibility of the proposed method is examined.