4
$\begingroup$

As a continuation of the question Least Squares Solution Using the DFT vs Wiener Hopf Equations raised by Dan Boschen.

The question is, given the model:

$$ \boldsymbol{y} = \boldsymbol{h} * \boldsymbol{x} + \boldsymbol{n} $$

Where $ \boldsymbol{y} $ and $ \boldsymbol{x} $ are known discrete signals (Here as a vectors) and $ \boldsymbol{n} $ is additive white noise.

We're after the Least Squares Estimation of $ \boldsymbol{h} $ under the following 2 convolution models:

  1. The $ * $ operator is the discrete convolution with zero boundary conditions.
    Also known as full convolution.
  2. The $ * $ operator is the discrete convolution with cyclic boundary conditions.
    Also known as periodic / cyclic convolution.

For reference, for the full convolution, the convolution matrix can be generated by (Code by Dan):

def convmtx(h,n):
    # creates the convolution (Toeplitz) matrix
    # h is input array, 
    # n is length of array to convolve 
    return linalg.toeplitz(np.hstack([h, np.zeros(n-1)]), np.hstack([h[0], np.zeros(n-1)]))
$\endgroup$
1

1 Answer 1

2
$\begingroup$

The solution to the models will be different as the model will be different.

We can rewrite the problem in a matrix form:

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

Where the matrix $ X $ is different between the 2 models:

  1. For the full case the convolution matrix is a Toeplitz Matrix.
  2. For the cyclic case the convolution matrix is a Circulant Matrix (A particular kind of Toeplitz matrix).

Remark: Pay attention that $ \boldsymbol{y} $ is also different for each case, so the real problem becomes:

$$ \boldsymbol{y}_{\text{Type}} = {X}_{\text{Type}} \boldsymbol{x} + \boldsymbol{n} $$

Where $ \text{type} \in \left\{ \text{Full}, \text{Cyclic} \right\} $.

To show the difference between the 2 models we'll build the matrix per model.
We'll use $ \boldsymbol{x} = {\left[ 1, 2, \ldots, n \right]}^{T} $ where $ n = 5 $, then the matrices are given by:

enter image description here

We can see the model matrix is different, hence the output is different (Not only in value but in size).
Also pay attention that there is some overlap. This is known that subset of the solution of the cyclic convolution matches the classic convolution.

Least Squares Solution

In Least Squares Solution Using the DFT vs Wiener Hopf Equations, my answer shows a derivation of the LS solution for the Cyclic case.

In this section we'll do the solution without noise when we try the 4 combinations for the values of:

  • $ \boldsymbol{x} = {\left[ 1, 2, 3, 4, 5, 6, 7 \right]}^{T} $.
  • $ \boldsymbol{h} = {\left[ 1, 2, 3, 4, 5 \right]}^{T} $.
    The vector $ \boldsymbol{h} $ is normalized to have sum of 1.
Data / Solution Full Cyclic
Full [0.0667, 0.1333, 0.2000, 0.2666, 0.3333]
Cyclic [0.0667, 0.1333, 0.2000, 0.2666, 0.3333]

How can we solve for the cases off diagonal?
Straight out from the code the dimensions won't match!

When we do simulations, we usually know the model hence we employ the method which matches our model.

In real world, choosing the model is one the crucial hyper parameters of solution.

In image processing, the assumption on the boundaries, by itself, makes a lot of difference.

But for 1D signals, we usually do a simple assumption. We assume we see the steady state of the convolution, which is the same for all boundary conditions.

For the case we choose the length of the filter, and it happens to be much smaller then the number of the sample we have, we can average many samples in order to have a better SNR for the estimation (Under the assumption of white noise). Yet in the case above, the number of samples is low, hence the number of samples with valid convolution is smaller than the number of coefficients, hence without a proper boundary conditions, the problem can not be solved.

So let's look at the case the length of the data is much longer than the length of the filter:

enter image description here

Now we can see that rows 5-17 in the cyclic case matches 5-17 on the full convolution. There are enough samples to estimate the coefficients. The more we have, the less important the boundary conditions become.

Remarks

  • Numerically it is better to solve the problem as stated in the Least squares problem. In the code of Dan he calculates the correlation matrix and try to solve the normal equations. This is bad numerically as it effectively doubles the condition number of the system.
  • It also better to solve using specialized solvers. For instance, do not use the inverse explicitly. For Circulant matrices use the DFT, for general Toeplitz, use Toeplitz Solver (Levinson Recursion, See my Levinson Recursion solver at my Projects page to see how faster it can be).
$\endgroup$
4
  • $\begingroup$ This is great Royi, and very nearly and answer to the referenced question I posted as well. What I don't understand, if it makes sense to comment here, is that the test case I showed (dsp.stackexchange.com/a/87352/21048) meets your description of having a data length much longer than the filter length, so as you explained here the effect of the boundary conditions should be negligible- yet it was the "numerically bad" approach that resulted in the better answer. $\endgroup$ Apr 6, 2023 at 18:46
  • $\begingroup$ This post as well as the ending remarks in particular are very interesting and helpful to me -- rather than long discussion here I will (once I come back up for air) I will post some additional supporting questions focused on that and this problem. Thank you! $\endgroup$ Apr 6, 2023 at 18:47
  • $\begingroup$ @DanBoschen, I didn't say the approach is bad, I said that the way it was implemented is less stable numerically wise. As creating the R matrix means you double the condition number of the system. Having large number of samples doesn't mean it is a free launch. You need to select properly the sub set of samples. It means also you need to adapt the usage of the DFT() and IDFT(). Not like suggested there. $\endgroup$
    – Royi
    Apr 6, 2023 at 19:48
  • $\begingroup$ Yes no problem - more numerically stable is opposite of bad to me, I have a lot to learn from you in that space and need some more time first to dig in and not waste anyone's time. I have my suspicions on the additional intuition I was looking for thanks to yours and the other responses/ discussions-but I need to confirm that, as well as understand better your points on best approach to even do the initial W-H approach I take and how I think I can make the FFT approach work-it may take me a while as I have to switch gears but I will post more related questions once I can make it clean. $\endgroup$ Apr 7, 2023 at 0:11

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.