Фильтрация в частотной области

Ваш текущий способ заполнения фильтра, чтобы он был того же размера, что и изображение, в основном правильный. Мы часто небрежно говорим о фильтрации сигнала длины M с помощью фильтра длины 3, но неявное предположение состоит в том, что мы дополняем оба сигнала до длины M или, возможно, до длины M+3-1.

Некоторые детали вашего подхода усложняют ситуацию:

1) Умножение на (-1) ^ (x + y) просто переводит ДПФ и не требуется. (См. Основы обработки сигналов Таблицу 3.7 "Круговой сдвиг частоты" для одномерного случая. В этом обозначение, вы позволяете k_0 быть N/2, поэтому член W_N в левом столбце просто переключается между -1 и 1.)

2) Поскольку фильтр Превитта имеет ненулевую поддержку только 3x3, ваш вывод должен иметь размер N+2 на N+2. Здесь нужно запомнить формулу длина (сигнал) + длина (фильтр) - 1.

Вот как бы я подошел к этому:

clear
x = im2double(imread('cameraman.tif'));
[M, N] = size(x);

h = [-1 0 1;
    -1 0 1;
    -1 0 1];

P = M + size(h,1) - 1;
Q = N + size(h,2) - 1;

xPadded = x;
xPadded(P, Q) = 0;

hPadded = h;
hPadded(P,Q) = 0;

hShifted = circshift(hPadded, [-1 -1]);

H = fft2(hShifted);
X = fft2(xPadded);

Y = H .* X;
y = ifft2(Y);

yCropped = y(1:M, 1:N);
imshow(yCropped,[]);

Вот как я решил свою проблему. Сначала я удалил шаги 2 и 7 из алгоритма. Затем отцентрируйте преобразование, поменяв местами первую половину индексов со второй половиной как в горизонтальном, так и в вертикальном направлении. Я сделал это, чтобы центрировать преобразование изображения. Затем я отменил это после вычисления обратного ДПФ результирующей матрицы. Я не уверен, почему мой описанный выше метод не работает, но теперь он работает.

1) Преобразуйте матрицу изображения NxN в матрицу 2*Nx2*N, добавив нули

2) Вычислить ДПФ матрицы изображения

3) Центрируйте преобразование изображения, поменяв местами первую и вторую половину строк и столбцов.

4) Создайте фильтр размерами 2Nx2N и центром в координатах (N,N)

5) Умножить матрицу изображения на матрицу фильтра

6) Рассчитайте обратное ДПФ и извлеките действительную часть результата.

7) Децентрализуйте результат, повторно применив шаг 3 к результирующей матрице.

8) Наконец, извлеките верхнюю левую часть NxN результирующей матрицы.

Выше приведена модифицированная версия шагов, которые я выполнил при применении моей фильтрации. Вот мой код (отредактированная/новая версия)

function res=myFreqConv(signal,mask)
    signal=double(signal);

    l=size(signal,1);

    % padding the image matrix with zeros and making it's size equal to
    % 2Nx2N
    pad_signal=zeros(2*l,2*l);
    pad_signal(1:l,1:l)=signal;
    m=size(mask,1);
    mask_f=zeros(2*l,2*l);

    % Creating a mask of 2Nx2N dims where the prewitt filter values are  
    at 
    % the center of the mask i.e. the indices are like this 
    % [(N-1,N-1), (N-1,N), (N-1,N+1);(N,N-1), (N,N), (N,N+1); (N+1,N-1), 
    (N+1,N), (N+1,N+1)] 
    for i=-1:1
       mask_f(l+i,l-1)=-1;
       mask_f(l+i,l+1)=1;
    end

    % calculate DFT of mask
    mask_f=myDFT(mask_f);

    signal_dft=myDFT(pad_signal);

    % shifting the image transform to center
    indices=cell(1,2);
    indices{1}=[2*l/2+1:2*l 1:2*l/2];
    indices{2}=[2*l/2+1:2*l 1:2*l/2];
    signal_dft=signal_dft(indices{:});

    %multiply mask with image
    res=mask_f.*signal_dft;
    res=real(myIDFT(res));

    % shifting the image transform back to original
    res=res(indices{:});
    res=res(1:l,1:l);
end