We present a very fast method for calculating an inverse filter for audio reproduction system. The proposed method of FFT-based inverse filter design, which combines the well-known principles of least squares optimization and regularization, can be used for inverting systems comprising any number of inputs and outputs. The method was developed for the purpose of designing digital filters for multi-channel sound reproduction. It is typically several hundred times faster than a conventional steepest descent algorithm implemented in the time domain. A matrix of causal inverse FIR (finite impulse response) filters is calculated by optimizing the performance of the filters at a large number of discrete frequencies. Consequently, this deconvolution method is useful only when it is feasible in practice to use relatively long inverse filters. The circular convolution effect in the time domain is controlled by zeroth-order regularization of the inversion problem. It is necessary to set the regularization parameter β to an appropriate value, but the exact value of β is usually not critical. For single-channel systems, a reliable numerical method for determining β without the need for subjective assessment is given. The deconvolution method is based on the analysis of a matrix of exact least squares inverse filters. The positions of the poles of those filters are shown to be particularly important.