Matlab: 2D дискретное преобразование Фурье и обратное - PullRequest
Новая
       30

Matlab: 2D дискретное преобразование Фурье и обратное

2 голосов
/ 28 февраля 2020

Я пытаюсь запустить программу в matlab для получения прямого и обратного ДПФ для изображения в оттенках серого, но я не могу восстановить исходное изображение после применения обратного. Я получаю комплексные числа в качестве моего обратного вывода. Как будто я теряю информацию. Есть идеи по этому поводу? Вот мой код: 

%2D discrete Fourier transform
%Image Dimension

M=3;
N=3;
f=zeros(M,N);
f(2,1:3)=1;
f(3,1:3)=0.5;
f(1,2)=0.5;
f(3,2)=1;
f(2,2)=0;

figure;imshow(f,[0 1],'InitialMagnification','fit')


%Direct transform


for u=0:1:M-1
   for v=0:1:N-1
       for x=1:1:M
           for y=1:1:N

             F(u+1,v+1)=f(x,y)*exp(-2*pi*(1i)*((u*(x-1)/M)+(v*(y-1)/N)));


            end
        end
    end
end


Fab=abs(F);

figure;imshow(Fab,[0 1],'InitialMagnification','fit')



%Inverse Transform

for x=0:1:M-1
    for y=0:1:N-1
       for u=1:1:M
            for v=1:1:N

                z(x+1,y+1)=(1/M*N)*F(u,v)*exp(2*pi*(1i)*(((u-1)*x/M)+((v-1)*y/N)));


            end
        end
    end
end

figure;imshow(real(z),[0 1],'InitialMagnification','fit')

Ответы [ 2 ]

1 голос
/ 28 февраля 2020

Во-первых, самая большая ошибка в том, что вы неправильно вычисляете преобразование Фурье. При вычислении F вы должны суммировать по x и y, чего вы не делаете. Вот как это исправить: 

F = zeros(M, N);

for u=0:1:M-1
   for v=0:1:N-1
       for x=1:1:M
           for y=1:1:N
             F(u+1,v+1)=F(u+1,v+1) + f(x,y)*exp(-2*pi*(1i)*((u*(x-1)/M)+(v*(y-1)/N)));
            end
        end
    end
end

Во-вторых, в обратном преобразовании ваш брекетинг неверен. Это должно быть 1/(M*N), а не (1/M*N).

Кроме того, за счет немного большего объема памяти вы можете ускорить вычисления, не вкладывая так много циклов. А именно, при вычислении БПФ вместо этого выполните следующие действия: 

x = (1:1:M)'; % x is a column vector
y = (1:1:N) ; % y is a row vector

for u = 0:1:M-1
    for v = 0:1:N-1
        F2(u+1,v+1) = sum(f .* exp(-2i * pi * (u*(x-1)/M + v*(y-1)/N)), 'all');
    end
end

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

x = (1:1:M)'; % x is in dimension 1
y = (1:1:N) ; % y is in dimension 2
u = permute(0:1:M-1, [1, 3, 2]); % x-freqs in dimension 3
v = permute(0:1:N-1, [1, 4, 3, 2]); % y-freqs in dimension 4
% sum the exponential terms in x and y, which are in dimensions 1 and 2. 
% If you are using r2018a or older, the below summation should be 
%    sum(sum(..., 1), 2) 
% instead of 
%    sum(..., [1,2])
F3 = sum(f .* exp(-2i * pi * (u.*(x-1)/M + v.*(y-1)/N)), [1, 2]); 
% The resulting array F3 is 1 x 1 x M x N, to make it M x N, simply shiftdim or squeeze
F3 = squeeze(F3);
1 голос
/ 28 февраля 2020

Есть несколько проблем с вашим кодом:

  1. Вы неправильно применяете определение DFT (или IDFT): вам нужно sum над исходной переменной (переменными), чтобы получить преобразование. Смотрите формулу здесь ; обратите внимание на сумму.

  2. В IDFT константа нормализации должна быть 1/(M*N) (не 1/M*N).

Обратите также внимание, что код можно сделать более компактным путем векторизации, избегая циклов; или просто с помощью функций fft2 и ifft2. Я предполагаю, что вы хотите вычислить его вручную и «низкоуровневый», чтобы проверить результаты.

Код с двумя исправлениями выглядит следующим образом. Изменения помечены комментариями. 

M=3;
N=3;
f=zeros(M,N);
f(2,1:3)=1;
f(3,1:3)=0.5;
f(1,2)=0.5;
f(3,2)=1;
f(2,2)=0;

figure;imshow(f,[0 1],'InitialMagnification','fit')

%Direct transform

F = zeros(M,N); % initiallize to 0
for u=0:1:M-1
   for v=0:1:N-1
       for x=1:1:M
           for y=1:1:N
               F(u+1,v+1) = F(u+1,v+1) + ...
                   f(x,y)*exp(-2*pi*(1i)*((u*(x-1)/M)+(v*(y-1)/N))); % add term
            end
        end
    end
end

Fab=abs(F);
figure;imshow(Fab,[0 1],'InitialMagnification','fit')

%Inverse Transform

z = zeros(M,N);
for x=0:1:M-1
    for y=0:1:N-1
       for u=1:1:M
            for v=1:1:N
                z(x+1,y+1) = z(x+1,y+1) + (1/(M*N)) * ... % corrected scale factor
                    F(u,v)*exp(2*pi*(1i)*(((u-1)*x/M)+((v-1)*y/N))); % add term
            end
        end
    end
end

figure;imshow(real(z),[0 1],'InitialMagnification','fit')

Теперь исходное и восстановленное изображение отличаются только очень маленькими значениями, порядка eps, из-за обычной плавающей точки неточности: 

>> f-z
ans =
   1.0e-15 *
  Columns 1 through 2
  0.180411241501588 + 0.666133814775094i -0.111022302462516 - 0.027755575615629i
  0.000000000000000 + 0.027755575615629i  0.277555756156289 + 0.212603775716506i
  0.000000000000000 - 0.194289029309402i  0.000000000000000 + 0.027755575615629i
  Column 3
 -0.194289029309402 - 0.027755575615629i
 -0.222044604925031 - 0.055511151231258i
  0.111022302462516 - 0.111022302462516i
...