We present a new method for harmonic inversion in multi-dimensions, i.e., extracting the wave vectors omega(k) and amplitudes d(k) from a signal c(n)= Sigma(k)d(k)e(-in omega lambda), where n defines the multi-index. The method is an extension of the filter-diagonalization method for 1D signals. As such it enables the harmonic inversion in any small wavevector domain D-omega l by solving a small generalized eigenvalue problem. The computed omega(k) and d(k) can then be used to create a high resolution image F(omega) for omega is an element of D-omega. The method greatly overperforms the conventional Fourier analysis for a model 2D signal containing as many as 10 000 damped sinusoids with moderate amount of noise.