# Estimate Filter Coefficients from the Result of Linear Convolution with a Known Signal

If I have samples of input say x(1:500) and it passes through FIR filter with 9 taps and some unknown coefficients. The output y(1:508) is also known. The aim is to estimate the filter response in frequency domain. I know using Convolutional theorem $\hat{H} = (X^H*X)^{-1}(X^H*Y)$, but this requires fFT of full 500 samples, However I can only do fft of 128 samples at a time. How to estimate the filter response in frequency domain ?

• Couldn't you simply run a 128 sample impulse at the input, then do the FFT at the output? – a concerned citizen Dec 18 '17 at 7:30
• yep, assuming the taps aren't distributed over delays > 118, I don't see why you'd need all your input for classification. It would, however, be very interesting why you can't do a larger FFT. That is an unusual problem. – Marcus Müller Dec 18 '17 at 12:30
• @Munish, Welcome to our community. Could you please review my answer? Could you point if anything is missing? Thank You. – Royi Sep 19 '20 at 10:51

The idea here is to build the problem in its Matrix Form.
We have the filter $$h$$ which is represented by the matrix $$H$$ to represent Linear Convolution operation:

$$H = \begin{bmatrix} {h}_{1} & 0 & 0 & \ldots & & 0 \\ {h}_{2} & {h}_{1} & 0 & \ldots & & 0 \\ \vdots & & \ddots & & & \vdots \\ 0 & \ldots & {h}_{m} & {h}_{m - 1} & \ldots & {h}_{1} \end{bmatrix}$$

This means that the Linear Convolution is given by:

$$\boldsymbol{y} = H \boldsymbol{x}$$

This is the equivalent of vY = conv(vX, vH, 'full') in MATLAB.
It is clear that for $$\boldsymbol{x} \in \mathbb{R}^{n}$$ and $$\boldsymbol{h} \in \mathbb{R}^{m}$$ (Namely the filter has $$m$$ taps) then $$H \in \mathbb{R}^{ \left( n + m - 1 \right) \times n }$$.

The above matrix $$H$$ is Toeplitx Matrix.
You may be able to generate it using GenerateToeplitzConvMatrix() I added to my answer repository.

Bu we're after the coefficients of $$\boldsymbol{h}$$. The nice thing is we can write the above as:

$$\boldsymbol{y} = X \boldsymbol{h}$$

Where $$X \in \mathbb{R}^{\left( n + m - 1 \right) \times m}$$ is still a Toeplitz Matrix since the convolution (When shape is full in MATLAB notation) is commutative.

So once we have the matrix $$X$$ we can solve the above equation like any other linear system:

$$\hat{ \boldsymbol{h} } = \arg \min_{\boldsymbol{h}} {\left\| X \boldsymbol{h} - \boldsymbol{y} \right\|}_{2}^{2} = {\left( {X}^{T} X \right)}^{-1} {X}^{T} \boldsymbol{y}$$

But since $$X$$ has a special structure, being Toeplitz, there must be a way to exploit this, right?
Indeed there is, solving it using the Levinson Recursion. You'll be able to find a solver in the code repository of this question as well (Also see my Levinson Recursion Implementation). The idea is similar to Hilmar's solution.

But I had another idea. Toeplitz Matrix are a generalization of Circulant Matrices. If $$X$$ was a Circulant Matrix we could use the highly efficient fft() function to solve the problem.

So the question is how to implement Circular Convolution using Linear Convolution.
Namely, we want to adjust the matrix $$X$$ and the output $$\boldsymbol{y}$$ so the will represent Circular Convolution between $$\boldsymbol{x}$$ to $$\boldsymbol{h}$$.

Let's say we know how build those (Let's call them $$\bar{X}$$ and $$\bar{ \boldsymbol{y} }$$) then we have:

$$\bar{ \boldsymbol{y} } = \bar{X} \boldsymbol{h}$$

Yet since $$\bar{X}$$ it diagonalizable by the DFT Matrix then (See Eigen Values of Circulant Matrix):

$$\bar{X} = {F}^{H} \Omega F \Rightarrow \bar{ \boldsymbol{y} } = {F}^{H} \Omega F \boldsymbol{h} \Rightarrow \boldsymbol{h} = {F}^{H} {\Omega}^{-1} F \bar{\boldsymbol{y}}$$

Where $$F$$ is the DFT Matrix and $$\Omega$$ is a Diagonal Matrix of the Eigen Values of the matrix $$\bar{X}$$. Those Eigen Values are basically the DFT Transform of $$\boldsymbol{x}$$.

So basically we can solve this using the DFT (fft()). In MATLAB code it will be something in the lines of:

vH = ifft(fft(vYBar) ./ fft(vX), 'symmetric');


The full MATLAB code is available on my StackExchange Signal Processing Q45879 GitHub Repository.
It includes implementation of all 3 methods and one can see that each yields a perfect solution as long the system is well defined.

### Resources

You may wish to easily do sliding window fft, and to get a set of outputs. That means that for the first run you will do fft on samples 1:128, on the second run fft on samples 2:129 and so on. In this way, you will get 500-128 fft outputs, each one of 128 length. After this, you may average all fft sample outputs, and get what you are looking for. The second and a better way is possible if you have control over the filter, i.e, you can send input to your filter, the best way is to have just a single pulse running through it, and calculating fft after that.

If your system is noise free, you really only need the first nine samples of input and output.

$$h[0] = y[0]/x[0]$$ $$h[1] = (y[1]-h[0] \cdot x[1])/x[0]$$ $$...$$ If there is noise, numerical problems, or if the initial state of the filter is unknown, you can set this up as a least square error problem in the time domani.

Check LMS in Wikipedia, it is an iterative way, works for any size.