%% x is a sequence of "time" values, e.g. in seconds
% The last sample is removed so that if repeated, there 
% will not be two times in a row with least-significant digit 0.
% As an exa5mple, [x, x] would be..
% [0, .1, .2, ...  99.8, 99.9, 0, 0.1, 0.2, 0.3, ..., 99.9]
x = 0:.1:100;
x = x(1:end-1);
%% Y is a pure sine wave (single frequency)
y = sin(x*pi/10);
plot(y,'.-')
%% frequency domain of sin-wave (y)
% A single spike is seen at a very low frequency.
% A second spike is seen at a very high frequency due to aliasing.
freq = fft(y);
stem(0:length(freq)-1,abs(freq))
%% Create square-wave signal
s = sign(y);
plot(x,s,'.-')
%% Create an arbitrary signal repeating 5 times

% Pick one of these or make your own. Must have 200 samples.
% s = [1:100, 100:-1:1]; % Triangle
% s = (1:200)/200 -0.5;    % Saw-tooth
% s = [ones(1,50), zeros(1,50), ones(1,50), zeros(1,50)]; % square wave shown spring 2016 (higher frequency, harder to see)
s = [ones(1,100), zeros(1,100)]; % better square wave (lower frequency, easier to see)
% s = [rand*ones(1,50), rand*ones(1,50), rand*ones(1,50), rand*ones(1,50)];

% Always include this to make the signal repeats
s = [s s s s s];

%% frequency domain of s (square wave) signal
sf = fft(s);
%% Show frequency spectrum
% We see sine-waves at frequencies 5, 15, 25, etc.
% These are multiples 1, 3, 5, 7, etc. of the base frequency (odd
% multiples)
stem(0:length(sf)-1,abs(sf))
xlim([0 30])
%% Pull out the most significant peaks.
dft = abs(sf(1:5:end))';
aft = angle(sf(1:5:end))';
dft = dft(1:end/2);
aft = aft(1:end/2);

% DC (0-frequency value) has to be handled as a special case
dft(1) = dft(1)/2;
%% Add the sine-waves one at a time to reconstruct
% the original signal
tot = zeros(size(y));
for i = 1:length(dft)
    y = dft(i)*cos((i-1)*x*pi/10+aft(i));    
    tot = tot+y;

    figure(1) % Time domain
    clf
    plot(x,y,'b-');
    hold on;
    plot(x,tot,'r--');
    plot(x,s*500,'k');
    hold off
    disp('Paused. press any key to continue');
    
    figure(2) % Frequency domain
    clf
    tf = fft(tot/500);
    range = 1:5:length(x);
    stem(0:length(range)-1,abs(sf(range)),'k');
    hold on
    stem(0:length(range)-1,abs(tf(range)),'r');
    stem(i-1,dft(i),'b');
    xlim([0, 20]);
    hold off
    pause
end
