function w = fft_wavedec(im,wtype,J,parent)
%FFT_WAVEDEC: FFT- basova implementace wavelet transformace
%   w = fft_wavedec(im,wtype,J,parent) vypocte wavelet transformaci na
%   vstupni obraz a to Fourierovou metodou. Pouziva podminky okraju.
%   Velikost vstupniho signalu 2D muze byt libovolna
% 	
%   Vstup:
%   - im : vstupni signal o rozmerech [Nx,Ny].
%   - (volitelne) wtype : wavelet filtr, zakladni je 'sym8'.
%   - (volitelne) J = [Jx Jy] : urovne dekompozice. Zaklad je zvolen
%   maximalni pocet
%   - (volitelne) parent : pokud je rodic roven 1, pak jsou vlnky filtru
%   pouzivane k rozkladu. Zakladni hodnota je 0
% 	
%   Vystup:
%   - w : obraz stejne velikosti jako je "im" v obvyklem zastoupeni. "W" se
%   sklada ze ctyr kvadrantu [LowLow HighLow; LowHigh HighHigh], vlevo
%   nahoře je low-pass filtrovany obraz, na kterém dale může být použita.
% 	
%   dale fft_waverec, fft_wavefilters, wfilters, fft_gdc_filter.
% 
%   Author: Florian Luisier and Thierry Blu, March 2007
%   Biomedical Imaging Group, EPFL, Lausanne, Switzerland.
%  Toolbox stazen z http://bigwww.epfl.ch/
% 
%   Reference:
%   [1] T. Blu and M. Unser, "The fractional spline wavelet transform: 
%   definition and implementation," Proc. IEEE International Conference on 
%   Acoustics, Speech, and Signal Processing (ICASSP'2000), Istanbul,
%   Turkey, 5-9 June 2000, vol. I, pp. 512-515.
%   [2] F. Luisier, T. Blu, M. Unser, "A New SURE Approach to Image
%   Denoising: Interscale Orthonormal Wavelet Thresholding," 
%   IEEE Transactions on Image Processing, vol. 16, no. 3, pp. 593-606, 
%   March 2007.

[Nx,Ny] = size(im);
Jmax = aux_dyadic_max_scales([Nx,Ny]);
if(nargin<2)
    wtype = 'sym8';
end
if(nargin<3)
    J = Jmax;
end
if(nargin<4)
    parent = 0;
end

%% zjisteni maximalni poctu iteraci
J(1) = min(Jmax(1),J(1));
J(2) = min(Jmax(2),J(2));

%% Wavelet filter

if(J(1)>0)
    FaX = fft_wavefilters(Nx,wtype);
else
    FaX = ones(2,Nx);
end
if(J(2)>0)
    FaY = fft_wavefilters(Ny,wtype);
else
    FaY = ones(2,Ny);
end
if(parent)
    if(J(1)>0)
        Wx = fft_gdc_filter(FaX(2,:),wtype);
        Wx2 = conj([Wx Wx]);
        FaX(2,:) = FaX(1,:).*Wx2;
    end  
    if(J(2)>0)
        Wy = fft_gdc_filter(FaY(2,:),wtype);
        Wy2 = conj([Wy Wy]);
        FaY(2,:) = FaY(1,:).*Wy2; 
    end
end
Gx = conj(FaX(1,:));
Hx = conj(FaX(2,:));
Gy = conj(FaY(1,:));
Hy = conj(FaY(2,:));

%% provedeni 2D wavelet transformace
w = im;
L = sort(J);
for j=1:L(1)
	for i=1:Nx,
		X = fft(im(i,:));
		Y = Gy.*X;
		Z = Hy.*X;
		Y = 1/2*(Y(1:Ny/2)+Y(Ny/2+(1:Ny/2)));
		Z = 1/2*(Z(1:Ny/2)+Z(Ny/2+(1:Ny/2)));
		im(i,:) = [ifft(Y) ifft(Z)];
	end;
	for i=1:Ny,
		X = fft(im(:,i).');
		Y = Gx.*X;
		Z = Hx.*X;
		Y = 1/2*(Y(1:Nx/2)+Y(Nx/2+(1:Nx/2)));
		Z = 1/2*(Z(1:Nx/2)+Z(Nx/2+(1:Nx/2)));
		im(:,i) = [ifft(Y) ifft(Z)].';		
	end;
	w(1:Nx,1:Ny) = im;

	Nx = Nx/2;
    Ny = Ny/2;
    
	im = w(1:Nx,1:Ny);
	Gx = Gx(1:2:end);
	Hx = Hx(1:2:end);
	Gy = Gy(1:2:end);
	Hy = Hy(1:2:end);    
end

% Perform the 2D wavelet transform for the remaining iteration(s) (if any)
%--------------------------------------------------------------------------
for j=1:(L(2)-L(1))
    if(L(2)==J(1))
        for i=1:Ny,
            X = fft(im(:,i).');
            Y = Gx.*X;
            Z = Hx.*X;
            Y = 1/2*(Y(1:Nx/2)+Y(Nx/2+(1:Nx/2)));
            Z = 1/2*(Z(1:Nx/2)+Z(Nx/2+(1:Nx/2)));
            im(:,i) = [ifft(Y) ifft(Z)].';
        end;
        w(1:Nx,1:Ny) = im;
        Nx = Nx/2;
        Gx = Gx(1:2:end);
        Hx = Hx(1:2:end);       
    else
        for i=1:Nx,
            X = fft(im(i,:));
            Y = Gy.*X;
            Z = Hy.*X;
            Y = 1/2*(Y(1:Ny/2)+Y(Ny/2+(1:Ny/2)));
            Z = 1/2*(Z(1:Ny/2)+Z(Ny/2+(1:Ny/2)));
            im(i,:) = [ifft(Y) ifft(Z)];
        end;
        w(1:Nx,1:Ny) = im;
        Ny = Ny/2;        
        Gy = Gy(1:2:end);
        Hy = Hy(1:2:end);          
    end    
	im = w(1:Nx,1:Ny);   
end
w = real(w);