  function Show_FFT_Factors()
% The DFT matrix F_{n} is dense, but it can generally be factored into a product
% of very sparse matrices. For example
%      F_{64} = A_{6}*A_{5}*A_{4}*A_{3}*A_{2}*A_{1}
% This script displays the structure of these very sparse factor matrices.
% They have two nonzeros per row. Thus F_{64}x can be computed very fast
% via A_{6}(A_{5}(A_{4}(A_{3}(A_{2}(A_{1}x))))). Each A_{q}*vector is O(n)
% and we have log2(n) such products. That is where the Onlog(n)) comes
% from.

close all
t = 6;
n = 2^t;
Fn = exp(-2*pi*sqrt(-1)/n).^((0:n-1)'*(0:n-1));
spy(Fn)
title(sprintf('The %1d -by- %1d DFT Matrix',n,n),'fontsize',16)
shg
pause
for q=1:t
    L = 2^q;
    m = n/L;
    omega = exp(-2*pi*i/L);
    D = diag(omega.^(0:L/2-1));
    I = eye(L/2,L/2);
    Ak = kron(eye(m,m),[I D ; I -D]);
    fill([.5 n+.5 n+.5 .5],[.5 .5 n+.5 n+.5],'y')
    hold on
    spy(Ak)
    hold on
    for k=(0:L:n)+.5
        plot([k k],[.5 n+.5],'k')
        plot([.5 n+.5],[k k],'k')
    end
    hold off
    axis off
    title(['F_{64} = A_{6} A_{5} A_{4} A_{3} A_{2} A_{1}    ' sprintf('This is A_{%1d}',q)],'fontsize',16) 
    shg
    pause
end

