clear;
pathin='D:\ВШЭ\filtr\Lab7'
filename=sprintf('%s/system_out.dat',pathin);
fin=fopen(filename,'rt');
A=fscanf(fin,'%f',[5 inf]);% A - array of data
fclose(fin);

N=size(A,2);
N_max=10000;
year=A(1,1:N_max);
phi1_in=A(2,1:N_max);
phi2_in=A(3,1:N_max);
m_1=A(4,1:N_max);
m_2=A(5,1:N_max);

dt=year(2)-year(1);

 % part 2

  
  cwt(phi1_in,years(dt),'amor');


inp=phi1_in+1j*phi2_in;
m=m_1+1j*m_2; 

[spectr, omega] = ampl_fft(inp,dt)
plot(year,m_1)
[spectr_m, omega] = ampl_fft(m,dt)
  plot(omega/2/pi,abs(spectr),omega/2/pi,abs(spectr_m))

Q = 100;
FC=365/433;

 [ out_signal ]=PMInversion(year,m,dt,FC,Q)

 noise=10*randn(2,N_max);

 
plot(year,phi1_in,year,real(out_signal) )


 m_noisy=m+(noise(1,:)+1j*noise(2,:));
 plot(year,real(m_noisy),year,real(m))
 
 [ out_noisy_signal ]=PMInversion(year, m_noisy,dt,FC,Q)

plot(year,real(out_noisy_signal),year,phi1_in )




f_c = 0;

%filter center frequency year^-1
 f_0=300%filter parameter year^-1
 
 inv=1
 outfilename=sprintf('CW_TRF_%3d_%2d_1.dat',round(1/f_c),round(1/f_0));
  
 
[ filtered_signal ]=ChandPantFreqFilter(year,m_noisy,f_0,f_c,dt,FC,Q,inv,outfilename)


plot(year, real(m_noisy),year, real(inp),'black',year,real(filtered_signal),'green')

%inp=interp_signal+1j*interp_signal;
%year=time_new;
[ spectr_in, freq_in] = spect_fftn(year,inp);
 [ spectr, freq] = spect_fftn(year, m);
 plot(freq_in,abs(spectr_in),freq,abs(spectr)/5);
 


  [ spectr_fltr, freq_fltr] = spect_fftn(year, filtered_signal );
  plot(freq_fltr,abs(spectr_fltr),freq,abs(spectr));
  
  
  noise=(randn(N_max,1)+i*randn(N_max,1))';
mn=m+noise  
  plot(year,real(noise),year,real(m));
  f_c=0
  f_om=10;
 inv=0
 outfilename=sprintf('CW_TRF_%3d_%2d_1.dat',round(1/f_c),round(1/f_om));

[ inverted_signal ]=ChandPantFreqFilter(year,mn,f_om,f_c,dt,FC,Q,inv,outfilename);

plot(year, real(mn),year, real(inp),year,real(inverted_signal))

