A high-order harmonic spectral filter (HSF) is implemented to smooth out a field variable defined on a spherical surface using the double Fourier series (DFS) as orthogonal basis functions. The filter consists of the solution of the high-order harmonic diffusion equation with the implicit method, where the high-order harmonic operator is split into second-or lower-order harmonic operators. The second-order harmonic operator is replaced by a pentadiagonal matrix whose elements are the spectral coefficients. First, a biharmonic spectral filter (BiHSF), the prototype of the high-order HSF, is developed where only the second-order harmonic operator is included. It is found that the computational error for the inversion of a pentadiagonal matrix remains in the order of machine rounding. Compared to the mixed-order HSF with DFS used in the previous study, which contains the second-and third-order harmonic operators, the BiHSF can provide a sharper cutoff of high wavenumbers as well as improved computational efficiency. These advantages come from the fact that for each set of spectral coefficients the BiHSF needs only a single inversion of the pentadiagonal matrix whereas the mixed-order HSF requires triple inversions and an auxiliary operation of the tridiagonal matrix. Based on the BiHSF, the highorder HSF up to the sixth order, which is composed of a multiple inversion of tri-or pentadiagonal matrices, is designed. Tests with the rotated Gaussian fields revealed that the HSF with DFS is isotropic. Application to various problems, including a time-dependent strongly nonlinear case and the observed flow, indicates that the high-order HSF is well capable of smoothing out selectively high-wavenumber components without significantly affecting the conserved quantity, such as total (kinetic) energy.