Picking first arrival of a particular wave in a waveform
8 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
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:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/162559/image.jpeg)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/162560/image.jpeg)
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:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/162561/image.jpeg)
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 Kommentare
Antworten (1)
Greg Dionne
am 6 Apr. 2017
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 Kommentare
Greg Dionne
am 19 Apr. 2017
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
Siehe auch
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!