clear;

load('signal.mat');

N_signal = length(signal);
dt = 1/365;




time_old=0:1/12:(N_signal-1)*1/12

plot(time_old,signal)
time_new=time_old(1):1/365:time_old(N_signal)

interp_signal = interp1(time_old,signal,time_new);
plot(time_new, interp_signal);
N_new=size(time_new,2)


ar = zeros(1, N_new);
ar(1)=0.5*randn(1);
ar(2)=-0.8*ar(1)+0.5*randn(1);
for i=3:1:N_new
ar(i)=-0.8*ar(i-1)+0.1*ar(i-2)+0.5*randn(1);
end

ar = 20*ar;

figure;
hold on;
plot(ar);
title('Шум');
grid on;
hold off;


noisy_signal = interp_signal + ar;


figure;
hold on;
plot(time_new, noisy_signal);
title('Сигнал с шумом');
xlabel('year');
ylabel('amplitude');
grid on;
hold off;

syms alpha  beta  f_c;
% если хотите аналитзировать выражения символьно,
%то их надо определить
%затем придать значения параметрам
%и подставить их
%G=double(subs(G))
%F=double(subs(F))

Q = 100;
f_c = 365/433;
alpha = 2 * pi * f_c;
beta = pi * f_c / Q;

G = [[-beta, -alpha]; [alpha, -beta]];
F = -G;
C = [[1,0];[0,1]];

sys = ss(G, F, C, 0);
bode(sys)
%Initial conditions
x0 = [0 0];
input=vertcat(noisy_signal,noisy_signal);

figure;
lsim(sys, input, time_new, x0);


[m, tOut, x]=lsim(sys, input, time_new, x0);

 outfilename='system_out.dat';
 fout=fopen(outfilename, 'wt');
 for (j=1:1:N_new)
   fprintf(fout,'%10.6f ',time_new(j));
    fprintf(fout,'%10.6e ',(input(1,j))); 
   fprintf(fout,'%10.6e ', (input(2,j)));
   fprintf(fout,'%10.6e ', (m(j,1))); 
   fprintf(fout,'%10.6e ', (m(j,2)));
   fprintf(fout,'\n'); 
 end;
 fclose(fout);
 
 