# FIR filter design using DFT matrix inversion

In my previous question I use a system identification case as an example. Here I came across another question. I want to design a FIR filter of length $$L$$ to equalize an LTI system whose impulse response is $$h(n)$$ of length $$N$$, and my target response is $$d(n)$$ which is maybe a delayed $$\delta(n-D)$$. The problem can be formulated as finding the unknown $$g(n)$$ that satisfies $$h(n)*g(n)=d(n)$$ Frequency domain deconvolution can be used. Calculate $$N+L-1$$ point FFT $$H(k) \cdot G(k) = D(k)$$ Impose regularization we have $$G(k) = \frac{D(k)H(k)^*}{H(k)H(k)^*+\epsilon(k)}$$

Now I have to transform $$G(k)$$ back to time domain by $$N+L-1$$ point IFFT. But I can't guarantee that the IFFT result will end up with $$L-1$$ point zeros. I can simply truncate the result as @Hilmar suggests, but I don't think this is mathematically correct.

I came up with a solution with the DFT matrix. Let $$\mathbf{g}= \big[g(0), g(1), \ldots, g(N-1)\big]^T$$ be the time-domain vector of $$g(n)$$, and $$\mathbf{G} = \big[G(0), G(1), \ldots, G(N+L-2)\big]^T$$ be the frequency-domain vector of $$G(k)$$, the DFT matrix $$T$$ represents the linear transform $$\mathbf{Tg=G}$$ However, due to $$g(n)$$ and $$G(k)$$ have different lengths, this DFT matrix is not a square matrix, it is a truncated DFT matrix.

Finally, the FIR coefficients can be derived by the least-squares method $$\mathbf{g}=\big(\mathbf{T}^H\mathbf{T}+\lambda \mathbf{I}\big)^{-1}\mathbf{T}^H \mathbf{G}$$

Since the final solution is a least-square solution, I can't guarantee that $$N+L-1$$ point FFT of $$g(n)$$ equals to $$G(k)$$.

Is it a good way to design FIR filters? Does it have any disadvantage? Any better way?

• Inverting an LTI system is a bit of an art. The best method will depend A LOT on the properties of your system (minimum phase, bandpass vs full range, latency), the noise conditions (SNR as a function of frequency, time domain artifacts, clipping, ...) and the quantitative requirements of your application. It would help if you could narrow it down a bit? As it stands, the questions feels too broad. This will work in some cases and totally bomb in others Feb 8 at 10:53
• @Hilmar My goal is to design FIR filters to achieve some acoustic tasks, such as SISO task impulse response shortening, MIMO task like crosstalk cancellation and sound field reproduction. Most of them can be defined as least-squares problems in the frequency domain and should be transformed back to a time-domain FIR filter. Feb 8 at 11:09
• @Hilmar Using LTI inversion as an example, my system is the linear part of the electroacoustic transfer path from a loudspeaker to a microphone in a reverberant environment. The IR is measured with exponential swept sine. I want to reshape the impulse response of the whole system to improve its transient response, so the target response is either an unit delta function or the minimum-phase part of the original IR. Feb 8 at 11:16
• However, I would like to get an answer that is not limited to this example. Because designing filters in the frequency domain is often very simple, especially in MIMO systems, the final response of the system can be expressed as a multiplication of matrices and vectors. The frequency domain approach is studied in many papers, but the authors avoid to tell how it is eventually converted to a time-domain FIR filter. Feb 8 at 11:23
• Well, depending on your frequency range of interest, the bandwidth of your transducers and the reverb time of your environment, you are going to need A LOT of FIR samples. Also a room impulse response is by nature IIR (exponential decays) so your initial assumption of "Length N" doesn't really hold. Feb 9 at 8:19

Given that we have information already in the time domain, do the least squares solution in that domain, which is the approach used for least squares equalizers: This works quite well to correct for a linear channel distortion when the sounding signal uses has rich spectral occupancy (since the accuracy of the solution for any given point in frequency will be proportional to the SNR for that point). I suspect (but haven't confirmed) that equalizing the channel in the frequency domain using the DFT would suffer the same drawbacks as using the "Frequency Sampling Method" for FIR filter design, which is inferior due to time domain aliasing as demonstrated at this link.

The graphic below shows the solution we seek: we know the sounding sequence $$t[n]$$ and the received output of the channel $$r[n]$$. We desire a corrected $$r[n]$$ to match as close as possible the original sounding sequence, using the equalizer filter with coefficients $$h[n]$$. The corrected sequence at the output of the equalizer is shown as $$t'[n]$$, which we desire $$t'[n]=t[n]$$ after allowable delay. The output of the equalizer FIR is the convolution of $$r[n]$$ with the equalizer coefficients $$h[n]$$. Convolution can be written as a matrix equation by using $$r[n]$$ in a Toeplitz matrix as demonstrated in the graphic below. Thus the resulting matrix equation is given as:

$$t' = A \cdot h$$

We wish to solve for $$h$$, but $$A$$ is not a square matrix so can't be inverted. We can multiply $$A$$ by its conjugate transpose to make it invertible, multiplying both sides to not change the equation, resulting in the classical form for the least squares solution to an overdetermined equation:

$$A^Tt' = A^TA \cdot h$$

$$h = (A^TA)^{-1}A^T t'$$

I further detail this including MATLAB code and how intentional delay is importantly included at this link.

The same formulation can instead be used to estimate the channel rather than invert it as demonstrated in this link.

• This is a fully time-domain method using convolution matrix. Any difference to the frequency-domain method? BTW, is it LS Equalizer or LMS Equalizer in the graph? Feb 8 at 11:36
• Wouldn't the frequency domain method would suffer the same degradations as designing filters using the "frequency sampling method": due to time domain aliasing many more coefficients would be needed to properly equalize the channel to the desired accuracy. As the time span increases, the noise also increases so we wish to minimize the time span relative to the actual impulse response time of the channel. Feb 8 at 11:41
• @DSPnovice I believe LS and LMS are the same: An actual LMS equalizer uses an adaptive algorithm that iteratively converges to the least squares solution (used for changing channel conditions)---what is shown here is the least squares solution that such an algorithm would converge to but computed directly (when computing resources allow). Feb 8 at 11:44
• Do you mean that the time domain LS method is better than the frequency domain method at least for channel equalization? Does my method called "frequency sampling method"? Indeed another method use very long FFT point to calculate the FIR filter, and truncate it to a shorter length, which I think is called "windowing method". However in my method no truncation is performed. Feb 8 at 11:59
• @DSPnovice Why don't you try it both ways and compare the two approaches (under lower SNR conditions) both in performance and processing required? I think that would be a very interesting result...sometimes "better" doesn't justify the cost. Although at first glance here I don't see any difference in processing required, do you? Actually perhaps there is more work to do it in the frequency domain as you need to take DFT's to get to the same starting point? Feb 8 at 12:11