|
%% 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)
|