function [M] = fa_demo(signal, num_harmonics)
%function [out, dc_term] = fa_demo(signal, num_harmonics)
%
% FA_DEMO: Graphically performs Fourier series analysis on one period
%          of a signal, displaying results as it calculates.
%
% Input Parameters:
%    signal        -- One period of a periodic signal to be analyzed.
%                     Must be real.  For signals longer than 100, the
%                     animation becomes very slow (depends on the speed
%                     of the machine it's being run on).
%    num_harmonics -- The number of positive-frequency coefficients to
%                     calcuate.  This parameter must be less than or 
%                     equal to one-half the length of signal. (Default: 10)
%
% Output Parameters:
%    out     -- the coefficients resulting from the Fourier series analysis.
%               If the signal is real, out will be a vector containing the 
%               complex Fourier coefficients from 1 to num_harmonics.  In
%               this case, the negative coefficients are equal to the positive.
%   dc_term  -- the DC offset term (or "zeroth Fourier coefficient")
%

if ~exist('num_harmonics','var') | isempty(num_harmonics)
    num_harmonics = 10;
end

if size(signal,2) > size(signal,1) 
    signal = signal';
end

if num_harmonics > length(signal)/2
    error(['Cannot extract more than ' num2str(floor(length(signal)/2)) ' harmonics.']);
end

if ~isreal(signal)
    error('Signal must be real for the demo.  Use fourier_analysis.m instead.');
end

figure(1);
clf;

N = length(signal);
n = (0:(N-1));
out = zeros(1,num_harmonics);

dc_term = sum(signal)/N;

frame = 1;

for k = 1:num_harmonics;
    ce = exp(-j*k*n/N*2*pi)';
    corr = ce.*signal;
    corr_sum = cumsum(corr);
    out(k) = corr_sum(end)/N;
    if abs(out(k)) < 1e-14
        out(k) = 0;
    end
    ax2 = [-num_harmonics num_harmonics 0 max(abs(out))+.001];
    csi = corr_sum./(1:N)';
    s = signal; s(end+1:end+2) = [1 -1];
    sig_ax = [0 N-1 min(s) - abs(min(s))*.1 max(s) + abs(max(s))*.1]; 
    re_bar = [.5 2.5 min([0; real(csi)]) max([real(csi); abs(csi)]) ];
    im_bar = [.5 2.5 min([0; imag(csi)]) max([imag(csi); abs(csi)]) ]; 
  
        for i = 1:N
            subplot(3,4,[1 2 3]); 
            hold off;
            h = plot(n,signal,'b-');
            set(h,'LineWidth',2);
            hold on;
            h = plot(n,real(ce),'r-');
            set(h,'LineWidth',2);
            stem(0:i-1,real(corr(1:i)),'g');
            axis(sig_ax);
            ylabel('Real Part');
            
            subplot(3,4,[5 6 7]); 
            hold off;
            subplot(3,4,[5 6 7]);
            h = plot(n,signal,'b-');
            hold on;
            set(h,'LineWidth',2);
            h = plot(n,imag(ce),'r-');
            set(h,'LineWidth',2);
            stem(0:i-1,imag(corr(1:i)),'m');
            axis(sig_ax);
            ylabel('Imaginary Part');

            subplot(3,4,4);
            bar([real(csi(i)) 0 0 ; 0, abs(csi(i)), 0;],'stacked');
            axis(re_bar);
            xlabel(num2str([real(csi(i)) abs(csi(i))]));
            title(['Signal Correlation';
                   'Real/Imag      Abs']);
            
            subplot(3,4,8);
            bar([0 0 imag(csi(i)); 0, abs(csi(i)) 0],'stacked');
            colormap([0 1 0; 0 1 1; 1 0 1]);
            axis(re_bar);
            xlabel(num2str([imag(csi(i)) abs(csi(i))]));

            drawnow
        end
     
        dout = [dc_term out];
        subplot(3,4,[9 10]);
        stem(-k:k, abs([fliplr(out(1:k)) dc_term out(1:k)]));
        ylabel('|C_k|');
        xlabel('Coefficient Number');
        ax = axis;
        axis([-k-1 k+1 ax(3:4)]);
        subplot(3,4,[11 12]);
        stem(-k:k, angle([fliplr(out(1:k)) dc_term out(1:k)]));
        xlabel('Coefficient Number');
        ylabel('\angle{}C_k');
        ax = axis;
        axis([-k-1 k+1 ax(3:4)]);
        
        pause;
end

