Ядро 2D Свертка сигнала в MATLAB

Я хочу реализовать двумерную свертку ядра сигнала в MATLAB. По сути, это нелинейный квадратичный фильтр сигнала (не изображения), описываемый следующей формулой:

где k2(t1,t2) — двумерное ядро ​​свертки. x представляет собой вектор сигнала 1d [N,1], а y представляет собой выходной сигнал [N,1].

До сих пор я делал это очень грубой силой, неэлегантным способом. Мне было интересно, могу ли я использовать функции MATLAB filter2/conv2, чтобы сделать это более эффективно! Я знаю, что эти функции предназначены для обработки изображений, о которых я не слишком много знаю, поэтому я надеялся, что кто-то может помочь!


matlab image-processing convolution signal-processing
person DankMasterDan    schedule 10.02.2014    source источник
comment
M и N одно и то же?   -  person Matt J    schedule 11.02.2014
comment
Нет, М - это размер ядра (памяти системы). N — длина сигнала. Примерные цифры, M = 30 ячеек, N = 10 000 ячеек.   -  person DankMasterDan    schedule 11.02.2014

Ответы (3)


arrow_upward
1
arrow_downward

Проблема при использовании diag(conv2(x,x,k)); заключается в том, что вы вычисляете нечто гораздо большее (всю 2d-матрицу), а затем сохраняете только диагональ. Это может быть дорого в зависимости от размера ваших сигналов. Вы можете попробовать с

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)

Он работает быстрее, чем другой вариант во многих случаях, которые я пробовал. Один вывод может, например,

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

ans =

   1.5059e-12

Я не знаю, насколько велики ваши сигналы или ваше ядро, поэтому вы можете попробовать.

person sebas    schedule 10.02.2014
comment
Примерные цифры, M = 30 ячеек, N = 10 000 ячеек. - person DankMasterDan; 11.02.2014
comment
@DankMasterDan Не пытайтесь использовать conv2 с такими размерами. Это, вероятно, никогда не закончится в обычном компьютере. - person sebas; 11.02.2014
comment
@DankMasterDan С небольшой модификацией он может быть более эффективным. Просто добавьте круглую скобку, чтобы умножить матричный вектор на вектор, который не транспонируется перед другим. x(imax:-1:imin)'*(k2(jmin:jmax,jmin:jmax)*x(imax:-1:imin)); Это всегда быстрее, чем без скобок. - person sebas; 11.02.2014
comment
Спасибо себас! Это не такое элегантное решение, как я себе представлял, но оно сократило мое вычислительное время вдвое! - person DankMasterDan; 11.02.2014
comment
Добро пожаловать. В следующий раз вы также можете пометить свой вопрос словом «алгоритм», и вы, вероятно, привлечете внимание более хороших программистов, чем я. - person sebas; 11.02.2014

arrow_upward
1
arrow_downward

Возможно следующим образом

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

or 

y=diag(conv2(x,x,k));
person Matt J    schedule 10.02.2014
comment
Настоящая оптимизация этого потребует знания относительных размеров x и k. - person Matt J; 11.02.2014

arrow_upward
0
arrow_downward

Для упомянутых вами размеров может быть целесообразно использовать БПФ, хотя это эквивалентно циклической свертке, поэтому вам нужно убедиться, что N включает достаточное заполнение нулями. Я предполагаю это ниже.

Из-за разделимости вам нужен только одномерный БПФ на x

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

  XX=X*X.';

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

Как указывает sebas, вы отбрасываете все, кроме диагонали ytmp, поэтому некоторые вычисления тратятся впустую. Если вы собираетесь повторять это для многих сигналов x длины N, вы можете сократить эти потери, предварительно вычислив матрицу IDFT,

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

Затем, вместо вычисления ytmp, вы можете вычислить диагональ ifft напрямую, выполнив

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

Даже если вы не повторяете вычисление для нескольких x, эффективность NlogN БПФ плюс отсутствие петель могут компенсировать потери.

person Matt J    schedule 11.02.2014
comment
вам нужно получить реальный (ytmp) после обратного fft? - person lennon310; 11.02.2014
comment
Я не знаю, следует ли ожидать реального результата. Это зависит от симметрии k, и ОП нам этого не сказал. Но если ожидается, что выходные данные будут действительными, то да, вам придется отбросить любой остаток с плавающей запятой в мнимой части, как вы это делаете с ifft() и др. - person Matt J; 11.02.2014