% eda10_08
% illustrate two-dimensional fourier transform
% Part 1, FFT of cosine wave, reordered so zero us in center
% Part 2, Illustration of use of symmetry

% standard distance/wavenumber set up
Nx=32;
Nkx=Nx/2+1;
Dx=1;
x=Dx*[0:Nx-1]';
xmax = Dx*(Nx-1);
kxmax=2*pi/(2*Dx);
Dkx=kxmax/(Nx/2);
kx=Dkx*[0:Nx/2,-Nx/2+1:-1]';

Ny=32;
Nky=Ny/2+1;
Dy=1;
y=Dy*[0:Ny-1]';
ymax = Dy*(Ny-1);
kymax=2*pi/(2*Dy);
Dky=kymax/(Ny/2);
ky=Dky*[0:Ny/2,-Ny/2+1:-1]';

% Part 1: cosine wave
theta = (pi/180)*30;
t1 = [cos(theta), sin(theta)]';;
kr1 = kxmax/4;
F = zeros(Nx,Ny);
for n = [1:Nx]
for m = [1:Ny]
    F(n,m) = cos( kr1*t1(1)*x(n) + kr1*t1(2)*y(m) );
end
end

% amplitude spectral density
Ftt = fft2(F);
Fasd = ((Dx^2)/xmax) * ((Dy^2)/ymax) * sqrt((abs(Ftt).^2));

% put into more reasonable order, with zero in the center
Fasd2=zeros(Nx,Ny);
Fasd2(Nx-Nkx+1:Nx,Ny-Nky+1:Ny) = Fasd(1:Nkx,1:Nky);
Fasd2(1:Nkx-2,Ny-Nky+1:Ny) = Fasd(Nkx+1:Nx,1:Nky);
Fasd2(1:Nkx-2,1:Nky-2) = Fasd(Nkx+1:Nx,Nky+1:Ny);
Fasd2(Nkx-1:Nx,1:Nky-2) = Fasd(1:Nkx,Nky+1:Ny);

eda_draw('  ', F, 'caption F',' ',' ',' ',  Fasd, 'caption Ftt',' ',' ',' ', Fasd2, 'caption Ftt' );

% Part 2: example of the use of symmetry

% new function F is noise, which has a complicated Fourier Transform
F=random('Normal',0,1,Nx,Ny);
Ftt = fft2(F);

% for real data, use symmetries to build right hand columns of
% 2D Fourier Transform from left hand columns
% First copy left Ny/2+1 columns; they are treated as known
Ftt2(:,1:Ny/2+1)=Ftt(:,1:Ny/2+1);
% then rebuild right from left columns
for m = [2:Ny/2]
    mp=Ny-m+2;
    Ftt2(1,mp) = conj(Ftt2(1,m));
    Ftt2(2:Nx,mp) = flipud(conj(Ftt2(2:Nx,m)));
end

% invert back to F
F2 = ifft2(Ftt2);

% compute error, as a check
Ett=sum(sum(abs(Ftt2-Ftt)));
E = sum(sum(abs(F2-F)));
disp('Symmetry test:');
disp(sprintf('    error in transform %f', Ett));
disp(sprintf('    error in function  %f', E));

