cutting up the signal into repeating parts

I was wondering how can I cut up the signal from the data into repeating parts(each containing a single action potential) and then plot all these action potentials on a single graph, at the end shift all these action potentials in time so they are aligned along the x-axis (time ). data is attached

[per isakson]

"a single action potential" where does it start and end? Start: The signal is increasing and passes -40 or something of that kind?

Answer 1

[Star Strider]

I seem to have seen this waveform somewhere else recently!

Try this:

D = xlsread('Star11.xls');
tv = D(:,1);
ap = D(:,2);

[pks, locs] = findpeaks(ap, 'MinPeakHeight',100, 'MinPeakDist',20);     % Determine Peaks & Indices

figure(1)
plot(tv, ap)
hold on
plot(tv(locs), pks, '+r')
hold off
grid

for k1 = 1:numel(locs)-1
    apc{k1} = ap(locs(k1)-10:locs(k1+1)-10);                            % Define MUAP Frames
    tvc{k1} = tv(locs(k1)-10:locs(k1+1)-10);
end

figure(2)
hold all
for k1 = 1:numel(apc)
    plot(tvc{k1}-tvc{k1}(1), apc{k1})                                   % Plot MUAP Frames
end
hold off
grid

That should get you started.

Experiment to get the result you want.

EDIT — Added plot:

[Jan]

+1. This answer is compact and clear. It contains all details to understand the solution and no confusing additions. The code is clean: Spaces, indentation, efficient indexing, even the comments start in the same column for a good readability, no unnecessary clutter. A pre-allocation would be useful for education. The diagram clarifies the result on first view. An experienced Matlab programmer can understand the solution in 10 seconds, a beginner might need a minute. Therefore this is an efficient answer.

@sara bitarafan: The signal contains sine waves with a wavelength of about 0.1. Is this the nature of the signal or an artifact of too strong filtering?

[Star Strider]

@Jan — Thank you!

[Star Strider]

My code produces the plot I posted. The peaks are aligned correctly.

The current Question is a follow-on to a previous post in which I did all the necessary initial signal analysis and processing. This is the same signal that my earlier code produced, as the desired result. I recognised it immediately. My Answer to the previous Question was Accepted.

I measured the distance between the peaks using the peak distances returned by using 'MinPeakHeight' to isolate only the peaks above 100, and derived my 'MinPeakDist' values from that.

As a physician and biomedical engineer, I know where MUAPS begin in time. A brief calculation gave me the necessary index offsets that I then used in my code.

My code produces the correct result.

[John BG]

Hi Start Strider

this is John BG <mailto:[e-mail] [e-mail]>

ok, you say you measured the distance between peaks.

But the code you have supplied in this answer doesn't have anything that measures anything between peaks, does it?

In your answer you just read the data, and jump straight on findpeaks with the perfect squelch and perfect distance between peaks.

Everyone knows that findpeaks is a tricky function, powerful indeed, yet out of the blues as you have done in your answer, it's not immediate to get the right MinPeaksDist.

If it's the case that you did the fore-work to find the right parameters for findpeaks in another answer, you haven't supplied any link to such previous work.

Please have a look at my plotting of just 3 cycles: the peaks suffer low frequency drift.

You did not supply any plotting in your original answer.

Only after I supplied my answer with all relevant plots you copied the plotting centered on the peak locations, as it's evident by the

EDIT added plot

Please don't take my comment the wrong way. I am just highlighting the obvious.

Appreciating time and attention

John BG

<mailto:[e-mail] [e-mail]>

[John BG]

Ok, why don't we ask the opinion of Sara, the question originator should always have the last word, don't you agree?

As a power user with many more points than me do you have a way to ask Sara to have a look and decide?

[Star Strider]

My Answers stand on their own. I do not include irrelevant details. If any OP has questions about my code and how I designed it, I address them and provide appropriate explanations.

I do not argue with others about the correctness or legitimacy of their Answers. To do so is patently inappropriate in a forum such as MATLAB Answers. ‘Flame wars’ such as you engage in here, absolutely do not belong in this forum.

I assume Sara has access to my previous Answer for this problem, since the attached file to the original post are the data that my code for that Accepted Answer produced.

[Star Strider]

If Sara wants an ensemble average of the MUAPS in these data, I provide the code and plot here:

% — CALCULATE & PLOT ENSEMBLE AVERAGE —                                                                    
minlen = min(cellfun(@numel, apc));                                     % Minimum Length Of MUAP Records
ens = zeros(minlen, numel(apc));                                        % Preallocate
for k1 = 1:numel(apc)
    ens(:,k1) = apc{k1}(1:minlen);                                      % Trim MUAPs To Shortest Length
end
ensavg = mean(ens,2);                                                   % Calculate Ensemble Average
ci95 = 1.96*std(ens,[],2)/sqrt(numel(apc));                             % Calculate 95% Confidence Intervals
eatv = mean(diff(tv))*(0:minlen-1);                                     % Time Vector For Ensemble AVerage

figure(3)
plot(eatv, ensavg, '-r', 'LineWidth',1)
hold on
plot(eatv, ensavg+ci95, ':g', 'LineWidth',1.5)
plot(eatv, ensavg-ci95, ':g', 'LineWidth',1.5)
hold off
grid
legend('Ensemble Average', '95% Confidence Intervals')

[sara bitarafan]

thank you for your answer. it is very well explained and complete.

[Star Strider]

As always, my pleasure!

[Star Strider]

To be clear, my code identifies and isolates the action potentials beginning 100 ms prior to the identified peak of each. The first plot my code creates is an ensemble of the identified potentials, not an ‘eye diagram’, since this is physiology, not communications engineering. An ‘eye diagram’ is completely irrelevant here.

With respect to the ensemble average that I calculated from the ensemble data in the first plot and presented as an ensemble average with appropriate confidence intervals in the second plot, the 95% confidence intervals are correct for the standard normal distribution, and are calculated as:

CI95 = norminv(0.975, 0, 1)

The corresponding value for tinv with 16 degrees-of-freedom is 2.12, slightly more accurate for a small sample, but not different from the large sample estimate.

[sara bitarafan]

Thank you :)

[Star Strider]

As always, my pleasure!

[Jan]

Interested readers might proceed here:

• https://www.mathsisfun.com/data/confidence-interval.html
• https://en.wikipedia.org/wiki/Confidence_interval

Answer 2

[John BG]

Hi Sara Bitarafan

this is John BG <mailto:[e-mail] [e-mail]>

please have a look at the attached script. .

1.-

Acquiring signal:

clear all;clc;close all
A=xlsread('Star11.xls')
t=A(:,1)'
s=A(:,2)'

figure(1)
plot(t,s);grid on

Removing DC

min_s=min(s)
s=s+abs(min_s)

center signal, closest to dominant cycle to a sin cos signal

max_s=max(s)
s=s-max_s/2
hold on
plot(t,s);grid on;axis tight

For this sample there's not much DC, but it may the case that without removing DC the sought eye diagram is difficult to get.

2.

Calculating the dominant cycle with fft:

.

S=fft(s)
absS=abs(S);
figure(2);
plot(absS);grid on;

max_absS=max(absS);
n_maxS=find(absS==max_absS)

.

2nd value is just FFT mirror

n_maxS=n_maxS(1)

.

If n_maxS = max_amount_cycles this would be the highest discernible requency, with the FFT: it would be just 2 time samples per cycle.

max_amount_cycles=floor(numel(t)/2)

.

n_maxS: amount of cycles

.

hold on;plot(n_maxS,max_absS,'ro')

.

nT: amount of samples per dominant cycle:

nT=floor(numel(t)/n_maxS)

.

dominant cycle T

t1=t([1:nT])
T=t1(end)

amount of lost samples ignoring .2

mod(nT,n_maxS)  % samples lost, not a crisis

.

instead of a for loop:

% sc=zeros(n_maxS,T)
% for ,,

s2=s([1:1:end-(numel(t)-n_maxS*nT)])

.

But the dominant cycle is not that constant.

figure(3);
for k=1:1:3   % n_maxS
    plot(s2(k*[1:nT]))
    hold on;
end
% 003

.

For the same basic cycle t(nT) we see 1,2 and 3 peaks

This is caused by a 2nd tone almost half way the dominant tone in amplitude and really close to the dominant.

fft approach is more reliable when there's a clear frequency peak but no high tones anywhere, particular near the dominant.

Yet the variable nT is going to be useful in next point.

Also, assuming the single action potential refers to a single peak

3.-

When the dominant cycle is not that constant we can use squelch:

Set a threshold that includes all peaks above a given amplitude.

close all

s=A(:,2)'
th1=100
figure(1);plot(t,s);grid on

[pks,locs]=findpeaks(s,'MinPeakHeight',th1)
hold on
plot(t(locs),pks,'bd')

avoid false peaks imposing min distance between found peaks, that may be for instance: nT

.

[pks,locs]=findpeaks(s,'MinPeakHeight',th1,'MinPeakDistance',floor(nT/2))
hold on
plot(t(locs),pks,'rd')
% 004

.

So far what can we assert about this signal and the 'potential' events you have mentioned in the question?

How often do the events take place

mean(diff(locs))

.

How much jitter suffer such events understanding jitter ~ standard deviation of locs

var(diff(locs))^.5

4.-

Plotting the centered diagram with all events:

Let be dt the + - span left and right of each peak event

dt=floor(.5*mean(diff(locs)))

And the eye diagram:

sc2=zeros(numel(locs),2*dt+1)   % + - window span around each peak 
for k=1:1:numel(locs)
    if locs(k)<dt        % 1st cycle, with 1st peak closer than dt to beginning of signal.
       s0=s([1:locs(k)+dt])
       sc2(k,:)=[zeros(1,2*dt+1-numel(s0)) s0]
    end

      if locs(k)>numel(t)-dt  % last cycle with peak closer than dt to end of signal.
          sc2(k,:)=s([locs(k)-dt:end])
      end

      if locs(k)<numel(t)-dt && locs(k)>dt
          sc2(k,:)=s([locs(k)-dt:locs(k)+dt])
      end

end

% 4. and the eye diagram is:
figure(2);
for k=1:1:numel(locs)
    plot(sc2(k,:))
    hold on
end
grid on
% 005

.

.

Sara

if you find this answer useful would you please be so kind to consider marking my answer as Accepted Answer?

To any other reader, if you find this answer useful please consider clicking on the thumbs-up vote link

thanks in advance for time and attention

John BG

<mailto:[e-mail] [e-mail]>

[sara bitarafan]

thank you for your answer. it it very well explained and neat.

[John BG]

Thanks Sara