%http://ivoronline.com/Projects/Speech/Cepstrum/Cepstrum.pdf
%1.5.     MATLAB functions 
% 
%1.5.1.  Figure 1.1. 
% 
% 
%% MAIN SIGNAL PARAMETERS.----------------------------------------------------------------------------------- 
clear;   
T       = 0.01;         %Period of signal in seconds. 
f       = 1/T;          %Frequency in periods per second. 
w       = 2*pi/T;       %Frequency in radians per second. 
Ts      = (T/6/2);      %Sample period in seconds is 2 samples per smallest period. 
Ts_cont = (Ts)/100;     %To draw continuous signal we take 100 points between two samples. 
  
%% CONTINUOUS SIGNAL.---------------------------------------------------------------------------------------- 
  
% DEFINE. 
t_cont    = (0:Ts_cont:T);   %Time in seconds. Show 1 period of signal. 
x_cont_0  = 0.5; 
x_cont_1  = sin(2*w*t_cont); 
x_cont_2  = sin(4*w*t_cont); 
x_cont_3  = sin(6*w*t_cont); 
x_cont    = x_cont_0+x_cont_1+x_cont_2+x_cont_3; 
  
% DRAW. 
figure(1);  hold off;  
plot(t_cont, x_cont_0, 'r'); grid on; hold on;  
plot(t_cont, x_cont_1, 'g');  
plot(t_cont, x_cont_2, 'b');  
plot(t_cont, x_cont_3, 'm');  
plot(t_cont, x_cont,   'k','LineWidth',2);  
  
% TITLE & LABELS. 
h  = xlabel('t [s]');  set(h,'FontSize',15); 
h  = ylabel('x(t)');   set(h,'FontSize',15); 
h  = title('x(t)=0.5+sin(2\omegat)+sin(4\omegat)+sin(6\omegat) , Ts=(T/6)/2'); set(h,'FontSize',15); 
  
%% DISCRETE SIGNAL.------------------------------------------------------------------------------------------ 
  
% TAKE SAMPLES FROM CONTINUOUS SIGNAL. 
step = round(Ts/Ts_cont); 
start= 50;  %To avoid that samples of fastest signal are always 0. 
t    = t_cont  (start:step:length(t_cont)-1);    
x_1  = x_cont_1(start:step:length(x_cont_1)-1);    
x_2  = x_cont_2(start:step:length(x_cont_2)-1);    
x_3  = x_cont_3(start:step:length(x_cont_3)-1);    
x    = x_cont  (start:step:length(x_cont)-1);    
  
% DRAW. 
figure(1);   
plot(t, x,   'k.','MarkerSize',30); grid on;   hold on;  
plot(t, x_1, '.' ,'MarkerSize',15);  
plot(t, x_2, '.' ,'MarkerSize',15); 
plot(t, x_3, '.' ,'MarkerSize',15);  
plot(t, x_3, '.' ,'MarkerSize',15);  
  
%% DFT. ----------------------------------------------------------------------------------------------------- 
[a0,ak,bk]= DFT_Get_a0akbk_FromSamples(x);  
[A0,Ak]   = DFT_Get_A0Ak_From_a0akbk(a0,ak,bk); 
freq=(1:length(Ak))/Ts/length(x); 
figure(2); hold off; 
bar([0 freq],[A0 Ak],0.2);  
h = xlabel(strcat('Frequency [Hz] \rightarrow t [s]')); set(h,'FontSize',15); 
h = ylabel(strcat('Amplitude')); set(h,'FontSize',15); grid on; hold on; 
h = title('Components: Asin(\omegat+\phi) \rightarrow  0.5-0.5cos(\omegat)'); set(h,'FontSize',15); 
  
%% DEFINE CONTINUOUS SIGNAL.--------------------------------------------------------------------------------- 
T       = 200;            %Period of signal in seconds. 
f       = 1/T; 
w       = 2*pi/T;          %Frequency in radians per second. 
Ts_cont = T/100;           %Sampling period is 100 points in smallest period. 
t_cont  = (0:Ts_cont:3*T); %Time in seconds. Show 2 periods of signal. 
x_cont  = 0.5-0.5*cos(w*t_cont); 
figure(2);   
plot(100+t_cont, x_cont,'r'); grid on; hold on; 
x2   =Ak(1:length(Ak));         %So that we take full 3 periods of so-called signal. 
freq2=freq(1:length(freq)); 
plot(freq2,x2,'r.','MarkerSize',20); grid on; 
  
%% CEPSTRUM.-------------------------------------------------------------------------------------------------  
 
[ac0,ack,bck]= DFT_Get_a0akbk_FromSamples(x2);  
[Ac0,Ack]    = DFT_Get_A0Ak_From_a0akbk(ac0,ack,bck); 
freq=(1:length(Ack))/100/length(x2); 
figure(4); bar([0 freq],[Ac0 Ack],0.1,'LineWidth',2);  
h = xlabel(strcat('Frequency [Hz] \rightarrow Quefrency ')); set(h,'FontSize',15); 
h = ylabel(strcat('Amplitude')); grid on; hold on;           set(h,'FontSize',15); 
h = title('Components: Asin(\omegat+\phi)');                 set(h,'FontSize',15);  