clear all
close all
%----------------------------------------------------------------------
% Load des données
[FileName,PathName]  = uigetfile(pwd,'sélection du fichier à traiter');
X=rdmseed([PathName,FileName]);
%----------------------------------------------------------------------
% Création des vecteur temps et signal
t = cat(1,X.t);				% En temps Matlab
d = double(cat(1,X.d));		% En count

Signal=double(d(1:2:end));
Tsignal=(t(1:2:end));

NfqSel=100;					%à faire correspondre à la fréquence d'acquisition
NyF=NfqSel/2;
%----------------------------------------------------------------------
% Deconvolution de la réponse intrsumentale
Sensi=4.344928536E+17*3.021E+8;  			% normalization facteur * sensitivité du numériseur (à trouver dans le firmware du numériseur)                      
CalibX(1:3) = [1 Sensi 11];					%11:nb de poles    
CalibX(33)=6;								%nb de zeros                 

% paramètre du capteur Zeros et Poles : à modifier
%à chercher les infos sur le site https://www.passcal.nmt.edu/content/instrumentation/sensors/broadband-sensors/t120-bb-sensor
% ou %https://ds.iris.edu/NRL/sensors/nanometrics/nanometrics_trillium_sensors.html
%instrument : trilium compact 120s
%numériseur : centaur3X

  ZZeros=[0.000000E+00  0.000000E+00
 0.000000E+00  0.000000E+00
-3.920000E+02  0.000000E+00
-1.960000E+03  0.000000E+00
-1.490000E+03  1.740000E+03
-1.490000E+03 -1.740000E+03];               

  Ppoles=[                                  
-3.691000E-02  3.702000E-02  
-3.691000E-02 -3.702000E-02  
-3.430000E+02  0.000000E+00  
-3.700000E+02  4.670000E+02  
-3.700000E+02 -4.670000E+02  
-8.360000E+02  1.522000E+03 
-8.360000E+02 -1.522000E+03  
-4.900000E+03  4.700000E+03  
-4.900000E+03 -4.700000E+03  
-6.900000E+03  0.000000E+00 
-1.500000E+04  0.000000E+00];                
       
CalibX(4:14) = Ppoles(:,1)+1i*Ppoles(:,2);  	%de 4 à 4+nbpoles-1
CalibX(34:39) = ZZeros(:,1)+1i*ZZeros(:,2); 	%de 34 à 34+nb zeros-1
% dcorr=corr_inst(Signal, CalibX',1/200,1,'v',[1/200 1/190 99 100]);

dcorr=corr_inst(Signal, CalibX',1/NfqSel,1,'v',[1/300 1/290 NfqSel/2-1 NfqSel]); % valeurs ok en 10-5 ou 10-6

%-----------------------------------------------------------------------
% affichage du signal déconvolué de la réponse instrumentale

figure;
subplot(2,1,1)
hold on;grid on;box on;
plot(Tsignal(1:end),Signal(1:end),'k')			%travail entre 1 et 4319879 (ou 1 et end) -> pour une journée entière
datetick('x','keepticks')
ylabel('Amplitude (m/s)','fontname','arial','fontsize',14)
xlabel('Time','fontname','arial','fontsize',14)
title('Raw signal','fontname','arial','fontsize',14)

subplot(2,1,2)
hold on;grid on;box on;
plot(Tsignal(1:end),dcorr(1:end),'k')
datetick('x','keepticks')
ylabel('Amplitude (m/s)','fontname','arial','fontsize',14)
xlabel('Time','fontname','arial','fontsize',14)
title('Deconvolved signal','fontname','arial','fontsize',14)
%-----------------------------------------------------------------------

% Filtres 

% figure;

% passe haut
% [b,a] = butter(4,5/NyF,'high'); %féquences à modifier en fonction de ce qu'on veut voire
% Sigfilt =filter(b,a,dcorr);
% 
% subplot(3,1,1)
% hold on;grid on;box on;
% plot(Tsignal(1:4319879),Sigfilt(1:4319879),'k')
% datetick('x','keepticks')
% ylabel('Amplitude (m/s)','fontname','arial','fontsize',14)
% xlabel('Time','fontname','arial','fontsize',14)
% title('Filtered signal [>5Hz]','fontname','arial','fontsize',14)

%passe bande
fb=0.008;							% basse fréquence de coupure pour le filtre passe bande
fh=49;								% haute fréquence de coupure pour le filtre passe bande
[b,a] = butter(4,[fb/NyF fh/NyF],'bandpass');
Sigfilt =filter(b,a,dcorr);

subplot(3,1,2)
hold on;grid on;box on;
plot(Tsignal(1:end),Sigfilt(1:end),'k')
datetick('x','keepticks')
ylabel('Amplitude (m/s)','fontname','arial','fontsize',14)
xlabel('Time','fontname','arial','fontsize',14)
title(strcat('Filtered signal [',num2str(fb),'Hz - ',num2str(fh),'Hz]'),'fontname','arial','fontsize',14)

% passe bas
% [b,a] = butter(4,1/(NyF),'low');
% Sigfilt =filter(b,a,dcorr);
% 
% subplot(3,1,3)
% hold on;grid on;box on;
% plot(Tsignal(1:4319879),Sigfilt(1:4319879),'k')
% datetick('x','keepticks')
% ylabel('Amplitude (m/s)','fontname','arial','fontsize',14)
% xlabel('Time','fontname','arial','fontsize',14)
% title('Filtered signal [<1Hz]','fontname','arial','fontsize',14)

%-------------------------------------------------------------------------------
% Spectres

fs=100;				% à modifier en fonction de la fréquence d'acquisition
n = 2^nextpow2(length(dcorr));
f = fs*(0:(n/2))/n;
FFTsig=abs(fft((dcorr),n))/n;
FFTsigFiltered=abs(fft((Sigfilt),n))/n;

figure;
subplot 211
loglog(f(1:end),(FFTsig(1:n/2+1)),'k')
hold on;grid on;box on;
loglog(f(1:end),filter(ones(100,1)./100,1,(FFTsig(1:n/2+1))),'r','linewidth',2)
set(gca,'fontname','arial','fontsize',14)
ylabel('Amplitude','fontname','arial','fontsize',14)
xlabel('Frequency','fontname','arial','fontsize',14)
title('FFT spectrum')

subplot 212
loglog(f(1:end),(FFTsigFiltered(1:n/2+1)),'k')
hold on;grid on;box on;
loglog(f(1:end),filter(ones(100,1)./100,1,(FFTsigFiltered(1:n/2+1))),'r','linewidth',2)
set(gca,'fontname','arial','fontsize',14)
ylabel('Amplitude','fontname','arial','fontsize',14)
xlabel('Frequency','fontname','arial','fontsize',14)
title(strcat('FFT spectrum sig. filtered ',num2str(fb),'Hz - ',num2str(fh),'Hz'))

%---------------------------------------------------------------------------
%% SPECTROGRAM HF 

% close all
% load COLORMAPWHITE.mat	

d=1;                			% début (en count) des données sur lesquelles va etre effectué le spectrogramme
f=length(Tsignal); 				% fin (en count) des données sur lesquelles va etre effectué le spectrogramme. indiquer "length (Tsignal)" pour calculer le spectrogramme sur l'intégralité des données
window=60000;      			 	% nombre d'échantillon dans la fenetre. indiquer 60000 pour une fenetre de 20 mn (20min*60sec*50Hz)

noverlap=0.9*window;			% recouverment entre 2 fenêtres de calcul
nfft=linspace(0,NfqSel/2,4*2560);

[b,a] = butter(4,[fb/NyF fh/NyF],'bandpass');
Sig =filter(b,a,dcorr(d:f));

[S,fr,t,p]=spectrogram(Sig,window,noverlap,nfft,NfqSel);
mycmap(1,:)=[0 0 1];
ENV=(abs(hilbert(Sig)))./max(abs(hilbert(Sig)));
ENVS=filter(ones(window,1)/window,1,ENV);%smoothing envelope by computing moving average

figure;
subplot 211
hold on;grid on;box on;
plot(Tsignal(d:f),Sig./max(abs(Sig)),'k')
plot(Tsignal(d:f),ENVS,'r','linewidth',2)
set(gca,'fontname','arial','fontsize',12)
datetick('x','HH:MM:SS','keepticks');
axis([Tsignal(d) Tsignal(f) -0.5 0.5])

subplot 212
imagesc(Tsignal(d:f),nfft,(20*log10(abs(S))));
axis xy
xlabel('Time','fontname','arial','fontsize',12)
ylabel('Frequency (Hz)','fontname','arial','fontsize',12)
datetick('x','HH:MM:SS','keepticks');
set(gca,'fontname','arial','fontsize',12)
axis([Tsignal(d) Tsignal(f) 0 50]) 			%modification des valeurs min et max des axes x et y
colormap('jet');
caxis([-115 -80]);							%modification des valeurs de la colorbar
colorbar;

%---------------------------------------------------------------------
%enregistrement des données

output = uigetdir(pwd,'sélectionner le fichier d enregistrement');
cd(output);
res=struct('d',{d},'f',{f},'x',{Tsignal},'nfft',{nfft},'S',{S});
save(strcat('spectro_',FileName(29:36),'.mat'),'-struct','res');

%--------------------------------------------------------------------
%load et affichage d'un spectro déja calculé

[SpectroName,SpectroPath]  = uigetfile(pwd,'sélection d un fichier.mat représentant un spectro');
spectro=load([SpectroPath,SpectroName]);

d=spectro.d;
f=spectro.f;
Tsignal=spectro.x;
nfft=spectro.nfft;
S=spectro.S;

imagesc(Tsignal(d:f),nfft,(20*log10(abs(S))));
axis xy
xlabel('Time','fontname','arial','fontsize',12)
ylabel('Frequency (Hz)','fontname','arial','fontsize',12)
datetick('x','HH:MM:SS','keepticks');
set(gca,'fontname','arial','fontsize',12)
axis([Tsignal(d) Tsignal(f) 0 50])
colormap('jet');
caxis([-115 -80]);
colorbar;

break

%-------------------------------------------------------------------
%% TRAITEMENT AUTO D'UNE SEQUENCE DE SISMO --- à finir de coder.
%
%clear all
%directory='C:\Users\augereau\Documents\manu\ricochet\traitement\sismo\SEED\*HHZ.mseed';
%List=dir(directory)
%
%for ii=1:length(List)
%  
%    % Load des données
%    X=rdmseed(List(ii).name);
%
%    % Création des vecteur temps et signal
%    t = cat(1,X.t);% En temps Matlab
%    d = double(cat(1,X.d));% En count
%
%    Signal=double(d(1:2:end));
%    Tsignal=(t(1:2:end));
%
%    NfqSel=100;
%
%    % Deconvolution de la réponse intrsumentale
%    Sensi=(536E+9)/800;
%    CalibX(1:3) = [1 Sensi 2]; 
%    CalibX(33)=2;
%      ZZeros=[+0.00000E+00  +0.00000E+00    
%              +0.00000E+00  +0.00000E+00]
%
%      Ppoles=[-0.1481  +0.1481  
%              -0.1481  -0.1481];
%
%    CalibX(4:5) = Ppoles(:,1)+1i*Ppoles(:,2);
%    CalibX(34:35) = ZZeros(:,1)+1i*ZZeros(:,2);
%    % dcorr=corr_inst(Signal, CalibX',1/200,1,'v',[1/200 1/190 99 100]);
%    dcorr=corr_inst(Signal, CalibX',1/NfqSel,1,'v',[1/300 1/290 NfqSel/2-1 NfqSel]);
%
%    % Spectres
%    fs=NfqSel;
%    n = 2^nextpow2(863870);
%    f = fs*(0:(n/2))/n;
%    FFTsig=abs(fft((dcorr),n))/n;  
%    FFTspectrum{ii}=filter(ones(100,1)./100,1,(FFTsig(1:n/2+1)));%sauvegarde de tts les spectres + lissage du spectre
%    FFTfrequencies{ii}=f(1:end);
%    Time0(ii)=t(1);
%    
%    display([num2str(ii),'/',num2str(length(List))])
%    clear d
%    clear dcorr
%    clear Signal
%    clear Tsignal
%    clear FFTsig
%    clear t
%
%end
%
%save('FFTBRestTest.mat','FFTspectrum','FFTfrequencies')
%
%%% PLOT
%
%CSCALE=jet(length(List));
%PSpectro=zeros(length(List),length(FFTfrequencies{1}));
%
%figure;
%subplot 211
%loglog(FFTfrequencies{1},FFTspectrum{1},'color',CSCALE(1,:),'linewidth',1.5)
%PSpectro(1,:)=FFTspectrum{1};
%hold on;grid on;box on;
%set(gca,'fontname','arial','fontsize',14)
%ylabel('Amplitude','fontname','arial','fontsize',14)
%xlabel('Frequency','fontname','arial','fontsize',14)
%
%for ll=2:length(FFTspectrum)
%    loglog(FFTfrequencies{ll},FFTspectrum{ll},'color',CSCALE(ll,:))
%    PSpectro(ll,:)=FFTspectrum{ll};
%
%end
%axis([1/200 50 1E-12 1E-5])
%
%subplot 212
%pcolor(Time0,FFTfrequencies{1},PSpectro')
%set(gca,'Yscale','log')
%shading interp
%colorbar;
%caxis([1E-10 5E-8]);
%datetick('x','keepticks')
%axis([Time0(1) Time0(end) 1E-2 10 ])
%set(gca,'fontname','arial','fontsize',14)
%xlabel('Date','fontname','arial','fontsize',14)
%ylabel('Frequency','fontname','arial','fontsize',14)
%
%
%------------------------------------------------------------------------------------------
%%% GEOPHONES ---- à finir de coder
%
%% Footsteps
%%filename='D:\Travail\Proc\EOST2015\Brest\geophonesPorsmilin\02\25\02\05\45_816126000_0001_0.asc'
%
%filename='D:\Travail\Proc\EOST2015\Brest\geophonesPorsmilin\02\27\10\25\46_515431000_0001_0.asc'
%
%TEST=dlmread(filename);
%plot(TEST(2:end,2))
%
%% Q: What is in those files? Triggering function? 3C?? (no sig....)







