%   Call: tdft_NxN.m.m
%   Demo code for calculating the NxN-point DFT by
%   the tensor transform: 3N/2 N-point DFT are used.
%   The case: N is a power of 2.
%
%   Dr. Artyom M. Grigoryan,         11/15/2001

       clear all;
        N=256;    
% The tree image can be taken from the library of the Signal and Image Processing
% Institute, University of Southern California, by address http://sipi.usc.edu/database/
        fid=fopen('tree.img','rb');           
        f=fread(fid,[N,N]); fclose(fid);  clear fid; f=f';           

        T_start=cputime; 
        % =============================================================        
        % 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;
        T_end=cputime; 
        fprintf(' cpu-working-time %6.4fs to perform 2-D DFT by TT \n',...
        T_end-T_start); 

        % 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 original image and inverse 2-D DFT:
        figure;
        colormap(gray(200));
        subplot(1,2,1); 
        imagesc(f);   
        axis('image'); axis('off'); 
        h_title1=title('tree image 256\times256');
        set(h_title1,'Color','k','FontName','Times','FontSize',9);             
        
        subplot(1,2,2); 
        image(abs(FT)/N*4);  
        axis('image'); axis('off'); 
        h_title=title('2-D DFT by TT');
        set(h_title,'Color','k','FontName','Times','FontSize',9);                 
        S_t=sprintf('calculation time is %6.4fs ',T_end-T_start);
        h_t=text(50,290,S_t);
        set(h_t,'Color','k','FontName','Times','FontSize',9);   
        
        subplot(1,2,1); 
        imagesc(Inv_image);
        axis('image'); axis('off');
        h_title2=title('Inverse 2-D DFT');
        set(h_title2,'Color','k','FontName','Times','FontSize',9);                 
        % =============================================================        
        % print -dps tree_byArttensor2DFFT.ps            
 % --------------------------------------------------------------------