Abstract: Sparse variational approximations are popular methods for scaling up inference in Gaussian processes to larger datasets. For $N$ training points, exact inference has $O(N^3)$ cost; with $M << N$ features, sparse variational methods have $O(NM^2)$ cost. Recently, methods have been proposed using harmonic features; when the domain is spherical, the resultant method has $O(M^3)$ cost, but in the common case of a Euclidean domain, previous methods do not avoid the $O(N)$ scaling and are generally limited to a fairly small class of kernels. In this work, we propose integrated Fourier features, with which we can obtain $O(M^3)$ cost, and the method can easily be applied to any covariance function for which we can easily evaluate the spectral density. We provide convergence results, and synthetic experiments showing practical performance gains.