% Signal Comparison and Correlation:
%
% @ Kefu Xue, Ph.D., June 2006
%
% three signals
%
T=10^-3; % pick a sampling frequency (the smaller T the better approximation)
Ds=2; % Duration of the sinusoids
De=3; % Duration of the exponential
Ns=Ds/T; % total number of samples for sinusoids
ts=[0:Ns-1]*T;
Ne=De/T; % total number of samples for exponential
te=[0:Ne-1]*T;
xa=cos(pi*ts);
ya=sin(pi*ts);
za=exp(-4*te);
Exa=sum(xa.*xa)*T; % Energy of sin(pi*t)
Eya=sum(ya.^2)*T; % Energy of cos(pi*t)
Eza=sum(za.^2)*T; % Energy of exponential 
disp(['The energy for cos() is ' num2str(Exa)]);
disp(['The energy for sin() is ' num2str(Eya)]);
disp(['The energy for exponential is ' num2str(Eza)]);
figure;subplot(2,1,1);plot(ts,xa,'b',ts,ya,'r',te,za,'g');axis tight;
title('Three signals in the comparison');xlabel('time in second');
%
% Correlation calculation:
% 
gamaxy=sum(xa.*ya)*T; % correlation of two sinusoids
gamaxz=sum(xa.*za(1:Ns))*T; % correlation of exponential and cos()
% For point multiplication, the vectors must be in the same dimension.
gamayz=sum(ya.*za(1:Ns))*T; % correlation of exponential and sin()
disp(['Correlation between sinusoids gamaxy ' num2str(gamaxy)]);
disp(['Correlation between cos() and exponential gamaxz ' num2str(gamaxz)]);
disp(['Correlation between sin() and exponential gamayz ' num2str(gamayz)]);
%
% Correlation functions
%
fgamaxz=xcorr(xa,za,1)*T;
fgamayz=xcorr(ya,za,1)*T;
disp(['Correlation function between cos() and exponential gamaxz ' num2str(fgamaxz)]);
disp(['Correlation function between sin() and exponential gamayz ' num2str(fgamayz)]);
fgamayz=fliplr(conv(fliplr(ya), za));
zerolag=length(za)+1;
fgamaxz=fliplr(conv(fliplr(xa), za));
ku=[(-length(za)+1):length(ya)-1]*T;
subplot(2,1,2);plot(ku,fgamayz,'r',ku,fgamaxz,'b');axis tight;
title('Exp. correlation to sin(in red) and cos(in blue)');xlabel('time in second');
[maxv,maxi]=max(fgamayz);
kshift=(maxi-(length(za)+1)); % Sgama=-(length(za)-1)=-Ex since Sx=0
tshift=kshift*T;
disp(['For the best correlation to sin(pi*t), it needs ' num2str(tshift) ' pi phase shift.']);
% end of file