% A small Fourier transform demo.
%
% A few points to note:
% - Frequency is periodic over the interval that corresponds to the
%   length of the signal.
% - Use the function 'fftshift' only for plotting.  Do not compute your
%   Fourier transform as F = fftshift(fft(f)).  Instead, for the Fourier
%   transform, use F = fft(f); for the inverse Fourier transform, use
%   fhat = ifft(F).  The reconstructed signal, fhat, should differ from f
%   only up to round-off errors.
%
%   Note that as fftshift is *not* idempotent if f has an odd number of
%   elements, i.e., fftshift(fftshift(F)) <> F.  This means that if f has
%   an odd number of elements and if you let F = fftshift(fft(f)) then,
%   later on when you compute fhat = ifft(fftshift(F)), you will find
%   that fhat <> f.
%
% Detailed comments are given below in the code.
%
% Du Huynh, March 2007.
% The University of Western Australia
% School of Computer Science and Software Engineering

% randomly generate 4 to 10 key points of the signal
Nkeypts = round(rand(1,1)*(10-4) + 4);
xkeypts = 0:(Nkeypts-1);
fkeypts = randn(1,Nkeypts);

% apply cubic interpolation to get a detailed signal f(x)
x = 0:0.3:Nkeypts;
f = interp1(xkeypts,fkeypts,0:0.3:Nkeypts,'cubic');
% get rid of those undefined terms (i.e. NaN - not a number) as
% interpolation usually can't handle end points (unless we extrapolate).
f = f(f ~= NaN);
N = length(f);

x = 0:N-1;
F = fft(f);               % F is the Fourier transform of f
omega = 0 : -1 : -(N-1);  % range of frequency values

% note that the frequency is periodic over the interval of length N, where
% N is the length (i.e. the number of elements) of our signal.  Thus, if we
% have N=32 and
%   omega = [0,-1,-2, ..., -15,-16,-17,-18, ..., -28,-29,-30,-31] then it
% would be the same as defining
%   omega = [0,-1,-2, ..., -15,-16, 15, 14, ...,   4,  3,  2,  1]
% i.e. the frequency values in the second half of the omega array can be
% redefined by the modulo 32 operation.
%
% Note that the frequency values wrap around.  So F(omega) = F(omega+N).
index = round(N/2)+1;   % index of the starting position of the second
                        % half of the omega array
omega(index:end) = mod(omega(index:end), N);

% Note that the function fft of Matlab computes the Fourier transform as
% follows:
%   F(omega) = sum( f(x).*exp(-i*2*pi*omega*x/N) )
% where x is assumed to be [0:N-1] if N points are sampled from function f.


figure('Position',[10,10,600,700]);
subplot(4,2,1);
plot(x, f, 'r-');
grid on
axis tight;
ax = axis;
xlabel('x');
ylabel('f(x)');
title('Original signal, f');

subplot(4,2,2);
plot(fftshift(omega), abs(fftshift(F)), 'g-');
xlabel('\omega');
ylabel('|F(\omega)|');
title('Fourier spectrum of f');
grid on
axis tight;

fprintf('Number of key points generated = %d.\n', Nkeypts);
fprintf('Length of signal = %d.\n', N);
fprintf('Frequency values vary from %d to %d.\n\n', min(omega), max(omega));
fprintf('Please press <Enter> in the Matlab command window to see\n');
fprintf('the successive partial reconstruction.\n');

for u=0:max(abs(omega))
    % Fu is the Fourier transform of f at frequency value u
    index = find(omega == u | omega == -u);

    for i=1:length(index)
        % Fu is obtained by zeroing out all the frequency components
        % except for -u and u.
        Fu = zeros(1,N);
        Fu(index(i)) = F(index(i));
        % reconsfu is the partial reconstruction using only the frequency
        % component pair -u and u.
        reconsfu = ifft(Fu);

        subplot(4,2,3+(i-1)*2);
        plot(x, real(reconsfu), 'b-');
        axis tight;
        grid on
        xlabel('\omega');
        title(['Info of signal (real part) at freq value ', ...
            num2str(omega(index(i)))]);

        subplot(4,2,4+(i-1)*2);
        plot(x, imag(reconsfu), 'b-');
        axis tight;
        grid on
        title(['Info of signal (imag part) at freq value ', ...
            num2str(omega(index(i)))]);
    end

    % Fuptou is obtained by zeroing out all the frequency components
    % except for those in the range [-u,u].
    Fuptou = zeros(size(F));
    allindex = find(abs(omega) <= u);
    Fuptou(allindex) = F(allindex);

    % reconsf is the partial reconstruction from Fuptou.
    reconsf = ifft(Fuptou);

    subplot(4,2,7); hold off
    % due to round-off error, reconsf may turn out to be complex
    % with small imaginary components.  So to ensure that the program
    % does not crash at the 'plot' function below, only the real
    % components are used for plotting.
    plot(x, real(reconsf), 'b-');
    axis(ax);
    grid on
    xlabel('x');
    ylabel('reconstructed f(x)');
    title('Partial reconstruction so far');
    drawnow;
    pause
end

subplot(4,2,7);
title('Reconstruction using all freq components');