Skip to content

Commit

Permalink
Cleaning function and solving bug on Q variations
Browse files Browse the repository at this point in the history
  • Loading branch information
@LucaDeSiena
LucaDeSiena committed Oct 8, 2024
1 parent 0227ab4 commit c690212
Show file tree
Hide file tree
Showing 7 changed files with 35 additions and 74 deletions.
Binary file modified Murat_inputMSH.mlx
View file
Binary file not shown.
Binary file modified Murat_inputRomania.mlx
View file
Binary file not shown.
Binary file modified Murat_inputToba.mlx
View file
Binary file not shown.
14 changes: 7 additions & 7 deletions bin/Murat_inversion.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -136,16 +136,16 @@
% Q inversion
rcQ_k = ray_crosses_Q(:,k);
rtQ_k = retain_Q(:,k);
A_k = A_i(rtQ_k,rcQ_k);
Q_k = Qm(rtQ_k,k);
luntot_k = luntot(rtQ_k);
time0_k = time0(rtQ_k);
rapsp_k = rapsp(rtQ_k,k);
tCm = tCoda(rtQ_k,k);
A_k = A_i(rtQ_k & rtQc_k,rcQ_k);
Q_k = Qm(rtQ_k & rtQc_k,k);
luntot_k = luntot(rtQ_k & rtQc_k);
time0_k = time0(rtQ_k & rtQc_k);
rapsp_k = rapsp(rtQ_k & rtQc_k,k);
tCm = tCoda(rtQ_k & rtQc_k,k);

[d1, const_Qc_k, ~, ~] = Murat_lsqlinQmean(tCm,tWm,Q_k,...
cf_k,sped,luntot_k,time0_k,rapsp_k);
const_Qc(rtQ_k,k) = const_Qc_k;
const_Qc(rtQ_k & rtQc_k,k) = const_Qc_k;

lCurveQ_k = lCurveQ(k);
FName = ['L-curve_Q_' fcName '_Hz'];
Expand Down
16 changes: 8 additions & 8 deletions bin/Murat_plot.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -187,16 +187,16 @@
% Then it plots first the logarithm of the energy ratio versus travel
% time.
storeFolder = 'Tests/Q';
energyRatio_k = energyRatio(rtQk,k);
energyRatio_k = energyRatio(rtQk & rtQck,k);
residualQ_k = residualQ(k);
Edirect_k =...
energyRatio_k./codaNoiseRatio(rtQk,k);
A_k = A_i(rtQk,rcQk);
luntot_k = luntot(rtQk);
time0_k = time0(rtQk);
rapsp_k = energyRatio(rtQk,k);
tCm = tCoda(rtQk,k);
Q_k = Qm(rtQk,k);
energyRatio_k./codaNoiseRatio(rtQk & rtQck,k);
A_k = A_i(rtQk & rtQck,rcQk);
Q_k = Qm(rtQk & rtQck,k);
luntot_k = luntot(rtQk & rtQck);
time0_k = time0(rtQk & rtQck);
rapsp_k = energyRatio(rtQk & rtQck,k);
tCm = tCoda(rtQk & rtQck,k);

CN_title =...
['Coda Normalization check ' fcName ' Hz'];
Expand Down
76 changes: 18 additions & 58 deletions bin/Murat_retainPeakDelay.m
View file
Original file line number Diff line number Diff line change
@@ -1,49 +1,8 @@
% function [pab,lpdelta_i,retain_pd_i,ray_crosses_pd_i]...
% =...
% Murat_retainPeakDelay(t_phase,l10pd_i,yesPD_i,Apd_i)
% % function [pab,lpdelta_i,retain_pd_i,ray_crosses_pd_i]...
% % =...
% % Murat_retainPeakDelay(t_phase,l10pd_i,yesPD_i,Apd_i)
% %
% % CREATES all constraints for peak delay inversion
% %
% % Input parameters:
% % t_phase: time relative to envelope
% % l10pd_i: logarithm of peak delay
% % yesPD_i: sets waveform where this must be computed
% % Apd_i: peak delay matrix
% %
% % Output parameters:
% % pab: coefficients of the linear relationship
% % lpdelta_i: residuals
% % retain_pd_i: keeps tab on which waveforms are kept for imaging
% % ray_crosses_pd_i: keeps tab on which rays are kept for imaging

% dataL = length(t_phase);
% outlierspd = false(dataL,1);
% fitrobust_i =...
% fit(t_phase(yesPD_i),l10pd_i(yesPD_i),'poly1','Robust','on');
% pab = [fitrobust_i.p1 fitrobust_i.p2];

% l10pdt = polyval(pab,t_phase);
% lpdelta_i = l10pd_i-l10pdt;
% I = abs(lpdelta_i) >= 2*std(lpdelta_i,'omitnan');
% outliers =...
% excludedata(t_phase(yesPD_i),lpdelta_i(yesPD_i),'indices',I(yesPD_i));
% outlierspd(yesPD_i) = outliers;
% retain_pd_i = yesPD_i & ~outlierspd;
% Apd_retain_pd_i = Apd_i(retain_pd_i,:);
% s_pd = sum(Apd_retain_pd_i);
% ray_crosses_pd_i = s_pd~=0;

% end


function [pab, lpdelta_i, retain_pd_i, ray_crosses_pd_i] ...
= ...
Murat_retainPeakDelay(t_phase, l10pd_i, yesPD_i, Apd_i)
= ...
Murat_retainPeakDelay(t_phase, l10pd_i, yesPD_i, Apd_i, stDevPD)
% function [pab,lpdelta_i,retain_pd_i,ray_crosses_pd_i] = ...
% Murat_retainPeakDelay(t_phase,l10pd_i,yesPD_i,Apd_i)
% Murat_retainPeakDelay(t_phase,l10pd_i,yesPD_i,Apd_i, stDevPD)
%
% CREATES all constraints for peak delay inversion
%
Expand All @@ -52,6 +11,7 @@
% l10pd_i: logarithm of peak delay
% yesPD_i: sets waveform where this must be computed
% Apd_i: peak delay matrix
% stDevPD: standard deviation around trending line
%
% Output parameters:
% pab: coefficients of the linear relationship
Expand All @@ -60,37 +20,37 @@
% ray_crosses_pd_i: keeps tab on which rays are kept for imaging

% Get the length of the input data
dataL = length(t_phase);
dataL = length(t_phase);

% Initialize outlier flag array
outlierspd = false(dataL, 1);
outlierspd = false(dataL, 1);

% Perform robust linear fitting on the data
fitrobust_i = fit(t_phase(yesPD_i), l10pd_i(yesPD_i), 'poly1', 'Robust', 'on');
fitrobust_i = fit(t_phase(yesPD_i), l10pd_i(yesPD_i), 'poly1', 'Robust', 'on');
pab = [fitrobust_i.p1 fitrobust_i.p2];

% Calculate fitted values and residuals
l10pdt = polyval(pab, t_phase);
lpdelta_i = l10pd_i - l10pdt;
l10pdt = polyval(pab, t_phase);
lpdelta_i = l10pd_i - l10pdt;

% Standard Deviation (3-sigma) method for outlier detection
mu = mean(lpdelta_i(yesPD_i), 'omitnan');
sigma = std(lpdelta_i(yesPD_i), 'omitnan');
lower_bound = mu - 1 * sigma;
upper_bound = mu + 1 * sigma;
mu = mean(lpdelta_i(yesPD_i), 'omitnan');
sigma = std(lpdelta_i(yesPD_i), 'omitnan');
lower_bound = mu - stDevPD * sigma;
upper_bound = mu + stDevPD * sigma;

% Flag outliers
I = (lpdelta_i < lower_bound) | (lpdelta_i > upper_bound);
I = (lpdelta_i < lower_bound) | (lpdelta_i > upper_bound);

% Update outlier flags
outlierspd(yesPD_i) = I(yesPD_i);

% Determine which waveforms to retain
retain_pd_i = yesPD_i & ~outlierspd;
retain_pd_i = yesPD_i & ~outlierspd;

% Process retained waveform data
Apd_retain_pd_i = Apd_i(retain_pd_i, :);
s_pd = sum(Apd_retain_pd_i);
ray_crosses_pd_i = s_pd ~= 0;
Apd_retain_pd_i = Apd_i(retain_pd_i, :);
s_pd = sum(Apd_retain_pd_i);
ray_crosses_pd_i = s_pd ~= 0;

end
3 changes: 2 additions & 1 deletion bin/Murat_selection.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -10,6 +10,7 @@
miPD = Murat.input.minimumPeakDelay;
fT = Murat.input.fitTresholdLinear;
QcM = Murat.input.QcMeasurement;
stDevPD = Murat.input.stDevPD;

Apd_i = Murat.PD.inversionMatrixPeakDelay;
A_i = Murat.Q.inversionMatrixQ;
Expand Down Expand Up @@ -91,7 +92,7 @@
l10pd_i = l10pd(:,i);
[pab,lpdelta_i,retain_pd_i,ray_crosses_pd_i]...
=...
Murat_retainPeakDelay(t_phase,l10pd_i,yesPD_i,Apd_i);
Murat_retainPeakDelay(t_phase,l10pd_i,yesPD_i,Apd_i,stDevPD);

% Qc
Qm_i = Qm(:,i);
Expand Down

0 comments on commit c690212

Please sign in to comment.