DFT and convolution, and fftshift, power spectrum

我们都知道傅里叶变换都有相应的卷积定理，对于限长离散信号的 DFT，有相应的 circular convolution, 但是DFT也可以用来求正常卷积。

本文要点都体现在最后的Matlab 例子。

1. g=ifft2(fft2(f).*fft2(h))
这样g 对应f和h的 circular convolution,
考虑一维信号，循环卷积通常是对两个长度为 N的信号说的，普通卷积有反折，平移，乘加三个步骤，循环卷积要把f或者h做周期为N的延拓，之后是同样的三个步骤。

Matlab 中 DFT 域循环卷积正确做法：(像素域可用 imfilter 完成,且h不进行如下的变换)
假设 f为mxn图像，h为size小于f的滤波器。f 不进行padding, 先将h移到mxn的中心，然后fftshift 移到四个角落。这样h会很怪，但是要想象把h进行mxn的周期延拓。
 g=ifft2(fft2(f).*fft2(h)) 结果为循环卷积
循环卷积对应矩阵相乘，矩阵为 circulant matrix.

2. Matlab 中 DFT域正常卷积正确做法：注：padding 为在尾部加zeroes


注意Matlab 中 fft2(f) 和 ifft2(F) 对应的 f 和 F 都是 centered at (0,0),亦即matlab的(1,1),因为 DFT 是 periodic，因此中心为 (0,0)(M-1,0)(0,N-1)...
如设计的频域滤波器中心在 M/2,N/2，要使用 fftshift 把中心移到4个角落

fftshift和ifftshift貌似作用是一样的，对于2D信号，将图像分为4部分，1、3平移互换，2、4平移互换。
也可以把信号和滤波器的中心都移到图像中心，但是反变换之前要移回来。

cropping: 假设filter 为 9x9，中心为(5,5)，用以上方法进行卷积计算后，最后步骤 cropping 要从 (5,5) 开始 crop，因为matlab 中图像中心为(1,1)，padding 不改变中心。

code:

clear 
clear 
f = imread('d:\Lib\DigitalImageProcessing\DIPUM\img\dipum_images_ch06\Fig0604(a)(iris).tif');
f = rgb2gray(f);
f = double(f)./255;
f = imresize(f,[128,128]);
[m n] = size(f);

%% Equivalence of padded DFT and normal convolution
%padding is appending zeros
addpath D:\Lib\DigitalImageProcessing\DIPUM\dipum_1.1.4

h2=zeros(9,9);
for i=-4:4
   for j=-4:4
      h2(i+5,j+5)= (1/(1+i*i+j*j));
   end
end
h2 = h2/sum(h2(:));

ps = paddedsize(size(f),size(h2));

out1 = real(ifft2(fft2(h2,ps(1),ps(2)).*fft2(f,ps(1),ps(2))));
out1 = out1(5:m+4,5:n+4);

figure
subplot 121
imshow(out1)

out2 = imfilter(f,h2,'same','conv'); %zero padding
subplot 122
imshow(out2)

% Test whether out1 and out2 are the same
find(abs(out1-out2)>10^-5)


%% Equivalence of non-padded DFT and circular convolution

middle = n/2 + 1;

h = zeros(size(f));
for i=-4:4
   for j=-4:4
      h(i+middle,j+middle)= (1/(1+i*i+j*j));
   end
end
h = fftshift(h); 
h = h/sum(h(:));

zz = real(ifft2(fft2(h).*fft2(f)));
zz2 = imfilter(f,h2,'same','conv','circular');

figure
subplot 121
imshow(zz,[])
subplot 122
imshow(zz2,[])

find(abs(zz-zz2)>10^-5)

FFTSHIFT 和 DFT/IDFT 的关系

注意箭头是双向的

Power Spectrum

对于实信号

注意是相关，而不是卷积。以上式子，对于普通相关和循环相关都是成立的。

code
注意此处 imfilter 取 full，然后从 (m,n) 处截取。这是因为 imfilter 使用 same 的截取方法不保证 (1,1) 处是自相关的原点。

clear
f = rand(4,4);
r = [1    -1     1    -1
    -1     1    -1     1
     1    -1     1    -1
    -1     1    -1     1];

% DFT
ffft = fft2(f);

ffft2 = fft2(f.*r);
find(ffft~=fftshift(ffft2))

ffft3 = fft2(fftshift(f));
find(ffft~=ffft3.*r)

% IDFT
finv = ifft2(f);

finv2 = ifft2(f.*r);
find(finv~=fftshift(finv2))

finv3 = ifft2(fftshift(f));
find(finv~=finv3.*r)

% Power Spectrum
zz = imfilter(f,f,'full','corr','circular');
[m n]=size(f);
zz=zz(m:end,n:end);

zz2 = ifft2(abs(fft2(f)).^2);

find(abs(zz-zz2)>10^-4)