2
$\begingroup$

Suppose that I am given two $n \times m$ matrices $\bf{A}$ and $\bf{C}$, and let $\bf{B}$ be a matrix that is convolved with $\bf{A}$, such that:

$\bf{A} * B = C$

In the above, $*$ is the convolution operator (link).

I am searching for a numerical procedure (and ideally a reference book or paper on the mathematics of the operation) that will allow me to determine $\bf{B}$, given $\bf{A}$ and $\bf{C}$.

I know both $\bf{A}$ and $\bf{C}$ exactly. The $\bf{A}$ is the system input, whereas the $\bf{C}$ is the system output.

Although the convolution is a product in the Fourier domain, I might be able to determine $\bf{B}$ by division. However, using the 2D Fast Fourier Transform (FFT), do I have to pad the matrices in the spatial domain or deal with division by zero? How do I do this in an elegent manner?

Is there a way to do this in Matlab using the fft2 and ifft2 functions?

3 Answers 3

3

In your case, there is no noise in the output, then the Fourier method is a safe way for deconvolution: first zero-pad $\mathbf{A}$ to the size of $\mathbf{C}$ (without zero-padding will leads to result corresponding to circulant convolution); then compute the spectra using fft2; divide the spectrum of $\mathbf{C}$ by the one of $\mathbf{A}$, then you can get the deconvolution result by using ifft2.

However, if the output contains noise, then the Fourier method may not a proper way for deconvolution, since division in frequency domain could amplify the noise and contaminate the results. The better way is, as Seyhmus Güngören mentioned, to convert the deconvolution problem into a system of linear equations or optimization problem with certain regularizations, then use iterative algorithms to solve the equations or optimization problem.

  • 0
    Whoops, this is actually my mistake; it was a bug in my test code! When calculating $\bf{C} ./ \bf{A}$, where the division is point-by-point (as in Matlab), and $\bf{C}$ and $\bf{A}$ are the same padded size, and $\bf{C}$ and $\bf{A}$ are complex, the deconvolution by division works well, provided that elements of $\bf{A}$ are not zero. In Matlab code, `ifft2(fft2(C,nn,mm) ./ fft2(A,nn,mm))`, where `nn` and `mm` are the padded size of $\bf{C}$, which is size $(n_A+n_B-1) \times (m_A+m_B-1)$. Thank you for all of your help, chaohuang.2012-09-04
2

Yes you need to make zero padding if you intend to do everything in the DFT domain. However there is an obvious risk of division with zeros. There are iterative algorithms to get rid of division with zero. These algorithms are not specific for the applications in the DFT domain.

I am not completely sure but I think it might be possible to convert the problem to a $1D$ convolution, reshaping the matrices as $1D$ vectors. Then one can use MATLAB with deconv command.

Another way is to write codes for your implementation. Basically you have a linear system of equations with some unknowns. If you can write them in a suitable form what remains is to solve the system. It is simply the multiplication of the inverse of the matrix at the left side with the matrix at the right side. But one should construct these matrices before.

  • 0
    @SeyhmusGüngören: Thanks for pointing this out; nothing is better than an implementation.2012-09-03
1

An extremely-readable paper on the mathematics of the deconvolution problem [1] discusses deconvolution for 1D vectors and 2D matrices. The paper is written from the perspective of an applied mathematician, and recasts the deconvolution problem as a special case of Fredholm integral equations. It discusses periodic and nonperiodic deconvolution problems, outlines a framework for deconvolution without using the FFT, and is probably one of the most exciting papers that I have ever read.

[1] P. C. Hansen, “Deconvolution and Regularization with Toeplitz Matrices,” Numerical Algorithms, vol. 29, no. 4, pp. 323–378, 2002.

  • 0
    Blind deconvolution is useful when the kernel $\bf{B}$ is not known exactly, and $\bf{A}$ needs to be recovered from $\bf{C}$.2012-09-03