Re: Signal: Least-squares deconvolution. Autoregressive model, impulse responses, yule-walker, etc.

Hi Gideon.

Thanks for the answer :)

As far as I can tell, the code you sent solves exactly for an impulse
response of the same length as x and y. Then discards the unwanted samples
of the learned IR afterwards. Solving directly can be preferable in some
cases for performance, but it solves a slightly different problem. It does
not guarantee the least squares solution.

To get the least squares solution, one assumes that the impulse response is
shorter than the signal from the start. Then solves the normal equation
instead of directly solving the convolution linear system.

It's the setting up of the normal equation that can be a bit cumbersome. If
it's useful enough to justify a scipy function you guys can decide :)

Here is the code + some  the math:
Code <https://gist.github.com/maka89/eec636033251d1466570f9989ca14a59>
Jupyter Notebook
<https://gist.github.com/maka89/a94f9c135374557746148dfae324074f>


Best regards,
Aleksander Kringstad



On Sun, 30 Jul 2023 at 12:23, Gideon Genadi Kogan <[email protected]> wrote:

> Dear Aleksander,
>
> I think that this functionality is already implemented in
> *linalg.solve_toeplitz*.
> See this <https://dsp.stackexchange.com/a/88853/34391> code for usage
> example.
>
> Best,
> Gideon
>
> On Sat, Jul 29, 2023 at 11:06 AM Aleksander Kringstad <[email protected]>
> wrote:
>
>>
>> Hi,
>> I have written some efficient code for learning impulse responses using
>> least-squares minimization.  First time trying to contribute. Wanted to
>> know if it seems interesting and if someone wants to review.
>>
>> *Problem*
>> I.e. given input x and output y (both N samples long), we wish to learn
>> an impulse response with M < N samples that minimizes the least squares
>> error in y.
>>
>> Given
>> yp[n] = x[n]h[0]+x[n-1]h[1]+x[n-2]h[2]+... + x[n-M]h[M],
>>
>> minimize sum( (y-yp)**2)  wrt. h[0],h[1],h[2],... ,
>>
>> Circular convolution is assumed.
>>
>> *Applications*
>> - Learning impulse responses.
>> - Learn parameters of autoregressive models.
>>
>> If you use x and y so that x[n] = y[n-1], t*he algorithm will be equal
>> to learning the parameters of an autoregressive model using the Yule-Walker
>> equation*.
>>
>> *Performance*
>> It will be fast for long time-series, since the A matrix and the b vector
>> in the normal equation are calculated using FFTs, instead of the regular
>> dot(X.T,X) and dot(X.T,b).
>>
>> Also, the toeplitz structure of the normal equation matrix is exploited
>> for performance. This allows for learning longer impulse responses. And for
>> using less memory to store the normal equation matrix.
>>
>> *Notes*
>>
>> Circular convolution is often not what we want, but can be overcome by
>> zero-padding.
>>
>> Alternatively, one wants to avoid that the first M samples of the output
>> to contributes to the error at all.This can be achieved by first using the
>> algorithm to learn the impulse response and then fine-tuning it using
>> scipy.optimize
>>
>> Can add functionality for this fine-tuning + zero-padding. Should also
>> add a few alternatives for linear solvers.
>>
>> Core functionality is only 10 lines of self-written code and uses only
>> numpy and scipy.linalg.solve_toeplitz.
>>
>>
>>
>>
>>
>> _______________________________________________
>> SciPy-Dev mailing list -- [email protected]
>> To unsubscribe send an email to [email protected]
>> https://mail.python.org/mailman3/lists/scipy-dev.python.org/
>> Member address: [email protected]
>>
> _______________________________________________
> SciPy-Dev mailing list -- [email protected]
> To unsubscribe send an email to [email protected]
> https://mail.python.org/mailman3/lists/scipy-dev.python.org/
> Member address: [email protected]
>

_______________________________________________
SciPy-Dev mailing list -- [email protected]
To unsubscribe send an email to [email protected]
https://mail.python.org/mailman3/lists/scipy-dev.python.org/
Member address: [email protected]