Kernel 2D Convolution of signal in MATLAB

Question

I want to implement 2d kernel convolution of a signal in MATLAB. This is essentially a nonlinear quadratic filter of a signal (not image) as described by the formula below:

where k2(t1,t2) is a 2-dimensional convolution kernel. x is a 1d [N,1] signal vector & y is the output [N,1] signal.

So far, I've been doing this by a very brute force, non-elegant fashion. I was wondering if I could use the MATLAB filter2/conv2 functions to do this more effectively! I know these functions are meant for image processing, which I don't know too much about, so I was hoping someone could help!

Answer 1

The problem when using diag(conv2(x,x,k)); is that you are calculating something much larger (the whole 2d matrix) and then you keep only the diagonal. It could be expensive depending on the size of your signals. You can try with

n = 500; m = 50;

x = rand(n,1);
k2 = rand(m,m);

tic; res1 = diag(conv2(x,x,k2)); toc;

tic;
res2 = zeros(n+m-1,1);
for k = 1:n+m-1
    imin = k+1-min(k,m); imax = min(k,n);
    jmin = max(1,k-n+1); jmax = min(k,m);
    res2(k) = x(imax:-1:imin)'*k2(jmin:jmax,jmin:jmax)*x(imax:-1:imin);
end
toc;
norm(res1-res2)

It works faster than the other option for many cases I have tried. One ouput can is, for instance

>> script
Elapsed time is 0.012753 seconds.
Elapsed time is 0.006541 seconds.

ans =

   1.5059e-12

I do not know how large are your signals or your kernel, so you can try.

Answer 2

Maybe as follows

y=diag(conv2(k,x*x.'));

or

y=diag(conv2(x,x,k));

Answer 3

For the sizes you mentioned, it could be worthwile using FFTs, although this is equivalent to cyclic convolution, so you would have to make sure N includes enough zero-padding. I assume this below.

Because of separability, you only need a 1D fft on x

  X=fft(x);
  K=fft2(k,N,N);

  XX=X*X.';

  ytmp=ifft2(XX.*K);
  y=diag(ytmp);

As sebas points out, you are discarding all but the diagonal of ytmp, so there is some wasted computation. If you are going to be repeating this for many signals x of length N, you can curb that waste by pre-computing a IDFT matrix,

  Afft=ifft2(eye(N)).'; 

Then, instead of computing ytmp, you can compute the diagonal of the ifft directly by doing

  y=sum(Afft2.*ifft( XX.*(K.')) );

Even if you aren't repeating the computation for multiple x, the NlogN efficiency of the FFT plus avoiding loops might still make up for the waste.