Picking first arrival of a particular wave in a waveform

3 次查看(过去 30 天)
I am trying to pick the first arrival of a particular wave (Shear Wave) in a wavetrain recorded by the sensors. Normally, the first arrival is strong, but at times it can be marred by the presence of noise or other waveforms arriving (Compressional Wave) before the Shear Wave arrives. There is no particular frequency at which the waveform is recorded. The input frequency is 6 MHz and the output signal is shown in the figure below:
The red line represents the first arrival of the shear wave. The wavetrain that can be seen starting at around 40-45 msec (in the x-axis) is compressional wave.
I am using 'smooth' filter to smoothen the waveform and then using findchangepts to highest changes in the smoothened waveform in the window 55-80 msec (x-axis) to detect the first arrival. It works on most occasion but fails when therer is higher change in the slope of the smooth waveform in more than the actual first arrival. Like in the figure below:
The first arrival in this case too is actually the peak seen around 70 msec. Please suggest what would be the best way to detect the first arrival using Matlab.
My piece of code is like this:
file_abs = abs(file(:,i)); % Revert all negative amplitudes to positive values
file_smooth = smooth(file_abs,0.04,'rloess'); % Apply 'Smooth' filter to the waveforms
ipt = (1250+(findchangepts(file_smooth(1250:1800),'Statistic','rms','MaxNumChanges',2)))*interval; % take y-axis values corresponding to 50 - 75 msec time frame (x-axis)
if(size(ipt,1)==1)
if(ipt(1)>min_tt)
output(i,iter) = ipt(1);
end
end
if(size(ipt,1)==2)
if(ipt(2)<max_tt)
output(i,iter) = ipt(2);
end
end
I have some 3.5 million input files for every measurements I make. So it is impossible to mark the first arrivals manually. Need to automate it as much as possible.
  4 个评论
Pritesh Bhoumick
Pritesh Bhoumick 2017-4-5
编辑:Pritesh Bhoumick 2017-4-5
I have shared 7 files each of them have data points for 127 waveforms separated by a 6 lines header. You can read them easily using 'textread' command and check. Some of the waveforms are have anomaly in detecting the shear wave first arrival and some are easy to deal with. The files are shared in the following link:
https://sooners-my.sharepoint.com/personal/pritesh_ou_edu/_layouts/15/guestaccess.aspx?folderid=1a1df138438ed48bc93d395142d085307&authkey=ASIflWzaKjsEu4N2unlWKIM&expiration=2017-05-05T19%3a00%3a36.000Z
You can use the following code to read thru the files that I have written.
%%Define initial parameter values
filename = {' 1.txt';' 4.txt';' 7.txt'}; % File name to be read
dtp = 2500; % No. of waveform data points for each wave form
hdr = 6; % No. of header lines
interval = 0.00000004e6; % Time interval for each data points (in usec)
min_tt = 55; % Anticipated minimum shear wave first peak travel time
max_tt = 75; % Anticipated maximum shear wave first peak travel time
output = zeros(int16(size(textread(filename{1},'%s'),1)/(dtp+hdr)),size(filename,1));
%%Generate time matrix
time = zeros(dtp,1); % Time (in usec)
for i = 1:(dtp-1)
time(i+1)=time(i)+interval;
end
%%Read input files and detect arrivals
for iter 1:3:7
input = textread(filename{iter},'%s'); % Input file reading
len = size(input,1); % Length of input file
wfn = len/(dtp+hdr); % No. of waveforms
%%Divide the waveform points for each waveforms
file = zeros(dtp,wfn); % Waveforms divided into columns
for i = 1:len
if(rem(i,dtp+hdr)>hdr)
file(rem(i,dtp+hdr)-hdr,fix(i/(dtp+hdr))+1)=str2double(input(i));
end
if(rem(i,dtp+hdr)==0)
file(dtp,fix(i/(dtp+hdr)))=str2double(input(i));
end
end
Lemme know if more details would help. Thanks for your time!
Pritesh Bhoumick
Pritesh Bhoumick 2017-4-5
@Joseph - I cannot assume a threshold value because it would be different for every waveform. Some waveforms have really strong amplitudes while some are weak. And moreover, some have higher noise amplitude and some dont. If even a single noise peak has high amplitude as 0.05, it will pick that as the first arrival then. Would be very tricky to take that way!

请先登录,再进行评论。

回答(1 个)

Greg Dionne
Greg Dionne 2017-4-6
I can (sort-of) get a feel when things are changing via AR modeling. Basically a sliding window approach that attempts to compare the window against (bi-directional) prediction. Nothing fancy here, just using FILLGAPS. Here's an animation using an order of 100 taking 120 samples on either side of the window. Maybe if you play with the order/overlap you can get something working.
for i=1:size(file,2)
analyze(file(:,i),100,1.2);
end
function analyze(wave,order,ovlp)
range = 1200:10:2200;
errdata = zeros(numel(range),2);
for j=1:numel(range)
i = range(j);
killwave = wave;
killwave(i:i+order) = nan;
recwave = fillgaps(killwave,round(ovlp*order),order);
errdata(j,1) = i;
errdata(j,2) = rms(wave(i:i+order)-recwave(i:i+order))/rms(wave(i:i+order));
subplot(2,1,1)
plot(errdata(1:j,1)+order/2,errdata(1:j,2));
axis([1 2500 0 3]);
subplot(2,1,2)
plot([wave recwave]);
drawnow
end
end
  7 个评论
Greg Dionne
Greg Dionne 2017-4-19
OK. I tweaked it a bit so that it tries to predict in the forward and reverse directions and compare against what actually is there. If you see a larger error then the greater the chance the model changed at that point. Feel free to play with the number of training points and the model order... Please sanity-check it (play with the order numbers and the length of the training). Good luck!
for i=1:size(data,2)
analyze(file(:,i),150,150);
end
function analyze(wave,ntrain,order)
range = 1200:10:2200;
errdata = zeros(numel(range),1);
for j=1:numel(range)
i = range(j);
subsetleft = wave(i+(1:ntrain));
subsetright = wave(i+ntrain+(1:ntrain));
guessright = arpredict(subsetleft, ntrain, order);
guessleft = flipud(arpredict(flipud(subsetright), ntrain, order));
errdata(j) = sum(abs(subsetleft - guessleft)) + sum(abs((subsetright - guessright)));
end
subplot(2,1,1)
plot(wave);
subplot(2,1,2);
plot(range,errdata);
xlim([1 numel(wave)]);
drawnow
end
function y = arpredict(x, n, order)
a = arburg(x, min(length(x)-1,order));
if any(isnan(a))
y = zeros(n, 1);
else
zi = filtic(1, a, x(end:-1:1));
y = filter(1, a, zeros(n,1), zi);
end

请先登录,再进行评论。

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by