How to find start and end points of signal in data file?

I have several hundred data files containing electric field over time measurements from a real system. The data has ~20 seconds of lead and trail time where the signal is off, so the data at the beginning and end of the file is just noise. I want to isolate the portion of the file containing the actual signal.

However, that 20 seconds is not exact so I can't just use that to index into the data file. I need some way to detect the start and stop of the signal. Visually, it's incredibly obvious, but I'm struggling to do it reliably programmatically. Any help is appreciated. I do have the signal processing toolbox.

The file I attached is for an applied signal at 1500 Hz, with a 10 kHz sample rate. I have files with applied signals all the way down to milliHertz, so it's a wide range of frequencies to deal with.

[Matt J]

Where would the signal start and end in your posted example?

[Daniel Flisek]

The region where the magnitude sharply increases, right in the middle of the data.

Answer 1

[dpb]

load sample_data
T=sample_data;                          % shorthand table name
T.Properties.VariableNames(1)={'Time'}; % ditto; remove the (sec)
T.t=T.Time-T.Time(1);                   % get the trace time stamp
dt=mean(diff(T.Time));                  % just checkin' the sample rate
T0=10;                                  % about how much pre-sample time is there
N=ceil(T0/dt);                          % bound around that time for smaller dataset
n1=floor(0.8*N); n2=height(T)-n1;
%fprintf('%10d%8d%8d%8d%8d\n',[N n1 n2 height(T) n2-n1+1])  % just checkin'....

    100001   80000  140101  220101   60102

plot(T.t,T.Data)
hold on
plot(T.t(n1:n2),T.Data(n1:n2),'r-')

S=T.Data(n1:n2);            % just grab the signal in range of interest
ipt=findchangepts(S,MaxNumChanges=2,Statistic='rms')  % see how findchangepts fares
figure
Z=zscore(S);                % standardize to get to uniform magnitude
M = movmad(Z,1/dt);         % do a 1-sec moving MAD
plot(T.t(n1:n2),Z)
hold on
plot(T.t(n1:n2),M,'r-','linewidth',3)

figure
findchangepts(M,MaxNumChanges=4,Statistic='linear')

The MAD seems to be pretty good indicator; some tuning with a broader sample of the data will probably be required.

[dpb]

load sample_data
T=sample_data; 
clear sample_data
Neven=round(height(T),-2);
R=range(reshape(zscore(T.Data(1:Neven)),[],10));
plot(R)
[mx,ix]=maxk(R,2);
hold on
plot(ix,mx,'xr','markersize',10)

An alternative more empirical approach could be something like the above -- look to bracket where the range of the data changes drastically and then screen in the ranges of the signal that bracket those locations. The above uses 10 groups; locating the two max values gives a group number that includes the ROI start/stop locations. One can then take those to sections and hone in on them.

i1=Neven/10*(ix(1)-1);
i2=Neven/10*(ix(2)+1);
figure
plot(T.Data(i1:i2))

The above could probably be turned into a recursive function where each pass takes the shorter vector from the above subsectioning. After a few passes, the two tails could be reduced significantly that would likely make the final selection much more robust and at a lot lower computational cost than more exotic filtering or processing.

% 2nd iteration...
Neven=round(i2-i1+1,-2);
R=range(reshape(zscore(T.Data(i1:i1+Neven-1)),[],10));
figure
plot(R)

Warning: Graphics acceleration hardware is unavailable. Graphics quality and performance might be diminished. See MATLAB System Requirements.

[Daniel Flisek]

This is excellent information, thank you for taking the time. It takes quite a while to run the findchangepts function on mine, but my computer is very slow to begin with. Is that normal? Is there a way to speed it up, considering I have so many files to process? I can probably downsample most of the low frequency runs, at least.

Edit: Also, I don't have the toolbox for zscore. Seems to work about the same without it. Why did you use it?

[dpb]

That it is such a long dataset without isolating at least a shorter section is the bottleneck and expected from findchangepts because it does an iterative calculation on the series a point at a time, splitting between two parts and computing the fitting error to minimize overall fit of the given decision model. That's a lot of work. That's the point of the first steps shown to be able to at least get rid of the bulk of the data known to be of no interest.

Effective generic ways to do that would greatly benefit in knowing something about just how variable the signals are over the population, not just a single case.

zscore is in the Statistics Toolbox; it simply subtracts the mean and divides by the standard deviation (the "z-statistic"). One ends up with a series with zero mean and standard deviation of 1.0 which puts all the data traces on the same scale for comparison of the change in amplitude instead of being dependent upon the actual magnitude. Makes the detection logic absolute magnitude independent and looking for specific pattern instead. It has the benefit that it doesn't change the relative shape, just scales everything to a common basis which should, as noted, make the decision logic in the end much easier to generalize.

A poor man's substitute could look something like--for this app, the second dim argument would likely not be used, treating everything by column.

function [z,mu,sigma] = myzscore(x,dim)
  if nargin==1
    mu = mean(x);
    sigma = std(x);
  else
    % User supplied dim
    mu = mean(x,dim);
    sigma = std(x,0,dim);
  end
  % Now standardize X to give Z, taking care when sigma is zero
  sigma0 = sigma + (sigma==0);
  z = (x - mu) ./ sigma0;
end

[Daniel Flisek]

Thanks a lot. I did something similar after looking up the zscore function.

This method works quite well for finding the signal bounds for most of the files in my dataset. Thanks again for your help!

[dpb]

Glad to help and was trying be as general as could without relying too much on the one specific data set. It's fundamentally the same idea as @Star Strider shows; if the time traces were symmetric (and without testing, they look as if would be reasonably close), his threshold could be directly converted into a z-score.

If there are some specific files that seem to be difficult, perhaps attaching one or two of those would lead to some other suggestions/ideas to make it more robust yet.

Answer 2

[Star Strider]

Another way would be to calculate a suitable threshold, find the index or time values at the points where the signal crosses them, then go from there. (This finds the start and stop index values, then the time values from them, you can work with only the time values if you want. That will likely be faster.) The comparisons used in the find calls are of course inequalities, since it is highly unlikely that the signal will actally equal the threshold values used in that calculation, so the derived time values are necessarily approximate.

Try something like this --

tic

load('sample_data.mat')

toc

Elapsed time is 0.049176 seconds.

format long
% sample_data
Fs = 10E+3;

t = sample_data.(1);
e = sample_data.(2);

[emin,emax] = bounds(e);
emdn = median(e)

emdn = 
    -3.310361842000000e+03

ethr_hi = (emax+emdn)/2;
ethr_lo = (emin+emdn)/2;

idx_hi = find(e >= ethr_hi);
[ehi_start,ehi_stop] = bounds(idx_hi);
t_hi_start = t(ehi_start);
t_hi_stop  = t(ehi_stop);

idx_lo = find(e <= ethr_lo);
[elo_start,elo_stop] = bounds(idx_hi);
t_lo_start = t(elo_start);
t_lo_stop  = t(elo_stop);

t_start = bounds([t_hi_start,t_lo_start])

t_start = 
     7.472572644899999e+04

t_stop = bounds([t_hi_stop,t_lo_stop])

t_stop = 
     7.472772864900000e+04

pulse_duration = t_stop - t_start

pulse_duration = 
   2.002200000002631

toc

Elapsed time is 0.154118 seconds.

figure
plot(t, e, DisplayName='Signal')
yline(emdn, '--', DisplayName='Median', LineWidth=1.5, Color=[1 0 1])
yline(ethr_hi, '--', DisplayName='High Threshold', LineWidth=1.5, Color=[1 1 1]*0.5)
yline(ethr_lo, '--', DisplayName='Low Threshold', LineWidth=1.5, Color=[1 1 1]*0.5)
xline(t_start, '-g', DisplayName='Start', LineWidth=1.5)
xline(t_stop, '-r', DisplayName='Stop', LineWidth=1.5)
grid
xlabel('Time (s)')
ylabel('E (units)')
legend(Location='best')
title('Original Data')

toc

Elapsed time is 0.536908 seconds.

figure
plot(t, e, DisplayName='Signal')
yline(emdn, '--', DisplayName='Median', LineWidth=1.5, Color=[1 0 1])
yline(ethr_hi, '--', DisplayName='High Threshold', LineWidth=1.5, Color=[1 1 1]*0.5)
yline(ethr_lo, '--', DisplayName='Low Threshold', LineWidth=1.5, Color=[1 1 1]*0.5)
xline(t_start, '-g', DisplayName='Start', LineWidth=1.5)
xline(t_stop, '-r', DisplayName='Stop', LineWidth=1.5)
grid
xlim([0.99999*t_start 1.00001*t_stop])
xlabel('Time (s)')
ylabel('E (units)')
legend(Location='eastoutside')
title('Zoomed Version')

toc

Elapsed time is 0.731735 seconds.

This may work generally for all your data, however since I don't have them I cannot determine that. It should be relatively straightforward to adapt it.

.

[Daniel Flisek]

I already accepted the other answer, but yours works just as well and is in fact faster. findchangepts() takes a long time on this dataset.

I did find a typo on this line:

[elo_start,elo_stop] = bounds(idx_hi);

It should use idx_lo instead of idx_hi.

Thanks a lot for taking the time to answer.

[Star Strider]

My pleasure!

Corrrect. I copied and pasted the previous line, however didn't change the argument.

My objective here was to create code that was as efficient as possible, since execution speed seemed to be important.