Picking first arrival of a particular wave in a waveform

Question

Pritesh Bhoumick 2017-4-4

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/333476-picking-first-arrival-of-a-particular-wave-in-a-waveform

评论： Greg Dionne 2017-4-19

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 个评论
显示 2更早的评论隐藏 2更早的评论

Pritesh Bhoumick 2017-4-5

编辑：Pritesh Bhoumick 2017-4-5

在 MATLAB Online 中打开

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

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Greg Dionne 2017-4-6

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/333476-picking-first-arrival-of-a-particular-wave-in-a-waveform#answer_262097

在 MATLAB Online 中打开

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 个评论
显示 5更早的评论隐藏 5更早的评论

Greg Dionne 2017-4-19

I don't see your attached files(?)

Greg Dionne 2017-4-19

在 MATLAB Online 中打开

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

请先登录，再进行评论。

Picking first arrival of a particular wave in a waveform

4 个评论
显示 2更早的评论隐藏 2更早的评论

回答（1 个）

7 个评论
显示 5更早的评论隐藏 5更早的评论

另请参阅

类别

标签

Community Treasure Hunt

Picking first arrival of a particular wave in a waveform

4 个评论 显示 2更早的评论隐藏 2更早的评论

回答（1 个）

7 个评论 显示 5更早的评论隐藏 5更早的评论

另请参阅

类别

标签

Community Treasure Hunt

4 个评论
显示 2更早的评论隐藏 2更早的评论

7 个评论
显示 5更早的评论隐藏 5更早的评论