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