Picking first arrival of a particular wave in a waveform

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.

[Greg Dionne]

No promises, but any chance you'd be willing to post a few waveforms for testing?

[Joseph Cheng]

With only a small sample size to look at, but could you take the 50~80msec window, threshold at ≈±0.05 then pick out "close" consecutive positive and negative peaks, else last peak detected (given example figure 1)?

[Pritesh Bhoumick]

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]

@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!

Answer 1

[Greg Dionne]

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

[Greg Dionne]

This seems to work okay-ish... not everything is caught though. (using 100 for order and 1.5 for ovlp).

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 .1]);
    subplot(2,1,2)
    plot([wave recwave]);
  end
  drawnow
end

[Pritesh Bhoumick]

This is working but with the findchangepts as well I am getting comparatively good result but not for all the waveforms (which is what I need to the maximum extent). I tried FFT as well, however, am not a signal processing guy and hence not able to do that effectively.

[Greg Dionne]

OK. I'll see if I can tweak this approach a little for your application. Maybe not FILLGAPS, but perhaps a combination of ARBURG and FILTER running in forward and reverse directions. If you have the Statistics/Machine Learning Toolbox, you may get SEGMENT to work for you.

[Pritesh Bhoumick]

I implemented FFT and modified the whole thing as in the attached files. Is there a way to optimize the whole thing further? Thanks for your inputs.

[Pritesh Bhoumick]

this code involves ARBURG a bit and smooth filter.

[Greg Dionne]

I don't see your attached files(?)

[Greg Dionne]

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