We consider the problem of computing a $k$-sparse approximation to the Fourier transform of a length $N$ signal. Our main result is a randomized algorithm for computing such an approximation (i.e. achieving the $\ell_2/\ell_2$ sparse recovery guarantees using Fourier measurements) using $O_d(k\log N\log\log N)$ samples of the signal in time domain that runs in time $O_d(k\log^{d+3} N)$, where $d\geq 1$ is the dimensionality of the Fourier transform. The sample complexity matches the lower bound of $Ω(k\log (N/k))$ for non-adaptive algorithms due to \cite{DIPW} for any $k\leq N^{1-δ}$ for a constant $δ>0$ up to an $O(\log\log N)$ factor. Prior to our work a result with comparable sample complexity $k\log N \log^{O(1)}\log N$ and sublinear runtime was known for the Fourier transform on the line \cite{IKP}, but for any dimension $d\geq 2$ previously known techniques either suffered from a polylogarithmic factor loss in sample complexity or required $Ω(N)$ runtime.