clc % emol
clearvars %emol
close all %emol

A=1; w=2*pi*1; 
fs=200;
T=100; % 100, 1000

  % generation of harmonic signal (emol)
t=0:1/fs:T-1/fs;
s=A*sin(w*t);

  % generation of broadband white noise signal, with
  % frequency band limited to half the sampling rate (emol)
randn('state',0);
n=randn(size(s));
  % conversion to a band-limited white noise signal, with
  % cut-off frequency equal to 20 Hz (emol)
fc=20;
    % ninth-order low-pass digital Butterworth filter (emol)
    % b: vector with coefs. of moving average part
    % a: vector with coefs. of auto-regressive part
[b,a] = butter(9,fc/(fs/2)); 
n=filtfilt(b,a,n);
  % noise power = twice the harmonic signal power (emol)
  % [Var(n) = 2 x Var(s)]; signal-to-noise ration (SNR) = -3dB
n=sqrt(2)*(std(s)/std(n))*n; 
  % addition of signal (emol)
y=s+n;

maxlags=4*fs; % maximum lag of 4s for autocorrelation function (emol)
[Ryy,tau]=xcorr(y,y,maxlags,'unbiased');
tau=tau/fs;

[Rss,tau2]=xcorr(s,s,maxlags,'unbiased'); %emol
tau2=tau2/fs; %emol

[Rnn,tau3]=xcorr(n,n,maxlags,'unbiased'); %emol
tau3=tau3/fs; %emol

figure(1)
subplot(2,2,1) %emol
plot(t(1:8*fs),y(1:8*fs),'b-',t(1:8*fs),s(1:8*fs),'g-')
xlabel('Time (s)')
ylabel('\its\rm(\itt\rm) and \ity\rm(\itt\rm)')
grid on %emol

subplot(2,2,2) %emol
plot(tau(maxlags+1:end),Ryy(maxlags+1:end),'r')
xlabel('Lag (\it\tau)'); ylabel('Autocorrelation of y(t)') %emol
axis([0 4 -1.5 1.5])
grid on %emol

subplot(2,2,3) %emol
plot(tau2(maxlags+1:end),Rss(maxlags+1:end),'r')
xlabel('Lag (\it\tau)'); ylabel('Autocorrelation of s(t)')
axis([0 4 -1.5 1.5])
grid on %emol

subplot(2,2,4) %emol
plot(tau3(maxlags+1:end),Rnn(maxlags+1:end),'r')
xlabel('Lag (\it\tau)'); ylabel('Autocorrelation of n(t)')
axis([0 4 -1.5 1.5])
grid on %emol

% Two time lengths are considered in line 7: T = 100 e T = 1000 s.
% The second time length supplies a better result because
% the signal of concern is not transient.