clear;

N_signal=1024;
% generating two-sin signal
for (k=1:1:N_signal)
          garm(k)=1*sin(2*pi/20*(k-1))+2*sin(2*pi/50*(k-1));
    
end;

plot(garm);

% ARMA process generating
ar(1)=0.5*randn(1);
ar(2)=0.9*ar(1)+0.5*randn(1);

for (i=3:1:N_signal)
    ar(i)=0.9*ar(i-1)-0.7*ar(i-2)+0.5*randn(1);
end;
Noisestd=1;
arnoise=Noisestd*ar;
wnoise=Noisestd*randn(1,N_signal);

signal=garm;
noisy_signal=signal+wnoise;
plot(1:N_signal,signal,1:N_signal,noisy_signal,'yellow');



signal_c=signal-mean(signal,2);

[acf1,lags]=xcorr(signal_c,signal_c);


for(tau=1:1:N_signal)
 acf(tau)=0;
 for(j=1:1:N_signal-tau)
    acf(tau)=acf(tau)+signal_c(j)*signal_c(j+tau-1);
 end;
  acf(tau)= acf(tau)/N_signal;
end;
plot(acf);


plot(lags,acf1/N_signal,1:N_signal,acf)


for(tau=1:1:N_signal)
 acfn(tau)=0;
 for(j=1:1:N_signal-tau)
    acfn(tau)=acfn(tau)+arnoise(j)*arnoise(j+tau-1);
 end;
  acfn(tau)= acfn(tau)/N_signal;
end;
plot(acf);

plot(1:N_signal,acfn)



%lets form covariance matrices
% of size N

N=1000;

Qss=zeros(N,N);
Qww=zeros(N,N);
tic;
for(k=1:1:N)
  for(m=k:1:N)
     tau=m-k+1;
  
      Qss(m,k)=acf(tau);
       Qss(k,m)=acf(tau); 
       Qar(m,k)=acfn(tau);
        Qar(k,m)=acfn(tau); 
  end;   
  Qww(k,k)=Noisestd^2;
end;  
toc

Qyy=Qss+Qww;


Qinv=inv(Qyy);

lsc_signal=transpose(Qss*Qinv*transpose(noisy_signal(1:N)));

plot(1:N,signal(1:N),1:N,noisy_signal(1:N),'green',1:N,lsc_signal(1:N),'red')



