Deconvolution and Polynomial factoring using the FFT

259 Views Asked by At

I've been trying to implement a general N dimensional deconvolver for various engineering applications and some math curiosities.

For speed and simplicity I've decided to try and do this with help of the fast fourier transform and the fact that for circular convolution ($*$), we have:

$$(f_1 * f_2)(t) = \mathcal{F}^{-1}(\mathcal{F}(f_1)\mathcal{F}(f_2))(t)$$

then assume we know $f_1*f_2$ and $f_2$ we try and find $f_1$ by a naiive divide strategy: $$ \mathcal{F}^{-1}\left(\frac{\mathcal{F}(f_1 * f_2)}{\mathcal{F}(f_2)+\epsilon}\right)$$

If we use circular convolution, how do we avoid effects like the one below:

$\left[\begin{array}{c} 1\\2\\1 \end{array}\right][1,2,1] = \left[\begin{array}{ccc}1&2&1\\2&4&2\\1&2&1\end{array}\right]$ divided/deconvolved by $[1,2,1]$ (should be $\left[\begin{array}{c} 1\\2\\1 \end{array}\right]$):

$$ans = \left[\begin{array}{cccccc} 5/6&1/6& -1/6& 1/6& -1/6& 1/6\\ 5/3& 1/3& -1/3& 1/3& -1/3& 1/3\\ 5/6& 1/6& -1/6& 1/6& -1/6& 1/6 \end{array}\right]$$

An "ordinary" convolution (with zero-padding) between $ans$ and [1,2,1] would not give us an correct answer: $$[1,2,1]*ans = \left[\begin{array}{cccccc} 2&1&0&0&0&0\\ 4&2&0&0&0&0\\ 2&1&0&0&0&0\\ \end{array}\right]$$

We can check that circular convolution between $ans$ and [1,2,1] would give us a shifted or rotated version of the 3x3 filter above ( and plenty of 0 entries ): $$[1,2,1](*_C)ans = \left[\begin{array}{cccccc} 2&1&0&0&0&1\\ 4&2&0&0&0&2\\ 2&1&0&0&0&1\\ \end{array}\right]$$

Is there an easy remedy to this? I have a few ideas for engineering solutions to this, but they are very ad-hoc.


Edit: Found one way that seemingly solves the problem above. There was one particular bug when deconvolving down to a singleton dimension, this seemed to be helped with the hint in the comments ( just to sum over the singleton dimension ).

Another example, to calculate the

$$f_1*f_2 = \left[\begin{array}{ccccc} 0&0&0&0&1\\ 0&0&0&0&0\\ 0&0&0&0&0\\ 0&0&0&0&0\\ -1&0&0&0&0 \end{array}\right] \hspace{0.5cm}f_2 = \left[\begin{array}{rr}0&1\\-1&0 \end{array}\right] \hspace{0.5cm} \text{gives:} \hspace{0.5cm} f_1 = \left[\begin{array}{rrrr}0&0&0&1\\0&0&1&0\\0&1&0&0\\1&0&0&0\end{array}\right]$$ Which corresponds to the $$(x^4-y^4) = (x-y)\sum_{i=0}^3x^iy^{k-i}$$

Here the second bug created two constant diagonal lines, but it was due to an error in how I handled the origo in the FFT-routines.