Positioning a cutted signal from a ECK Signal

Me and my Group need to filter a ECK Signal with Findpeaks, We need to then level the Signal so that the highest point(R Spike) are all on the same Position. Also we need to start and end the Signal at 0 so that we can compare the cutted Signals and draw a conclusion. We started and got some help from this link

https://de.mathworks.com/matlabcentral/answers/401320-cutting-up-the-signal-into-repeating-parts

now that we cutted the signal we wanted to allign the Signals bringing the R spike of every signal to the same spot.

Also we still need to make the signals start and end at 0 but that is not as important as alligning the signals.

This is the code we have till now (look code below)

We included the Text file with the signal and also our code

this is what we get atm. and we dont really now how to contiun and finish the task

as you can see the R spikes are not alligned at all. and we want to align them all so all the r spikes are at the same spot

Thanks for the Help.

clear all;
D = readmatrix("3_1_ecg60s.txt");
tv = D(:,1);

[pks, locs] = findpeaks(tv, 'MinPeakHeight',-0.95, 'MinPeakDistance',600);
locs1 = locs + 480;
for k1 = 1:numel(locs1)-1                          % Define MUAP Frames
    tvc{k1} = tv(locs1(k1)-1:locs1(k1+1)-1);
end

plot(tvc{k1})
%findpeaks(tvc{k1},2,'Annotate','extents','WidthReference','halfheight')
figure(2)
hold all
for k1 = 1:numel(tvc)

    plot(tvc{k1}-tvc{k1}(1)) 
    %axis([750 1500 -3.5 1]);          % Plot MUAP Frames
end
hold off
grid

[Ameer Hamza]

Can you explain, what do you mean by start and end at 0?

Also, when aligning peaks, what should be the location of all peaks. For example, if will align all peaks at x=400, then some of the signals will need to be truncated or moved to the negative x-axis.

[Samuel the Helpless]

Ths signal should start and end at 0. wich schould mean that the each signal needs to be reduced or raised by the first value and its difference to 0. (or atleast that how i think it should work) then each signal will be checked for its differnece to 0 and then the value will be added or subtracted for each value in the cutted signal so that the signal doesnt change only the values in the signal.

but if we align all the signals at x = 800 then we would move them all to the left wich would still make it able to compare the signals. ? or did i understand it wrong. ? the goal is to basicly move the signals all so that the R spike is for each signal at the same place.

Answer 1

[Star Strider]

You quoted part of my previous code (not entirely correctly) , so I’m ‘in for 1d, in for 1£ ', i guess.

Try this:

filename = '3_1_ecg60s.txt';
opts = detectImportOptions(filename);
D = readmatrix(filename,opts);
EKG = D(2:end,1);
Ts = 1E-3;                                                  % Sampling Interval (seconds)
tv = linspace(0, 1, numel(EKG))*Ts*numel(EKG);              % Time Vector
Start = 250;                                                % Ensemble Start: 250 ms Before R-Wave

figure
plot(tv, EKG)
grid
xlabel('Time (ms)')
ylkabel('Amplitude (mV)')
title('Original Signal')

[pks, locs, w, p] = findpeaks(EKG, 'MinPeakProminence',1.5);

for k1 = 1:numel(locs)-1                                   % Define EKG Frames
    EKGc{k1} = EKG(locs(k1)-Start:locs(k1+1)-Start);
    tvc{k1} = tv(locs(k1)-Start:locs(k1+1)-Start);
end

figure
hold on
for k1 = 1:numel(EKGc)
    plot(tvc{k1}, EKGc{k1})
end
hold off
grid
xlabel('Time (ms)')
ylkabel('Amplitude (mV)')
title('Segmented EKG')

figure
hold on
for k1 = 1:numel(EKGc)
    plot(tvc{k1}-tvc{k1}(1), EKGc{k1})
end
hold off
grid
xlabel('Time (ms)')
ylkabel('Amplitude (mV)')
title('Ensemble EKG')

It obviously still needs a good bit of filtering and other enhancements. I leave all that to you.

This should get you started.

[Samuel the Helpless]

Thx this is a great help. If a click the "answer accepted" button does it close the thred ?

there might be sth that would still come up.

but we can defenitly ty to expand on that.

Thank you!

[Star Strider]

My pleasure!

If you click the ‘answer acccepted’ it does not close the thread, and does mean that with my answer having been accepted, I will do my best to solve any other problems that may come up (if they have solutions, since not all do).

[Samuel the Helpless]

So we arent that good in matlab as we thought.

this is the signal that our professor used as an example how it should look at the end.

He even filtered the signal wich wasent in the task. (and is not needed)

What we wanted to do is take each cutted signal and start it at 0 ( on the y achses)

% the idea here was to take the first value in the vektor and calculate the difference to 0 
% that would be saved in another value and the difference would be added or subtracted to each value in the vektor(so the original signal stay but just starts at 0)
y = [-10:0]
x = [0:10]
c = 0

if y(first value of the vektor) = 10
    plot(y);
else 
    y ~= 10 
    for (y ~= 10)% this would check the value of y till its 0. then c(wich is the differnce to the original y and 0) we could use to applie to all the other vektors. 
        y + c 
        c = c +0.001
    end
end
% there must be someway to manipulate each value in a vektor but we dont really know it. 
%as you can see we hardly have matlab experience and so this should have a correct theroy but a totaly flawed programming
%we would applie this to all the Vektors we cutted so they would overlap on the same spot like the signal from my prof
% we dont really have any other clue how to do this.
% we tried this with a simple signal but failed. 
% also your code has 2 typos in the ylabel. where a k snuk itself in
%thank you very much ahead of time!

[Star Strider]

I am not certain what filter your professor used. I can get similar (although not exactly the same) results with an appropriate bandpass filter. Your professor also used a different starting point than I do. The 250 ms that I used is based on the fact that the normal P-R interval (beginning of the P-wave to the beginning of the QRS complex) is 100 to 120 ms (more than 120 ms is abnormal), and the normal QRS complex is less than 120 ms, although in a healthy heart, between 80 and 100 ms (it is consistently 100 ms in this trace). Taking half of that and adding 200 ms for the P-R interval gives 250 ms prior to the R-wave peak (that we are using to define each cycle here) for the beginning of each complete electrocardiographic cycle. The R-R intervals are not constant (Hering-Breuer reflex most likely), so the total lengths of each cycle are not the same.

If you want to start it with a y-value of 0, consider how you might do that for each cycle in the ensemble. (It is obvious, actually, although you will need to do the same relative calculation for each complete identified trace, in my code, each element of ‘EKGc’, the cell array of the electrocardiographic cycles.) When I change my code to do that, I get this result (including the bandpass filter):

As an aside, there is no obvious electrocardiographic pathology, so I would sign this out as: ‘Normal EKG’.

.

[Samuel the Helpless]

I dont know my partner thought we could use the same for funktion just for y values,

starty = 0
for k2 = 1:numel(pks)-1

 EKGc{k2} = EKG(pks(k2)-Starty:pks(k2+1)-Starty);
 tvc{k2} = tv(pks(k2)-Starty:pks(k2+1)-Starty); 
end

but we dont get it. we cant find the problem and i dont really know what the obvious answer should be. maybe ist just to much for our group. I though about manipulating the signal while we cut them into pieces. but other than that I have no clue.

[Star Strider]

It is not too much for your group. I am confident of that.

Consider the problem carefully.

Please note the difference (in my code) between ‘pks’ and ‘locs’.

The code in your latest post will throw an error because ‘pks’ (in my code) are not all integers (and I suspect none of them are, although I did not examine them closely). That aside, I suspect that using a ‘starty’ value of 0 in your code is going to start each vector of ‘EKGc’ in the centre of the R-wave (that assumes that you are subscripting ‘EKG’ correctly, and I believe you are not for the reason I already stated). The code you posted will throw the expected error. when you run it

If you want each record (element of ‘EKGc’) to begin at a y-value of 0, you will need to do the appropriate correction to each record individually. You can use the same loop to set the time start point and to correct the y-value, or you can correct the y-value in a separate loop. (I used two loops to create the image I posted most recently, because that is a new request, and I did not want to significantly alter my existing code. As you can see, it worked.)

.

[Samuel the Helpless]

figure
hold on
for k1 = 1:numel(EKGc) -1
   EKGcy{k1} = EKGc(locs(k1)-Starty:locs(k1+1)-Starty);
   
end
for k1 = 1:numel(EKGcy)
   
    plot(tvc{k1}-tvc{k1}(1), EKGcy{k1})
end

this is what we got. so we change EkGc to ekgcy to get the lines to 0 but how does he know how to change the value by to get to 0 ?

i just dont get it.

[Star Strider]

That code does not look correct to me, and I doubt if it is doing what you want it to do.

Think about it. If you want to make one element of a vector equal to 0 (or any number), and you do not want to change the additive relationships of that number with the other numbers in the vector, so that the vector ‘looks’ the same (when plotted) except for the shift, how would you do it?

[Samuel the Helpless]

yeah i would take the value of the difference to 0

if its 500 then it should add so long till it gets to 500 and then take the value and subtract or add it depending if its plus or minus

[Star Strider]

It is only necessary to do the addition or subtraction one time for each vector. (However, it will be necessary to do the same addition or subtraction on every element of the vector.)

[Samuel the Helpless]

Well then shouldnt it look similar to this

for k1 = 1:numel(EKGc)
    if EKGc <= 0
        for 
            EKGc = EKGc -0.0001 
        end
    else  EKGc >= 0
        EKGc = EKGc +0.001
    end
    plot(tvc{k1}-tvc{k1}(1), EKGc{k1})
end

i dont know any other way to programm it. and tbh if sb else wont show us the solution i dont think i can manage this task. and just quit....

[Star Strider]

You are seriously overthinking this!

Forgetting about the filter for the time being, arranging for all the elements of each element of cell array ‘EKGc’ to begin at 0 requires that you look specifically at the first element of every element of ‘EKGc’. Address the first element of each element of cell array ‘EKGc’ as:

FirstElement = EKGc{k1}(1);

since I suspect that one of the problems you are experiencing could be understanding cell array indexing conventions.

Note that if you want to address all the elements of a single element of cell array ‘EKGc’, that would be:

AllElements = EKGc{k1};

I know that ‘element’ here can be confusing. To clarify, each element of cell array ‘EKGc’ would be:

EKGc{k1}

and element ‘k2’ of each ‘EKGc{k1}’ cell array element would be:

EKGc{k1}(k2)

It is necessary to use cell arrays here because of the varying R-R intervals, and as the result, the varying times between consecutive R-waves.

[Samuel the Helpless]

well if i now that i still dont know how i should go about bringing the signal to 0.

c = EKGc{k1}(b)
j = 0
a = 0
if c == 0 % if c is zero the programm should check the next element of array. 
  b = 0 
  b = b+1
  a = a +1  % add 1 to check the next element and 1 for the other command. 
else if c > 0 % if c over 0 adds 0.001 till c = 0
    c = c +0.0001
else c < 0  % if c under 0 substracs 0.001 till c = 0
    c = c - 0.0001
end
for k1 = 1:numel(EKGc)  % checks the elements of EKCx 
    EKC{ka} + c  %first element of EKGc will be supptracted or added with c 
end
b = b+1 % add one for b to check the next EKGc element

i dont know this is really only way i can think of. its like you want to show me how to fish with a rod but the only way i know it is using my hands. thats how this feels to me. and it sucks cause my College should have showed me how to do this and they dindt or dont want to help.

[Star Strider]

The if blocks are not at all necessary.

Everything you need to do to get all of them to begin with a y-offset of 0 are in my previous Comment. I all but wrote the code, and since this is a homework assignment, we characteristically do not provide exact code unless there is no other way to answer the question. That is not the situation here.

[Samuel the Helpless]

OK thanks, ill try it tomorrow again. maybe our group is able to find the solution needed.

with a fresh mind.

[Star Strider]

My pleasure!

[Samuel the Helpless]

so we can solve the problem with all the excsting code we have? and we dont need to do code new methods or solutions. that we can take sth we got give it a little change and than we schould have our solution?

[Samuel the Helpless]

plot(EKGc{k1}-EKGc{k1}(1))

it seems that this piece of code changes the picture to this. wich is defenitfly better.

we just documentating our changes. so now we are just missing the end to get to 0.

[Samuel the Helpless]

so what the code basicly did was taking all elements of EKGc and subtract the first Element from all the elements. wich gives us 0 at the beginning.

plot(EKGc{k1}-EKGc{k1}(61))

i tried this to with the last element but that doesnt seem to work

[Star Strider]

That plot call is correct! Mine was slightly different:

plot(tvc{k1}-tvc{k1}(1), EKGc{k1}-EKGc{k1}(1))

There are several ways to remove the baseline wander. I chose to use a bandpass filter to get the image I posted in my earlier Comment. There are other options as well.

[Samuel the Helpless]

but that doesnt solve the endpoint that they all finish at 0. or does that have to do with the filter that all the signals end at 0 ? that if we filter the signal we get all the signals to end at 0 aswell

[Star Strider]

‘if we filter the signal we get all the signals to end at 0 aswell’

Yes. It may not be possible to entirely remove the trend, however it should be possible to remove almost all of it, as my plot image demonstrates.

[Samuel the Helpless]

so last week we desgined a filter and thought we could just copy it and integrate it

f = %signal we had
d = designfilt('lowpassfir', 'PassbandFrequency', 0.05, ...
               'StopbandFrequency', 150, 'PassbandRipple', 1, ...
               'StopbandAttenuation', 60, 'SampleRate', 10000);

fvtool(d);
dataOut = filter(d,f);
plot(dataOut)

this dindt really work as we expected and now we came up with 2 other filters wich dont reall work that great for us.

bandpass(EKG,[300 1300],Fs) %one time this 

b = fir1(48,[0.05 0.99],'bandpass'); %then this
EKG1 = filter(b,1,EKG);

so we just try and change things up to get a nice filtered signal.

we also tried to find help on the filter help page of matlab but that dindt really do the trick.

how do we decide what we need in a filter and wich values we should use ?

we filter this signal last weak.

e = readmatrix('3_2_ecg.txt')
f = e*5;


d = designfilt('lowpassfir', 'PassbandFrequency', 0.05, ...
               'StopbandFrequency', 150, 'PassbandRipple', 1, ...
               'StopbandAttenuation', 60, 'SampleRate', 10000);

fvtool(d);
dataOut = filter(d,f);
plot(dataOut);
hold on;
axis([0 10000 -700 -500]);
xlabel('Zeit in msek');
ylabel('Potential in mV');
xticks([961 2260]);
text(961, -680, 'Start');
text(2260, -680, 'Ende');
text(1455, -575, 'R');
text(1102, -620, 'P');
text(1320, -636, 'Q');
text(1543, -636, 'S');
text(2037, -610, 'T');
grid on;
title('Aufgaben-Skizze F8 & F9');
legend('Eingangssignal Ekg');
hold off;
saveas(gcf,'F7Skizze.jpeg');

to get this

wich is the most easiest and basic that we can use and how do we decide on the paramaters and wich funktion to use ?

[Star Strider]

A lowpass filter will only eliminate high-frequency noise. It will retain the baseline variations. The idea is to eliminate as much of the baseline variations and as much high-frequency noise as possible, without eliminating necessary frequencies that characterise the EKG signal itself.

Design your filter by doing a Fourier transform on the signal (to display its frequency content) to see how best to filter it. It is necessary to create an accurate frequency vector in order to correctly plot the one-sided Fourier transform of your signal, because you need to be able to determine what the passband and stopband frequencies should be.

With respect to filter design, I prefer elliptic filters, however there are several options.

I have contributed several posts on Fourier transforms and filter design for biomedical signal processing.

[Samuel the Helpless]

OK that seems like really much. and we will just say ok thats fine enough to work with and will use what we have. i bet there is a more easer solution that takes way less effort. cause we legit used 30 h for a 16 h and less task. well thanks for the help.

[Star Strider]

As always, my pleasure!

Doing the Fourier transform and then designing the filter is not difficult in MATLAB. This is my code for it:

Fs = 1/Ts;
Fn = Fs/2;
L = numel(EKG);
FT_EKG = fft(EKG - mean(EKG))/L;
Fv = linspace(0, 1, fix(L/2)+1)*Fn;
Iv = 1:numel(Fv);

figure
plot(Fv, abs(FT_EKG(Iv))*2)
grid
xlabel('Frequency (Hz)')
ylabel('Amplitude')
axis([0 100  0 0.05])

With the information from the Fourier transform, I designed this filter:

Wp = [ 5  40]/Fn;                                           % Normalised Passband (Passband = 5 Hz To 40 Hz)
Ws = [ 2  45]/Fn;                                           % Normalised Stopband (Passband = 2 Hz To 45 Hz)
Rp =  1;                                                    % Passband Ripple/Attenuation
Rs = 50;                                                    % Stopband Ripple/Attenuation
[n,Wp] = ellipord(Wp, Ws, Rp, Rs);                          % Calculate Elliptic Filter Optimum Order
[z,p,k] = ellip(n, Rp, Rs, Wp,'bandpass');                  % Elliptic Filter
[sos,g] = zp2sos(z,p,k);                                    % Second-Order-Section For Stability
EKG = filtfilt(sos, g, EKG);                                % Filtered EKG

I will post my entire code later if you want to see it.

.

[Samuel the Helpless]

yes thank you very much that would be great. thanks for all the help again.

[Star Strider]

This is simply for your information and reference. It would not be ethical to copy any of it and submit it as your own work.

My code for this problem (some of which I already posted):

filename = '3_1_ecg60s.txt';
opts = detectImportOptions(filename);
D = readmatrix(filename,opts);
EKG = D(2:end,1);
Ts = 1E-3;                                                  % Sampling Interval (seconds)
tv = linspace(0, 1, numel(EKG))*Ts*numel(EKG);              % Time Vector
Start = 250;                                                % Ensemble Start: 250 ms Before R-Wave

figure
plot(tv, EKG)
grid
xlabel('Time (ms)')
ylabel('Amplitude (mV)')
title('Original Signal')

Fs = 1/Ts;                                                  % Sampling Frequency (Hz)
Fn = Fs/2;                                                  % Nyquist Frequency (Hz)
L = numel(EKG);                                             % Signal Lengthj
FT_EKG = fft(EKG - mean(EKG))/L;                            % Fourier Transform
Fv = linspace(0, 1, fix(L/2)+1)*Fn;                         % Frequency Vector
Iv = 1:numel(Fv);                                           % Index Vector

figure
plot(Fv, abs(FT_EKG(Iv))*2)
grid
xlabel('Frequency (Hz)')
ylabel('Amplitude')
axis([0 100  0 0.05])

Wp = [ 5  40]/Fn;                                           % Normalised Passband (Passband = 5 Hz To 40 Hz)
Ws = [ 2  45]/Fn;                                           % Normalised Stopband (Passband = 2 Hz To 45 Hz)
Rp =  1;                                                    % Passband Ripple/Attenuation
Rs = 50;                                                    % Stopband Ripple/Attenuation
[n,Wp] = ellipord(Wp, Ws, Rp, Rs);                          % Calculate Elliptic Filter Optimum Order
[z,p,k] = ellip(n, Rp, Rs, Wp,'bandpass');                  % Elliptic Filter
[sos,g] = zp2sos(z,p,k);                                    % Second-Order-Section For Stability
EKG = filtfilt(sos, g, EKG);                                % Filtered EKG

figure
plot(tv, EKG)
grid
xlabel('Time (ms)')
ylabel('Amplitude (mV)')
title('Filtered Signal')

[pks, locs, w, p] = findpeaks(EKG, 'MinPeakProminence',1.5);

for k1 = 1:numel(locs)-1                                    % Define EKG Frames
    EKGc{k1} = EKG(locs(k1)-Start:locs(k1+1)-Start);
    tvc{k1} = tv(locs(k1)-Start:locs(k1+1)-Start);
end

figure
hold on
for k1 = 1:numel(EKGc)
    plot(tvc{k1}, EKGc{k1})
end
hold off
grid
xlabel('Time (ms)')
ylabel('Amplitude (mV)')
title('Segmented EKG')

figure
hold on
for k1 = 1:numel(EKGc)
    plot(tvc{k1}-tvc{k1}(1), EKGc{k1}-EKGc{k1}(1))
end
hold off
grid
xlabel('Time (ms)')
ylabel('Amplitude (mV)')
title('Ensemble EKG')

minLen = min(cellfun(@(x)numel(x), tvc));                   % Minimum EKG Segment Record Length
EKGce = cellfun(@(x)x(1:minLen), EKGc, 'Unif',0);           % Trim All Segments To Shortest Length
EKGensavg = mean(cell2mat(EKGce),2);                        % Calculate Ensemble AVerage

figure
plot((0:minLen-1)*Ts, EKGensavg)
grid
xlabel('Time (ms)')
ylabel('Amplitude (mV)')
title('Ensemble Averaged EKG')

You can also plot ‘EKGce’ if you want them all to be the same length. I calculated those and added the ensemble average to see what it would look like.

.