Télécharger (11 ko)

code de traitement de données sismo » Exemple_Process_trillium_BREST_data_F3_V1.m

Emmanuel Augereau, 16/10/2017 13:48

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







Télécharger le fichier comme

m