%% TRAITEMENT AUTO D'UNE SEQUENCE DE SISMO
clear all
directory='C:\Users\augereau\Documents\manu\ricochet\traitement\sismo\data_F3\seismoPorsmilin2502_03032017\*HHZ*.miniseed';
List=dir(directory)

for ii=1:2%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=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                   %
      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];                % paramètre du capteur : à modifier
                                      %à chercher sur le site https://www.passcal.nmt.edu/content/instrumentation/sensors/broadband-sensors/t120-bb-sensor
                                                %https://ds.iris.edu/NRL/sensors/nanometrics/nanometrics_trillium_sensors.html
      Ppoles=[                                  %instrument : trilium compact 120s
              -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];                %numériseur : centaur3X

    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]);

    % 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)));
    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(strcat('FFTBrest_',List(1).name(29:36),'.mat'),'FFTspectrum','FFTfrequencies','Time0')


%% 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
load CMAPBREST.mat
pcolor(Time0,FFTfrequencies{1},PSpectro')
set(gca,'Yscale','log')
shading interp
colorbar;
caxis([1E-13 5E-8]);
colormap (CMAP);
datetick('x','keepticks')
axis([Time0(1) Time0(end) 1E-2 50])
set(gca,'fontname','arial','fontsize',14)
xlabel('Date','fontname','arial','fontsize',14)
ylabel('Frequency','fontname','arial','fontsize',14)
