Estimate Filter Coefficients from the Result of Linear Convolution with a Known Signal
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} &amp; 0 &amp; 0 &amp; \ldots &amp;  &amp; 0 \\ 
{h}_{2} &amp; {h}_{1} &amp; 0 &amp; \ldots &amp;  &amp; 0 \\ 
\vdots &amp;  &amp; \ddots &amp;  &amp;  &amp; \vdots \\ 
0 &amp; \ldots &amp; {h}_{m} &amp; {h}_{m - 1} &amp; \ldots &amp; {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


Circulant Preconditioned Toeplitz Least Squares Iterations.
Conjugate Gradient Methods for Toeplitz Systems.
Scientific Applications of Iterative Toeplitz Solvers.
Iterative Methods for Linear Systems with Matrix Structure.
Least Squares and Total Least Squares Methods in Image Restoration.
Fast Iterative Methods for Least Squares Estimations.
A Weakly Stable Algorithm for General Toeplitz Systems.
Algorithms for Toeplitz Matrices with Applications to Image Deblurring.
A Superfast Toeplitz Solver with Improved Numerical Stability.
Super Fast Solution of Real Positive Definite Toeplitz Systems.
Fast Solution of Toeplitz and Cauchy Like Least Squares Problems (Code in Giuseppe Rodriguez Page - TPLS).
A Super Fast and Stable Toeplitz Direct Solver via Randomization and Hierarchically Semi Separable (HSS) Methods (Code in Jianlin Xia Page - HSS Software Package).