% file bonus20.m in http://www.cis.hut.fi/Opinnot/T-61.246/Laskarit/
%
% Demo why it is important to segment non-stationary signals
% in Fourier-analysis.

% Jukka Parviainen, 22.10.2004

L    = 85;
t    = [1 : L]/10000;   % time-axis
mt   = max(t);

x1   = cos(2*pi*500*t);   % 500 Hz cosine
x2   = cos(2*pi*1500*t);  % 1500 Hz cosine
x    = [x1 x2];           % concatenate two cosines in time
x1x2 = x1 + x2;           % sum of two cosine

figure(20); clf;

% time-domain plots, left column
subplot(4,2,1);           % subplot(#columns, #rows, which axis)
plot([1 : 2*L]/10000, x); % [x1 x2]
title('time-domain plots');
ylabel('[x1 x2]');
set(gca,'XLim',[0 2*mt]);
subplot(4,2,3);           % x1+x2
plot(t, x1x2);    
ylabel('x1 + x2');
set(gca,'XLim',[0 2*mt]);
subplot(4,2,5);           % 500 Hz cosine at 0 .. 0.01s
plot(t, x1);      
ylabel('x1, 500 Hz');
set(gca,'XLim',[0 2*mt]);
subplot(4,2,7);           % 1500 Hz cosine at 0.01 .. 0.02s
plot(t+max(t), x2);
set(gca,'XLim',[0  2*mt]);
ylabel('x2, 1500 Hz');
xlabel('time (s)');

% frequency-domain plots, right column
magxF    = abs(fft(x));     % [x1 x2]
magx1x2F = abs(fft(x1x2));  % x1+x2
magx1F   = abs(fft(x1));    % x1
magx2F   = abs(fft(x2));    % x2

ws = 10000*[0 : L-1]/100;    % frequency axis, 0 .. fs
wl = 10000*[0 : 2*L-1]/200; 

% spectra are 0 .. fs/2, if you want to get 0 .. fs
% comment the axis-lines with %
subplot(4,2,2);             % x1+x2
plot(wl, magxF);            
axis([0 5000 0 60]);
title('frequency-domain plots, fs = 10000 Hz');
subplot(4,2,4);             % [x1 x2]
plot(ws, magx1x2F);         
axis([0 5000 0 60]);
subplot(4,2,6);             % x1
plot(ws, magx1F);           
axis([0 5000 0 60]);
subplot(4,2,8);             % x2
plot(ws, magx2F);           
axis([0 5000 0 60]);
xlabel('frequency (Hz)');

% Questions:
%
% - see the spectra x1+x2 and [x1 x2],
%   is it possible to see if the signal
%   in time domain is x1+x2 or [x1 x2]?
% - what happens if you segment signal
%   [x1 x2] at some other point, e.g. at t = 0.008s?


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Demo showing: 
% "The longer the sequence, the better the spectral resolution."

figure(21); clf;

t  = [1 : 200]/10000;  
mt = max(t);

x  = cos(2*pi*495*t);    % cosine 500 Hz
x1 = x(1:16);             % very small part of it

magxF  = abs(fft(x)); 
w      = linspace(0, 10000, length(magxF));  
magx1F = abs(fft(x1)); 
w1     = linspace(0, 10000, length(magx1F));

subplot(2,2,1); 
plot(t*1000, x, '-o'); 
title('Cosine 500Hz, N=200');
set(gca,'XLim',[0 1000*mt]);
xlabel('time (ms)');

subplot(2,2,2); 
plot(w, magxF, '-o'); 
title('Spectrum 500Hz, DFT N=200');
axis([0 5000 0 130]);
grid on;
xlabel('Frequency (Hz)');

subplot(2,2,3); 
plot(t(1:16)*1000, x1, '-o');  
title('Cosine 500Hz, N=16');
set(gca,'XLim',[0 1000*mt]);
xlabel('time (ms)');

subplot(2,2,4);  
plot(w1, magx1F, '-o'); 
title('Spectrum 500Hz, DFT N=16');
xlabel('Frequency (Hz)');
axis([0 5000 0 8]);
grid on;
