7
$\begingroup$

I have an $n \times n$ asymmetric convolution kernel, $k(t_1,t_2)$. $k$ is zero everywhere except for in small regions near the corners.

I also have an $n \times n$ image, $f$.

Let $*$ denote circular convolution so that the convolution operator, $K(f) = k*f$, returns an $n \times n$ image.

What is a formula for the transpose of $K$? That is, the operator $K^T$ such that $\langle K(f) ,g \rangle = \langle f, K^T(g) \rangle$.

I worked it out to be $K^T(f) = k(-t_1,-t_2)*f$ but that is not holding up to my checks.

%test how to make convolution transpose

clear all;close all;clc;

n = 1024;

%my kernel, it is complex.
[t1 t2] = meshgrid(1:n,1:n);
k = fft2(exp(-t1/50).*exp(-t2/75));

%convolution functions
convk = @(u) ifftn(fftn(k).*fftn(u));
rk = rot90(k);
convkt = @(u) ifftn(fftn(rk).*fftn(u));

%my images
f = rand(n);
g = rand(n);

trapz(trapz(convk(f).*g)) %<Kf,g>
trapz(trapz(f.*convkt(g))) %<f,Kg>, not the same as the above line...

edit: $$ \langle k*f,g \rangle = \int_{x,y}g(x,y) \int_{t_1,t_2} k(x - t_1, y-t_2)f(t_1,t_2) dt_1dt_2dxdy $$ $$ = \int_{t_1,t_2}f(t_1,t_2) \int_{x,y} k(x - t_1, y-t_2)g(x,y) dxdydt_1t_2 $$ $$ = \langle f, k(-1)*g \rangle $$

edit: Ok I got it now. There was an off-by-one error. You can see this by examining the (square) circular convolution matrix. This 1D example explains how to do it.

n = 256;

%convolution kernel
k = randn(n,1) + 1j*rand(n,1);

%circular convolution matrix
K = zeros(n);
for j=1:n
    K(:,j) = circshift(k,j-1);
end

%circular convolution fft-based function
convk = @(u) ifft(fft(k).*fft(u));

%the transpose of circular convolution, time-reversed complex-conjugated
%shifted-by-one kernel.
kt = circshift(flipud(k),1);
convkt = @(u) ifft(conj(fft(k)).*fft(u));

%test that they are the same thing
x = randi(100,[n 1]) + 1j*randn(n,1);
y = randn(n,1) + 1j*rand(n,1);
fprintf('The right matrix for K? : %d \n', norm(K*x - convk(x)));
fprintf('The right matrix for Kt? : %d \n', norm(K'*x - convkt(x)));
fprintf('The dot product test : %d \n', dot(convk(x),y) - dot(x,convkt(y)));
fprintf('The dot product test : %d \n', dot(x,convk(y)) - dot(convkt(x),y));

Will output:

The right matrix for K? : 5.246206e-11 
The right matrix for Kt? : 3.582131e-11 
The dot product test : 3.492460e-10 
The dot product test : 9.313226e-10

edit: For 2D, I think I've got it right but I'm seeing errors around $10^{-8}$. I want to just blame round-off since if anything changes in the below code I get gigantic errors.

clear all;close all;clc;

%nxn array
n = 256;

%make convolution kernel
[t1 t2] = meshgrid(1:n,1:n);
k = fft2(exp(-t1/50).*exp(-t2/75));

%circular convolution fft-based function
convk = @(u) ifftn(fftn(k).*fftn(u));

%the transpose of circular convolution, time-reversed complex-conjugated
%shifted-by-one kernel.
kt = rot90(k,2);
kt = conj(kt);
kt = circshift(kt,[1 1]);
convkt = @(u) ifftn(fftn(kt).*fftn(u));

%test that they are the same thing
x = rand(size(k)) + 1j*randn(size(k));
y = randn(size(k)) + 1j*rand(size(k));

%check the output error
fprintf('The dot product test : %d \n', sum(sum(conj(convk(x)).*y)) - sum(sum(conj(x).*convkt(y))) );

edit: Oops, I should be looking at relative error.

abs(sum(sum(conj(convk(x)).*y)) - sum(sum(conj(x).*convkt(y))))/sum(sum(conj(convk(x)).*y))

is close to machine precision.

Improve this question
$\endgroup$
8
  • $\begingroup$ Change $rk = rot90(k);$ to $rk = k.'$ instead. $\endgroup$ – Spacey May 12 '12 at 15:22
  • $\begingroup$ Ok, can you please show your work on how you arrived to the $k(-t1, -t2)*f$ conclusion? $\endgroup$ – Spacey May 13 '12 at 1:47
  • $\begingroup$ it's there now. It's a lazy check since I did the work over all of $\mathbb{R}^2$ and I'm actually interested in a torus. I don't see what will change though. $\endgroup$ – dranxo May 13 '12 at 2:45
  • $\begingroup$ Compton, I follow, except on the second line, I do not think it should be interpreted as a convolution anymore. I did a quick 1-D example on paper, and it looks to be more similar to a cross-correlation of a time-reversed $k$ with $g$, instead of a time-reversed convolution with $g$. (The minus became a plus). Look at first formula here. (en.wikipedia.org/wiki/Cross-correlation). I think this is what is going on. $\endgroup$ – Spacey May 13 '12 at 4:17
  • $\begingroup$ Yes. Cross-correlation is time-reversed convolution. I suppose the question could be "Is cross-correlation the transpose of convolution on a torus?" $\endgroup$ – dranxo May 13 '12 at 4:42
11
$\begingroup$

The transpose (or more generally, the Hermitian or conjugate transpose) of a filter is simply the matched filter. This is found by time reversing the kernel and taking the conjugate of all the values.

Chapter 5 of Jon F. Claerbout - Earth Soundings Analysis: Processing versus Inversion has a great discussion of the topic of Adjoint operators (or transpose, or back projection, or matched filtering). It's in the context of Seismic imaging, but it really gets to the heart of what you're actually trying to do (in general) with these techniques.

Improve this answer
$\endgroup$
6
  • $\begingroup$ Good point. I missed the complex conjugate. Still, adding that in won't fix the code in the question. Also, I didn't quite see in the article where they explain how the transpose is a matched filter. $\endgroup$ – dranxo May 13 '12 at 23:53
  • $\begingroup$ They didn't, it just is! Think about what the transpose actually defines and it should be apparent. $\endgroup$ – Henry Gomersall May 14 '12 at 21:25
  • $\begingroup$ Well yes, I have the definition up there in my question. The problem is that conjugate time reversal isn't doing what it should if it is indeed the transpose. $\endgroup$ – dranxo May 14 '12 at 23:12
  • $\begingroup$ Well, that's what the transpose is ;) What do you think it should be doing? $\endgroup$ – Henry Gomersall May 16 '12 at 9:00
  • $\begingroup$ I added an excellent reference to the reply. $\endgroup$ – Henry Gomersall May 16 '12 at 9:06
1
$\begingroup$

If you have convolution given by:

$$ \boldsymbol{y} = \boldsymbol{h} \ast \boldsymbol{x} $$

You can transform it into Matrix Form:

$$ \boldsymbol{y} = H \boldsymbol{x} $$

Then the operator $ {H}^{T} $ means doing the correlation:

$$ {H}^{T} \boldsymbol{x} = \boldsymbol{h} \star \boldsymbol{x} $$

In MATLAB Code, if you use the convolution as black box then Correlation is Convolution with the flipped signal - conv(vH(end:-1:1), x).
For 2D signals you need to flip them both vertically and horizontally - conv(mH(end:-1:1, end:-1:1), mX).

Related question - Solving Deconvolution using Conjugate Gradient.

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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