% testxcov.m
% create a time series with 2Hz, lasting for 3sec. dt=0.001
close all
clear all
format compact
dt=0.001;
t=[0:dt:3];
f=2;
x=sin(2*pi*f*t);
figure(1)
plot(t,x)
xlabel('time')
ylabel('Amplitude')
title('sin(x) with 2Hz frequency')

%calculate the autocorrelation function of x
[c,tlag]=xcov(x,1500,'unbiased');
nc=length(c);
no=(nc-1)/2+1; %zero lag
tlag(no)
cn=c./c(no);
figure(2)
plot(tlag*dt,cn)
xlabel('time lag (sec)')
ylabel('auto-correlation r(tau)')
title('auto-correlation function for sin(x)')

% new signal with 2Hz and 5Hz embedded
f2=5;
x2=sin(2*pi*f2*t);
x3 = x+x2;
figure(3)
plot(t,x3)
xlabel('time')
ylabel('Amplitude')
title('2 periodic signals with 2Hz and 5Hz frequencies')

%calculate the autocorrelation function of x
[c3,tlag]=xcov(x3,1500,'unbiased');
nc=length(c3);
no=(nc-1)/2+1; %zero lag
tlag(no)
c3n=c3./c3(no)
figure(4)
plot(tlag*dt,c3n)
xlabel('time lag (sec)')
ylabel('auto-correlation r(tau)')
title('auto-correlation function for 2 periodic signals')

% create a red noise sequence
y =zeros(3001,1);
y(1,1)=0;
a=0.75;
for i =2: 3001
   y(i) = a*y(i-1) + sqrt(1-a^2)*randn(1,1);
end
figure(5)
plot(t,y)

%calculate the autocorrelation function of y
[cy,tlag]=xcov(y,1500,'unbiased');
nc=length(cy);
no=(nc-1)/2+1; %zero lag
tlag(no)
cyn=cy./cy(no)
figure(6)
plot(tlag*dt,cyn)
xlabel('time lag (sec)')
ylabel('auto-correlation r(tau)')
title('auto-correlation function for a red noise')

% x5 = y + x3
x5=y'+x3;

%calculate the autocorrelation function of x
[c5,tlag]=xcov(x5, 1500,'unbiased');
nc=length(c5);
no=(nc-1)/2+1; %zero lag
tlag(no)
c5n=c5./c5(no)
figure(7)
plot(tlag*dt,c5n)
xlabel('time lag (sec)')
ylabel('auto-correlation r(tau)')
title('auto-correlation function for a red noise and 2 periodicities')