Convolution where each filter is applied with matrix multiplication?

In my case it would be 2D convolution. So for example if the filter is 5x5 it would be applied to a 5x5 patch of the input and the two matrices would be multiplied using ordinary matrix multiplication to produce a 5x5 output matrix. Is there a nice way to implement this?

I think .unfold could be helpful.
Have a look at this post and let me know if you could adapt it to your use case.