Two fast algorithms of the 2-D DFT are presented here: the tensor transform-based 2-D DFT and the paired transform-based 2-D DFT. The paired 2-D DFT is the best algorithm for calculating the 2-D DFT of images, and it can be used for calculating other unitary transforms as well, such as the 2-D Hadamard transform, for the NxN case, when N is a power of 2. The main idea of these algorithms is in the unique respresentation of the 2-D image not in the form of 2-D matrix, but as the set of 1-D image-signals which split the 2-D DFT of the image by separate 1-D DFTs of the image-signals. Each such 1-D time signal is generated by the corresponding frequency. Art / February 10, 2009.

The following two codes are added for the the inverse 2-D DFT, by method of tensor transform, or complex splitting-signals of the 2-D DFT of the image: ifft_LxLT.m and ifft_NxNT.m, for inverting the 2-D DFT of the image NxN, when N is a prime number and the power of two, respectively. Art / December 8, 2009.

Fast 2-D DFT (tensor algorithm)

2-D DFT for image LxL when L>1 is a prime code: tdft_LxL.m
2-D DFT for image NxN when N is a power of 2 code: tdft_NxN.m
calculate the splitting-signal of the image code: ft_pst.m
calculate the cyclic group T(p,s) code: cyclic_group.m
example of processing an image 256x256 tree_byArttensor2DFFT.ps

Inverse 2-D LxL-point DFT, when L>1 is a prime code: ifft_LxLT.m
Inverse 2-D NxN-point DFT, when N is a power of 2 code: ifft_NxNT.m

Tensor representation of the image is used in the fast 2-D DFT. This representation allows for presenting the image in form of splitting-signals of the same length. The representation of the image is in the frequency-and-time domain. The set of frequencies in such a domain is incomplete; but they are selected to achieve a covering of all set of frequencies. The main idea of the tensor representation is in splitting the 2-D DFT of the image into a set of separate 1-D DFT of signals which we call splitting-signals, or image-signals. The splitting-signal is generated by a frequency (p,s) and carries the spectral information of the image in the frequency-points of the cyclic group with generator (p,s). The cardinality of the set of such generators for the image of size NxN equals N+1, when N is a prime >1, and 3N/2, when N is a power of two. The construction of this set is simple in the general case of N, too.

Examples of the Fast 2-D DFT
--------------------------------------------------- 2-D DFT for the image 5x5:
% Call: tdft_5x5.m
% Demo code for calculating the 5x5-point DFT by
% the tensor transform: six 5-point DFTs are used
%
% Dr. Artyom M. Grigoryan, 11/12/1997

N=5;
f=[ 1 2 1 3 1
2 4 1 1 2
1 3 2 2 1
4 2 1 1 3
2 4 1 2 1];
F2D=fft2(f);
% illustration of the method for one generator:
p=2; s=1;
Image_signal=zeros(1,N);
Image_signal=ft_pst(f,p,s);
% 7 14 6 12 9
TDFFT_signal=fft(Image_signal);
% 48,-0.4549-1.2286i,-6.0451-8.6453i,-6.0451+8.6453i,-0.4549+1.2286i
T_ps=cyclic_group(p,s,N);
% 0 0
% 2 1
% 4 2
% 1 3
% 3 4
FT=zeros(N);
for k=1:N
p1=T_ps(k,1)+1;
s1=T_ps(k,2)+1;
FT(s1,p1)=TDFFT_signal(k);
end;
% The whole 2-D DFT by the tensor transform
FT=zeros(N);
s=1;
for p=0:N-1
Image_signal=zeros(1,N);
Image_signal=ft_pst(f,p,s);
TDFFT_signal=fft(Image_signal);
T_ps=cyclic_group(p,s,N)+1;
for k=1:N
p1=T_ps(k,1);
s1=T_ps(k,2);
FT(s1,p1)=TDFFT_signal(k);
end;
end;
p=1; s=0;
Image_signal=zeros(1,N);
Image_signal=ft_pst(f,p,s);
TDFFT_signal=fft(Image_signal);
T_ps=cyclic_group(p,s,N)+1;
for k=1:N
p1=T_ps(k,1);
s1=T_ps(k,2);
FT(s1,p1)=TDFFT_signal(k);
end;
% to check if the 2-D DFT has been calculated correctly:
Inv_image=real(ifft2(FT));
% Inv_image=Inv_image.*(Inv_image>=0);
Inv_image=round(Inv_image);
fprintf(' %3d %3d %3d %3d %3d \n',Inv_image');
1 2 1 3 1
2 4 1 1 2
1 3 2 2 1
4 2 1 1 3
2 4 1 2 1

--------------------------------------------------- 2-D DFT for the image 8x8:
% Call: tdft_8x8.m
% Demo code for calculating the 8x8-point DFT by
% the tensor transform: 12 8-point DFTs are used
%
% Dr. Artyom M. Grigoryan, 11/12/1997

N=8;
f=[ 1 2 1 3 1 2 1 3
2 4 1 1 2 3 2 1
1 3 2 2 1 1 3 2
1 2 1 1 3 1 2 1
2 4 1 2 1 2 3 2
3 2 2 3 2 1 1 3
1 2 3 4 2 1 3 2
2 3 1 2 1 3 1 2];
F2D=fft2(f);
% 2-D DFT by the tensor transform
FT=zeros(N);
s=1;
for p=0:N-1
Image_signal=zeros(1,N);
Image_signal=ft_pst(f,p,s);
TDFFT_signal=fft(Image_signal);
T_ps=cyclic_group(p,s,N)+1;
for k=1:N
p1=T_ps(k,1);
s1=T_ps(k,2);
FT(s1,p1)=TDFFT_signal(k);
end;
end;
p=1; L=2;
for s=0:L:N-1
Image_signal=zeros(1,N);
Image_signal=ft_pst(f,p,s);
TDFFT_signal=fft(Image_signal);
T_ps=cyclic_group(p,s,N)+1;
for k=1:N
p1=T_ps(k,1);
s1=T_ps(k,2);
FT(s1,p1)=TDFFT_signal(k);
end;
end;
% to check if the 2-D DFT has been calculated correctly:
Inv_image=real(ifft2(FT));
%Inv_image=Inv_image.*(Inv_image>=0);
Inv_image=round(Inv_image);
fprintf(' %2d %2d %2d %2d %2d %2d %2d %2d \n',Inv_image');
1 2 1 3 1 2 1 3
2 4 1 1 2 3 2 1
1 3 2 2 1 1 3 2
1 2 1 1 3 1 2 1
2 4 1 2 1 2 3 2
3 2 2 3 2 1 1 3
1 2 3 4 2 1 3 2
2 3 1 2 1 3 1 2
--------------------------------------------------------------------------------------------------------

Fast 2-D DFT (paired algorithm)

2-D DFT for image NxN when N is a power of 2 code: ptdft_NxN.m
calculate the set of frequency points T'(p,s) code: Tps_set.m
calculate the set of all 3N-2 generators J'(N,N) code: Jset.m

Paired-transform based fast 2-D DFT uses a minimal number of multiplications and split the 2-D DFT into disjoint subsets of frequency-points. The 2-D paired transform is a 2-D frequency plus 1-D time representation of the image. The transform is not separable. The image is represented as a set of separate signals of different lengths. There is no redundancy in calculation of the 2-D DFT when using the paired transform. Below is the code for computing the 2-D NxN-point DFT, when N is a power of 2. The code for the NxN case, where N is a power of a prime L>2 will be posted later. / Art

Demo example of the fast 2-D DFT of the image 256x256
------------------------------------------------------------------------------------
2-D DFT for the tree image of size 256x256:
% Call: ptdft_NxN.m
% Demo code for calculating the NxN-point DFT by the 2-D paired transform
% This algorithm uses (3N-2) 1-D DFT of orders N/2, N/4, N/8,...,2, and 1.
% The case: N is a power of 2.
%
% Dr. Artyom M. Grigoryan, 1997 - 02/10/2009
% -----------------------------------------------------------------------------------
% Tree image can be found in the library of the Signal and Image Processing Institute
% University of Southern California by address http://sipi.usc.edu/database/
N=256;
fid=fopen('tree.img','rb');
f=fread(fid,[N,N]); fclose(fid); clear fid; f=f';
figure;
subplot(1,3,1);
image(f); title('Tree image');
colormap(gray(200));
axis('image'); axis('off');
% ===================================================
% 2-D DFT by the paired transform
FT=zeros(N);
[TT2,N_all]=Jset(N);
Ntwo=bitshift(N,-1);
t=0:Ntwo-1;
W=exp((-j*pi/Ntwo)*t);
r=log2(N);
N1=3*Ntwo;
L=Ntwo;
g=1;
k_start=1;
k_end=N1;
for r1=1:r
for kk=k_start:k_end
p=TT2(kk,1);
s=TT2(kk,2);
split_signalP=zeros(1,L);
Image_signal=ft_pst(f,p,s);
split_signalP=Image_signal(1:g:Ntwo)-Image_signal(Ntwo+1:g:N);
modify_signal=split_signalP.*W(1:g:Ntwo);
FFT_signal=fft(modify_signal);
T_ps=Tps_set(p,s,N,g)+1;
for k=1:L
p1=T_ps(k,1);
s1=T_ps(k,2);
FT(s1,p1)=FFT_signal(k);
end;
end;
N1=bitshift(N1,-1);
L=bitshift(L,-1);
g=bitshift(g,1);
k_start=k_end+1;
k_end=k_end+N1;
% display the process of filling the 2-D DFT ---
pause(.1);
subplot(1,3,2);
image(abs(FT)/N);
axis('image'); axis('off');
% ----------------------------------------------
end;
% FT(1,1)=sum(sum(f));
FT(1,1)=Image_signal(1)+Image_signal(1+Ntwo);
% end of calculation of the 2-D DFT =============================
title('2-D DFT by PT');
% display the 2-D DFT
subplot(1,3,2);
image(abs(FT)/N*4);
axis('image'); axis('off');
% -------------------------------------------------------
% to check if the 2-D DFT has been calculated correctly:
Inv_image=real(ifft2(FT));
Inv_image=Inv_image.*(Inv_image>=0);
Inv_image=round(Inv_image);
% display the inverse 2-D DFT image:
subplot(1,3,3);
image(Inv_image); title('Inverse 2-D DFT');
axis('image'); axis('off');
% see result of all calculation here:
% print -dps tree_byArtpaired2DFFT.ps
% ---------------------------------------------------------------------------------------