%c'est la sesion do 0908 96 pour sortir l compoante inertielle
%auparavant j'ai fait filtsh1 pour prendre les donnees
filtsh1
whos
%Grand total is 768507 elements using 8459992 bytes
%je fait un peu de place
clear agsh1 ansh1 hesh1 vesh1 vnsh1

%filtrage a frequence inertielle effective
[b,a]=fir1(600,[1/18 1/17]/15,kaiser(601,0));
[h,w]=freqz(b,a,2048,30);
subplot(211);
plot(w,abs(h));axis([0 0.1 0 1.1])


 
vifin=filtfilt(b,a,vif(:,1));    %regardons d'abord le bin 1
subplot(212);plot(vifin)
subplot(212);plot(abs(vifin))
title('module de vif filtre frequence inertielle')
subplot(211);title('module de reponse en frequence du filtre fir')
xlabel('frequence en hr-1')
% je prepare les fonctions rotations plus et moins

eiwt=exp((sqrt(-1))*2*pi*(1:3284)/494);eiwtmoins=exp((sqrt(-1))*(-2)*pi*(1:3284)/494);%pour evaluer les modules aplus et amoins je dois faire
%une moyenne sur une periode,probleme des bords a resoudre
%plus tard
moymob=ones(1,494)/494;
%vifinplus=vifin.*eiwt'; erreur il faut un point
%vifinmoins=vifin.*eiwtmoins';
%aplus=filtfilt(moymob,1,vifinplus);
%amoins=filtfilt(moymob,1,vifinmoins);
figure(2)
%subplot(211);plot(abs(aplus))
%subplot(212);plot(abs(amoins))
 %title('module, mm/s, composante circulaire anticyclo feffectif')
% title('module, mm/s, composante circulaire cyclo feffectif')
%figure(3)%subplot(212);plot(angle(amoins))
%subplot(211);plot(angle(aplus))
%ici ce sont les graphes de la phase de vifin*eiwt
%title('phase vifinmoins (bas)et vifinplus(haut)')
%normalement donc je peux reconstruire l'amplitude,aplus et amoins

vifinplus=vifin.*eiwt.';
vifinmoins=vifin.*eiwtmoins.';
amoins=filtfilt(moymob,1,vifinmoins);
aplus=filtfilt(moymob,1,vifinplus);
subplot(212);plot(abs(amoins))
 title('module, mm/s, composante circulaire anticyclo feffectif')
subplot(211);
plot(abs(aplus))
 title('module, mm/s, composante circulaire cyclo feffectif')
vinertm=aplus.*eiwtmoins.';
vinertp=amoins.*eiwt.';
%ces deux derniers ordes sont les bons car plus referent en fait a moinsplot(angle(vinertm))
plot(angle(vinertp))%il faudrait peutre prendre abs(plus)*eiwtmoins' pour onde inertielle
plot(abs(vif))
plot(abs(vif(:,1)))
hold on
plot(abs(vif(:,1)-vinertm),'w')
plot(abs(vinertm),'r')
subplot(212)
plot(angle(vinertm),'r')
hold on
plot(angle(vif(:,1)-vinertm),'w')
title('phase des courants inertiels et des residus')
subplot(211)
title('modules du courant filtre 2 heures,du courant inertie et du residu')
subplot(212)
xlabel('temps en pas de deux minutes dynaproc shet1')
figure
plot(real(vif(:,1)))
subplot(211)
plot(real(vif(:,1)))
hold onplot(real(vinertm),'w')
title('composante Est de courant filtre2 heure et partie inertielle') 
plot(real(vif(:,1))-real(vinertm),'g')
subplot(212)
plot(imag(vif(:,1)))
hold on
plot(imag(vinertm),'w'
)plot(imag(vif(:,1))-imag(vinertm),'g')
axis([0 3500 -200 200]
title('composante Nord de courant filtre2 heure et partie inertielle')
%save etuinert 