An efficient scheme is proposed for calculating the elements of the overlap and Hamiltonian matrices in filter diagonalization. This method is based on the cosine Fourier transform between the angle and order domains of the Chebyshev operator. The fast Fourier transform (FFT) implementation of this method renders a favorable scaling law in calculating the matrix elements. Our method represents a generalization of existing schemes and is applicable to matrices of any functions of the Hamiltonian operator. (C) 1997 Elsevier Science B.V.