%
%  fft_tests -- djm, 27 june 2005
%

clear

disp('')
disp('1)  timing test')
disp('2)  sin(sin(x))')
disp('3)  sin(sin(x)) coefficient')
disp('4)  interpolation')
disp('5)  differentiation')
disp('6)  integration')
disp('')

choice = input('choice -> ');

switch choice
    case{1}
        figure(choice);  clf;  hold on
        M0 = 8;  M = 12;  Ntests = 50;
        data = zeros(M,1);
        for m = 1:M;
            vect = rand(1,2^(M0+m));
            t1 = cputime;
            for j = 1:Ntests;
                junk = fft(vect);
            end
            data(m) = (cputime-t1)/Ntests;
        end
        pts = 2.^(M0 + (1:M))';
        plot(M0+(1:M),data./pts./log2(pts),'x')
        title('\bf fft timing test, N = 2^M')
        ylabel('\bf cputime/(N log2(N))')
        xlabel('\bf M = log2(N)')
    case{2}
        figure(choice);  clf;  hold on
        N = 32;  dx = 2*pi/N;  x = -pi:dx:pi-dx;
        y = sin(sin(x));  yf = fft(y)/N;
        k1 = fftshift(-N/2:N/2-1);
        subplot(2,1,1)
        plot(x,y,'x');  axis([-pi pi -1 1])
        title('\bf bessel series demo')
        xlabel('\bf x-axis')
        ylabel('\bf y = sin(sin(x))')
        subplot(2,1,2);  hold on
        plot(k1,log(abs(imag(yf))),'ro');  axis([-N/2 N/2-1 -40 6])
        kb = -N/2+1:2:N/2-1;
        plot(kb,log(besselj(kb,1)),'kx')
        xlabel('\bf k-axis')
        ylabel('\bf log |imag(c_k)|')
    case{3}
        figure(choice);  clf;  hold on        
        for N = 8:2:24
            dx = 2*pi/N;  x = -pi:dx:pi-dx;
            y = sin(sin(x));  yf = fft(y)/N;
            plot(N,log(abs(imag(yf(4)) - besselj(3,1))),'x')
        end
        title('\bf spectral accuracy of sin(sin(x))-coefficient')
        xlabel('\bf N-axis')
        ylabel('\bf log |imag(c_3) - J_3(1)|')
    case{4}
        figure(choice);  clf;  hold on
        N = 32;  dx = 2*pi/N;  x = -pi:dx:pi-dx;
        y = sin(sin(x));  yf = fft(y);
        k1 = fftshift(-N/2:N/2-1);
        
        N1 = 16*N;  dx1 = 2*pi/N1;  x1 = -pi:dx:pi-dx1;
        y1 = sin(sin(x1));  plot(x1,y1,'r');

        plot(x,y,'k.');  axis([-pi pi -1 1])
        title('\bf interpolation demo')
        xlabel('\bf x-axis')
        ylabel('\bf y = sin(sin(x))')
        yi = ifft(yf.*exp(i*(dx/4)*k1));
        plot(x+(dx/4),real(yi),'kx');
    case{5}
        figure(choice);  clf;  hold on
        N = 32;  dx = 2*pi/N;  x = -pi:dx:pi-dx;
        y = sin(sin(x));  yf = fft(y);  yp = cos(sin(x)).*cos(x);
        k1 = fftshift(-N/2:N/2-1);
        plot(x,yp,'bo');  axis([-pi pi -1 1])
        title('\bf differentiation demo')
        xlabel('\bf x-axis')
        ylabel('\bf y = (d/dx) sin(sin(x))')
        yi = ifft(i*yf.*k1);
        plot(x,real(yi),'kx');
    case{6}
        figure(choice);  clf;  hold on
        N = 32;  dx = 2*pi/N;  x = -pi:dx:pi-dx;
        y = sin(sin(x));  yf = fft(y);
        k1 = fftshift(-N/2:N/2-1);  k1(1) = 1;
        axis([-pi pi -1 1])
        title('\bf integration demo')
        xlabel('\bf x-axis')
        ylabel('\bf y = (d/dx) sin(sin(x))')
        yi = ifft(-i*yf./k1);
        plot(x,real(yi),'kx');

end